Problem 226: A Scoop of Blancmange

Project Euler Problem 226 Solution

EulerSolve provides an optimized solution for Project Euler Problem 226, A Scoop of Blancmange, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The blancmange function, also called the Takagi function, is the continuous nowhere-differentiable curve $$B(x)=\sum_{n=0}^{\infty}\frac{\phi(2^n x)}{2^n},\qquad \phi(u)=\operatorname{dist}(u,\mathbb{Z}).$$ Problem 226 asks for the area common to the region under this curve and the circle centered at $\left(\tfrac14,\tfrac12\right)$ with radius $\tfrac14$. The implementations do not search for a symbolic antiderivative. Instead, they turn the geometry into a one-dimensional integral and make every sample value of the blancmange function exact by exploiting its behavior on dyadic rationals. Mathematical Approach The key observation is that Simpson's rule samples points of the form $x=k/2^m$, and on exactly those points the blancmange series collapses to a finite binary-digit formula. That is why the numerical integration is efficient and reliable even though the curve itself is fractal. The Blancmange Function on a Dyadic Grid If $x=k/2^m$, then $2^n x$ is an integer for every $n \ge m$, so $\phi(2^n x)=0$ from that point onward. Therefore $$B\left(\frac{k}{2^m}\right)=\sum_{n=0}^{m-1}\frac{\phi\!\left(2^n\frac{k}{2^m}\right)}{2^n},$$ and the infinite series becomes a finite sum of $m$ terms. The implementations go one step further and avoid even that $m$-term evaluation by converting the value into a formula involving binary digit sums....

Detailed mathematical approach

Problem Summary

For a real number $u$ let $s(u)$ be the distance from $u$ to the nearest integer. The blancmange curve is the graph of the function

$$T(x)=\sum_{n=0}^{\infty}\frac{s(2^nx)}{2^n},\qquad 0\le x\le1,$$

and the region between the curve and the $x$-axis has area $\tfrac12$. Let $C$ be the circle with centre $\left(\tfrac14,\tfrac12\right)$ and radius $\tfrac14$. We must find the area of the part of this region that lies inside $C$, rounded to eight decimal places in the form 0.abcdefgh. ^[1] The curve is also called the Takagi curve, after Teiji Takagi, who used it as an example of a continuous function without a derivative. ^[2]

The whole solution in one sentence: every vertical line through the wanted region cuts it in a segment from the lower arc of the circle up to the curve, so the area is the integral of the segment's length over $0\le x\le\tfrac12$, which the programs evaluate with Simpson's rule on finer and finer grids of points $k/2^m$, where counts of binary digits give the curve exactly.

Mathematical Approach

We evaluate the curve exactly at the points $k/2^m$, describe the region, and find out how accurately Simpson's rule integrates a curve without a derivative.

1. The curve at a few points

The term with index $n$ is a zigzag: on every interval $[j/2^n,(j+1)/2^n]$ it rises with slope $1$ up to the height $2^{-n-1}$ and falls back to $0$ with slope $-1$. A few values by hand:

$T(\tfrac12)=\tfrac12$, because only the term $n=0$ is not zero;
$T(\tfrac14)=\tfrac14+\tfrac12s(\tfrac12)=\tfrac12$;
$T(\tfrac38)=\tfrac38+\tfrac12s(\tfrac34)+\tfrac14s(\tfrac32)$${}=\tfrac38+\tfrac18+\tfrac18=\tfrac58$;
$T(\tfrac13)=\tfrac13+\tfrac16+\tfrac1{12}+\dots=\tfrac23$, because $2^n/3$ always lies at distance $\tfrac13$ from the nearest integer.

Since $0\le s\le\tfrac12$, the terms are at most $\tfrac12,\tfrac14,\tfrac18,\dots$, so $0\le T\le1$; the terms after the first $M$ add up to at most $2^{-M}$ for every $x$, so the continuous partial sums converge uniformly and $T$ is continuous. Since $s$ is unchanged by integer shifts and by a change of sign, $T(1-x)=T(x)$. The zigzag of index $n$ consists of $2^n$ triangles with base $2^{-n}$ and height $2^{-n-1}$, of total area $2^{-n-2}$, and these areas add up to the given $\tfrac12$.

2. Exact values at the points $k/2^m$

Write $w(r)$ for the number of ones in the binary representation of $r$, its Hamming weight or population count, ^[3] for example $w(6)=w(110_2)=2$, and let $P(k)=w(0)+w(1)$${}+\dots+w(k-1)$. ^[4] Then, for all integers $m\ge0$ and $0\le k\le2^m$,

$$T\Bigl(\frac{k}{2^m}\Bigr)=\frac{mk-2P(k)}{2^m}.$$

To prove it, cut $[0,1]$ into the $2^m$ cells $[k/2^m,(k+1)/2^m]$. Inside cell $k$, the binary expansion $0.d_1d_2\ldots$ of $x$ begins with the $m$ digits of $k$. The terms with $n\ge m$ vanish at both ends of the cell, since $2^nk/2^m$ is an integer. A term with $n\lt m$ has its corners at multiples of $2^{-n-1}$, which are ends of cells, so it is straight on the cell; and since the fractional part of $2^nx$ is $0.d_{n+1}d_{n+2}\ldots$, it rises with slope $1$ if $d_{n+1}=0$ and falls with slope $-1$ if $d_{n+1}=1$. The first $m$ terms together therefore have the slope $(m-w(k))-w(k)=m-2w(k)$, so $T$ increases by $(m-2w(k))/2^m$ across cell $k$, and adding these steps from $T(0)=0$ gives the formula.

For $m=3$ the cells $k=000_2,\dots,111_2$ have $w(k)=0,1,1,2,1,2,2,3$ and the slopes $3,1,1,-1,1,-1,-1,-3$. Adding the steps, each divided by $8$, gives the values $T(k/8)$ for $k=0,\dots,8$:

$$0,\ \tfrac38,\ \tfrac12,\ \tfrac58,\ \tfrac12,\ \tfrac58,\ \tfrac12,\ \tfrac38,\ 0,$$

in agreement with Section 1; for instance $P(3)=0+1+1=2$ and $T(\tfrac38)=(9-4)/8=\tfrac58$. The totals $P(k)$ are counted digit place by digit place: the digit worth $2^b$ repeats $2^b$ zeros and $2^b$ ones, so $0,\dots,k-1$ contain $\lfloor k/2^{b+1}\rfloor$ complete periods with $2^b$ ones each, and the incomplete last period, of length $q=k\bmod2^{b+1}$, adds $\max(0,q-2^b)$. For $k=8$ this gives $4\cdot1+2\cdot2+1\cdot4=12$; only places with $2^b\lt k$ count.

3. A small copy of the curve in every cell

Split the sum at the index $m$:

$$T(x)=\sum_{n=0}^{m-1}\frac{s(2^nx)}{2^n}+\frac{T(2^mx)}{2^m}.$$

The tail consists of the terms $s(2^j\cdot2^mx)/2^{m+j}$ with $j\ge0$, where $T$ is read with period $1$ like $s$; this is the functional equation $T(x)=s(x)+T(2x)/2$ applied $m$ times. ^[2] On a cell of width $h=2^{-m}$ the first part is, by Section 2, the straight segment that joins the end values of $T$, and the second part is a copy of the whole curve, shrunk by the factor $h$ in both directions. Hence:

(a) on the cell, $T$ lies between the smaller end value and the larger end value plus $h$, since $0\le T\le1$;
(b) the area under $T$ over the cell exceeds the trapezoid $\tfrac h2(T_0+T_1)$ under the segment, with end values $T_0$ and $T_1$, by exactly $h\cdot h\cdot\tfrac12=\tfrac{h^2}2$.

4. The circle and the region inside it

Since $\tfrac1{16}-(x-\tfrac14)^2=x(\tfrac12-x)$, the disc bounded by $C$ lies over $0\le x\le\tfrac12$, and above such an $x$ it is the vertical segment between the lower and the upper arc,

$$y_\pm(x)=\frac12\pm\sqrt{x\left(\tfrac12-x\right)}.$$

The region under the curve meets the same line in $0\le y\le T(x)$, so the two segments overlap in a segment of length

$$\begin{aligned}f(x)&=\max\bigl(0,\ \min(T(x),y_+(x))\\&\qquad{}-\max(0,y_-(x))\bigr),\end{aligned}$$

and, as an area is the integral of the lengths of its vertical slices, the wanted area is $A=\int_0^{1/2}f(x)\,dx$. Both clips are mere safeguards. The lower arc never goes below $y=\tfrac14$, and the curve never rises above the upper arc. For $x\le0.499$, a check of all cells of width $2^{-22}$ with bound (a), the arc being monotone on each cell, leaves a clearance of at least $0.013$. For $t=\tfrac12-x\le0.001$, the functional equation and the symmetry give $T(x)-\tfrac12=T(2t)/2-t$, and with $2^{-m}\le2t\lt2^{1-m}$ the first $m$ terms of $T(2t)$ are at most $2t$ each, as $s(v)\le v$, and the tail at most $2^{-m}\le2t$; so $T(x)-\tfrac12\le mt\lt t\log_2(1/t)$, which lies below the arc's height $\sqrt{t(\tfrac12-t)}$ above $\tfrac12$, as $\sqrt t\log_2(1/t)\le0.32\lt\sqrt{\tfrac12-t}$. So $f=\max(0,T-y_-)$.

5. Where the curve enters the circle

At $x=0$ the curve is below the circle, at $x=\tfrac14$ it passes through the centre, and at $x=\tfrac12$ it passes through the rightmost point $(\tfrac12,\tfrac12)$ of the circle, where $y_-=y_+$. Bisection on $T-y_-$ finds a point where the curve crosses the lower arc:

$$\begin{aligned}x_0&\approx0.078908,\\T(x_0)=y_-(x_0)&\approx0.31772.\end{aligned}$$

For $\tfrac14\lt x\lt\tfrac12$ the first term gives $T(x)\ge x$, and $x-y_-(x)$${}=\sqrt{\tfrac12-x}\,\bigl(\sqrt x-\sqrt{\tfrac12-x}\bigr)\gt0$, so the curve is above the arc there. On $[0,\tfrac14]$, where the lower arc falls, the check with bound (a) on all cells of width $2^{-22}$, refined around $x_0$ down to width $2^{-100}$, puts every cell left of $x_0$ below the arc and every cell right of it above; only an interval of width below $10^{-29}$ around $x_0$ stays undecided. As $|T-y_-|\le1$, it changes the area by less than $10^{-29}$, so up to this amount

$$A=\int_{x_0}^{1/2}\bigl(T(x)-y_-(x)\bigr)\,dx,$$

a scoop between the lower arc and the curve. The programs never need $x_0$; it only explains their accuracy.

6. Simpson's rule on the dyadic grid

Simpson's rule cuts $[a,b]$ into an even number $N$ of steps of width $h$ and integrates, on each pair of steps, the parabola through the three points. A pair with the values $f_0,f_1,f_2$ gives $\tfrac h3(f_0+4f_1+f_2)$, and adding the pairs, whose shared end points get the weight $2$, gives ^[5]

$$\begin{aligned}\frac h3\bigl(f_0&{}+4f_1+2f_2+4f_3\\&{}+\dots+4f_{N-1}+f_N\bigr),\end{aligned}$$

where $f_i$ is the value at the $i$-th point. On $[0,\tfrac12]$ with $h=2^{-m}$ there are $N=2^{m-1}$ steps, an even number for $m\ge2$, and the points $x_i=i/2^m$ are points of Section 2, so $T(x_i)$ is exact, and only the arithmetic with the circle rounds. We call the result $S_m$. Both end values are $0$, because the disc is a single point at $x=0$ and at $x=\tfrac12$. For $m=3$, with $y_-(\tfrac18)=y_-(\tfrac38)=\tfrac12-\tfrac{\sqrt3}8$ and $y_-(\tfrac14)=\tfrac14$, the inner values are $f(\tfrac18)=\tfrac{\sqrt3-1}8\approx0.0915$, $f(\tfrac14)=\tfrac14$ and $f(\tfrac38)=\tfrac{\sqrt3+1}8\approx0.3415$, so

$$\begin{aligned}S_3&=\frac1{24}\Bigl(4\cdot\frac{\sqrt3-1}8+2\cdot\frac14\\&\qquad{}+4\cdot\frac{\sqrt3+1}8\Bigr)\\&=\frac{\sqrt3+\frac12}{24}\approx0.0930.\end{aligned}$$

7. How fast Simpson's rule converges here

For a function with four continuous derivatives the error of composite Simpson's rule is $-(b-a)h^4f^{(4)}(\xi)/180$, so halving $h$ divides it by about $16$. ^[5] The blancmange curve has no derivatives. On a panel $[a,a+2h]$ the rule is a combination of two trapezoid rules:

$$\begin{aligned}&\frac h3(f_0+4f_1+f_2)=\\&\frac{4\cdot\frac h2(f_0+2f_1+f_2)-h(f_0+f_2)}{3}.\end{aligned}$$

Take $f=T$, $h=2^{-m}$ and a panel that is a cell of width $2h$. By Section 3(b) the trapezoid rule with two cells of width $h$ misses $2\cdot h^2/2=h^2$, and the one with a single cell of width $2h$ misses $(2h)^2/2=2h^2$. So Simpson's rule misses exactly $(4h^2-2h^2)/3=2h^2/3$ on every panel, that is, $h/3$ per unit length. By hand on $[0,1]$ with $h=\tfrac14$: the values $0,\tfrac12,\tfrac12,\tfrac12,0$ give $\tfrac1{12}(4\cdot\tfrac12+2\cdot\tfrac12+4\cdot\tfrac12)=\tfrac5{12}$, exactly $\tfrac1{12}$ below $\tfrac12$.

For the area $A$: on the panels left of $x_0$ we have $f=0$, and the rule is exact. On the panels right of the one containing $x_0$ we have $f=T-y_-$; their total length differs from $\tfrac12-x_0$ by at most $2h$, so $T$ causes the shortfall $(\tfrac12-x_0)h/3$ up to $O(h^2)$. The smooth arc adds far less, except near $x=\tfrac12$, where, with $t=\tfrac12-x$, it contains the square root $\sqrt{t(\tfrac12-t)}$: on the last panel the integral and the rule are both of order $h^{3/2}$, and on a panel at distance $t$ the usual error bound is of order $h^5t^{-7/2}$, which summed over $t=2h,4h,\dots$ is again of order $h^{3/2}$. On the panel containing $x_0$, $|T-y_-|$ is of order at most $mh$, as the first $m$ terms of $T$ change by at most $m$ times the distance and the tail by at most $h$, so this panel adds only $O(mh^2)$. Hence

$$\begin{aligned}A-S_m&=\frac{\frac12-x_0}{3}\,h+O\bigl(h^{3/2}\bigr),\\h&=2^{-m}.\end{aligned}$$

The measured $(A-S_m)\cdot2^m$, with the reference value of the verification section, is $0.14039$ at $m=22$ and $0.14037$ at $m=27$, approaching $(\tfrac12-x_0)/3\approx0.14036$. At every level computed, from $m=3$ to $28$, $S_m$ lies below $A$, and on fine grids each new level only halves the error.

8. When to stop

With $c=(\tfrac12-x_0)/3$ and $A-S_m\approx c\,2^{-m}$, the improvement from one level to the next is

$$\begin{aligned}S_m-S_{m-1}&\approx2c\,2^{-m}-c\,2^{-m}\\&=c\,2^{-m}\approx A-S_m,\end{aligned}$$

so the last improvement estimates the remaining error. The programs compute $S_{22},S_{23},\dots$ and stop at the first level whose value differs from the previous one by less than $2\cdot10^{-9}$. The differences for $m=23,\dots,27$ are $1.67\cdot10^{-8}$, $8.37\cdot10^{-9}$, $4.18\cdot10^{-9}$, $2.09\cdot10^{-9}$ and $1.05\cdot10^{-9}$, so they stop at $m=27$, about $1.05\cdot10^{-9}$ below $A$. This could change the eighth decimal only near a rounding boundary, and the reference value shows that $A$ is more than three times as far from the nearest one. Richardson extrapolation ^[6] would remove the leading error: $2S_m-S_{m-1}$ is within $6\cdot10^{-12}$ of $A$ already at $m=22$; the programs do not use it.

How the Code Works

Python and Java follow the C++ program ^[7] step by step, Java with camel-case names and the loop of integrate_chunk inside simpsonLevel. Only C++ has takagi_dyadic, the checks, --skip-checkpoints and the options --eps=, --min-level= and --max-level=, by default $2\cdot10^{-9}$, $22$ and $28$, with levels from $3$ to $30$.

prefix_popcount(n) computes $P(n)$ by the count over digit places of Section 2, for $b=0,1,\dots$ as long as $2^b\lt n$. takagi_dyadic(k, level) evaluates the formula of Section 2 for the checks. overlap_height(x, curve_y) computes $f(x)$ of Section 4: it returns $0$ unless inside, which is $\tfrac1{16}-(x-\tfrac14)^2$, is positive; then its square root dy is the half-length of the circle's vertical segment around $y=\tfrac12$, and the result is the length of the overlap of that segment with the segment from $0$ to curve_y.

integrate_chunk(level, begin, end) walks through the points $i$ from begin to end${}-1$ of level $m$ and keeps the numerator $mi-2P(i)$ of $T(i/2^m)$, a whole number that Python stores as an integer and C++ and Java as a floating-point value. prefix_popcount(begin) starts it, and after each point it grows by $m-2w(i)$, with $w(i)$ from __builtin_popcountll, int.bit_count or Long.bitCount. ^[3] Each point costs one call of overlap_height, added to odd_sum or even_sum by the parity of $i$; no series is ever summed.

simpson_level cuts the inner points $1,\dots,2^{m-1}-1$ into consecutive blocks, sums each with integrate_chunk, adds the block sums and returns h*(f0+fn+4*odd_sum+2*even_sum)/3 with $h=2^{-m}$, f0 $=0$ and fn $=$ overlap_height(0.5, 0.5) $=0$. solve computes $S_{22},S_{23},\dots$ and returns the first $S_m$ with $|S_m-S_{m-1}|\lt$eps, or the last value if none qualifies up to the maximum level. C++ prints the value, converted to double, with eight fixed decimals, Python with the format .8f and Java with String.format and %.8f.

Unless --skip-checkpoints is given, the C++ program first runs run_checkpoints: $P(8)=12$ and $P(16)=32$; $T(\tfrac12)=T(\tfrac14)=\tfrac12$, $T(\tfrac18)=\tfrac38$ and $T(\tfrac38)=\tfrac58$ from takagi_dyadic; and $S_{20},S_{21}\gt0$ with $S_{21}\lt0.2$ and $|S_{21}-S_{20}|\le10^{-4}$. A failure ends the program with exit code $2$.

As $0\le T\le1$, the numerator is an integer between $0$ and $2^m$, so every $x_i$ and $T(x_i)$ is an exact binary fraction; rounding enters only in the circle and the sums. Python and Java use 64-bit floating point, float and double, and C++ uses long double: the 80-bit x87 format with many x86-64 compilers, a 128-bit format on some ARM64 systems, and just double with MSVC ^[8] and, as measured here, on Apple ARM64.

Complexity and Verification

Level $m$ evaluates $2^{m-1}-1$ inner points with constant work each, a population count, a square root and a few operations, so its time is proportional to $2^m$ and the last level dominates: levels $22$ to $27$ evaluate $132\,120\,570$ points, and the C++ checks add $1\,572\,862$. Apart from the block sums the memory is constant. By Section 7 the levels up to an error $\varepsilon$ evaluate about $0.14/\varepsilon$ points, where a smooth function would need about $\varepsilon^{-1/4}$. On a desktop computer the C++ and Java programs take well under a second, the Python program a few seconds.

All three programs print the same eight decimals. Checks run for this article, not part of the published solutions:

with exact fractions: the formula of Section 2 against the definition for every $m\le12$ and $0\le k\le2^m$, prefix_popcount against direct counting for $n\le20000$, and the shortfalls $h^2/2$ per trapezoid cell and $2h^2/3$ per Simpson panel for every $m\le8$;
the cell checks of Sections 4 and 5 on all $2^{21}$ cells of width $2^{-22}$ and their refinement around $x_0$;
an independent vectorised computation of $S_3,\dots,S_{28}$, which agrees with the C++ functions to within $2\cdot10^{-14}$;
a reference value $A_0=I(\tfrac12)-I(x_0)-\int_{x_0}^{1/2}y_-$, computed with $60$-digit arithmetic, where $I(a)=\int_0^aT$ comes from the rules $I(a)=a^2/2+I(2a)/4$ for $a\le\tfrac12$ and $I(a)=\tfrac12-I(1-a)$ for $a\ge\tfrac12$, ^[2] the integral of $y_-$ is taken in closed form and $x_0$ by bisection; by Section 5 it differs from $A$ by less than $10^{-29}$. It agrees with $2S_{27}-S_{26}$ to within $4\cdot10^{-14}$, confirms the error law and the stopping level, and rounds to the same eight decimals as the output of the programs.

Footnotes and References

The references below supply the background facts used above. The exact values at the points $k/2^m$, the shortfall of Simpson's rule for the blancmange curve, the error law and the stopping analysis are derived in this article.

Project Euler 226 — A Scoop of Blancmange. The official statement defines the curve by the sum of s(2ⁿx)/2ⁿ for 0 ≤ x ≤ 1, with s the distance to the nearest integer, gives the area 1/2 under it, and asks for the area under it enclosed by the circle with centre (1/4, 1/2) and radius 1/4, to eight decimals.
Wikipedia — Blancmange curve. Named after Teiji Takagi, whose paper is titled "A Simple Example of the Continuous Function without Derivative". The article gives T(x) = s(x) + T(2x)/2, T(1 − x) = T(x) and, from the integral 1/2 over [0, 1], the rules I(x) = I(2x)/4 + x²/2 for x ≤ 1/2 and I(x) = 1/2 − I(1 − x) for x ≥ 1/2 for I(x) = ∫₀ˣ T.
Wikipedia — Hamming weight. The number of ones in the binary representation of an integer, also called its population count; GCC has __builtin_popcount, Java Long.bitCount, and Python int.bit_count() since version 3.10.
OEIS A000788 — Total number of 1's in binary expansions of 0, ..., n. In the notation of this article a(n) = P(n + 1). The entry remarks that the graph of the sequence is a version of the Takagi curve.
Wikipedia — Simpson's rule. The composite rule for an even number of subintervals has the weights 1, 4, 2, 4, …, 4, 1 times h/3 and the error −(b − a)h⁴f⁽⁴⁾(ξ)/180; for functions that lack derivatives it may give very poor results.
Wikipedia — Richardson extrapolation. If A(h) = A* + Chⁿ + O(hⁿ⁺¹), then (tⁿA(h/t) − A(h))/(tⁿ − 1) approximates A* to order hⁿ⁺¹; for t = 2 and n = 1 this is 2A(h/2) − A(h).
C++ — Euler226.cpp. The linked immutable C++ revision contains prefix_popcount, overlap_height, integrate_chunk, simpson_level, solve and run_checkpoints. It is the implementation record for this article.
cppreference — Fundamental types. The type long double need not follow IEEE-754: some ARM64 implementations use binary128, many x86 and x86-64 ones the x87 80-bit extended format, and MSVC the format of double.

Problem 226 source code

C++

#include <algorithm>
#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <string>
#include <thread>
#include <vector>

namespace {

using u64 = std::uint64_t;

struct Options {
    long double eps = 2e-9L;
    int min_level = 22;
    int max_level = 28;
    int threads = 0;
    bool run_checkpoints = true;
};

bool parse_ld_after_prefix(const std::string& arg, const std::string& prefix, long double& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }

    try {
        value = std::stold(tail);
    } catch (...) {
        return false;
    }
    return value > 0;
}

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }

    int parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<int>(c - '0');
    }
    value = parsed;
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_ld_after_prefix(arg, "--eps=", options.eps) ||
            parse_int_after_prefix(arg, "--min-level=", options.min_level) ||
            parse_int_after_prefix(arg, "--max-level=", options.max_level) ||
            parse_int_after_prefix(arg, "--threads=", options.threads)) {
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }

    return options.eps > 0.0L && options.min_level >= 3 && options.max_level >= options.min_level &&
           options.max_level <= 30 && options.threads >= 0;
}

u64 prefix_popcount(u64 n) {
    u64 total = 0;
    for (int bit = 0; bit < 63; ++bit) {
        const u64 half = 1ULL << bit;
        if (half >= n) {
            break;
        }
        const u64 cycle = half << 1U;
        const u64 full = n / cycle;
        const u64 rem = n % cycle;
        total += full * half;
        if (rem > half) {
            total += rem - half;
        }
    }
    return total;
}

long double takagi_dyadic(const u64 k, const int level) {
    const long double numer = static_cast<long double>(level) * static_cast<long double>(k) -
                              2.0L * static_cast<long double>(prefix_popcount(k));
    const long double denom = static_cast<long double>(1ULL << level);
    return numer / denom;
}

long double overlap_height(const long double x, const long double curve_y) {
    const long double dx = x - 0.25L;
    const long double inside = 0.0625L - dx * dx;
    if (inside <= 0.0L) {
        return 0.0L;
    }

    const long double dy = std::sqrt(inside);
    const long double circle_low = 0.5L - dy;
    const long double circle_high = 0.5L + dy;
    const long double lo = std::max(0.0L, circle_low);
    const long double hi = std::min(circle_high, curve_y);
    return hi > lo ? hi - lo : 0.0L;
}

struct PartialSums {
    long double odd_sum = 0.0L;
    long double even_sum = 0.0L;
};

PartialSums integrate_chunk(const int level, const u64 begin, const u64 end) {
    PartialSums out;
    if (begin >= end) {
        return out;
    }

    const long double inv_denom = 1.0L / static_cast<long double>(1ULL << level);
    long double numer = static_cast<long double>(level) * static_cast<long double>(begin) -
                        2.0L * static_cast<long double>(prefix_popcount(begin));

    for (u64 i = begin; i < end; ++i) {
        const long double x = static_cast<long double>(i) * inv_denom;
        const long double curve_y = numer * inv_denom;
        const long double fx = overlap_height(x, curve_y);
        if ((i & 1ULL) != 0ULL) {
            out.odd_sum += fx;
        } else {
            out.even_sum += fx;
        }
        numer += static_cast<long double>(level - 2 * static_cast<int>(__builtin_popcountll(i)));
    }

    return out;
}

long double simpson_level(const int level, int threads) {
    const u64 denom = 1ULL << level;
    const u64 steps = denom >> 1U;
    if (steps < 2) {
        return 0.0L;
    }

    if (threads <= 0) {
        threads = static_cast<int>(std::thread::hardware_concurrency());
        if (threads <= 0) {
            threads = 1;
        }
    }

    const u64 interior = steps - 1;
    if (interior == 0) {
        return 0.0L;
    }
    threads = std::max(1, std::min<int>(threads, static_cast<int>(interior)));

    std::vector<PartialSums> partials(static_cast<std::size_t>(threads));
    std::vector<std::thread> workers;
    workers.reserve(static_cast<std::size_t>(threads));

    const u64 chunk = (interior + static_cast<u64>(threads) - 1ULL) / static_cast<u64>(threads);
    for (int t = 0; t < threads; ++t) {
        const u64 begin = 1ULL + static_cast<u64>(t) * chunk;
        const u64 end = std::min(steps, begin + chunk);
        if (begin >= end) {
            continue;
        }
        workers.emplace_back([&, t, begin, end]() { partials[static_cast<std::size_t>(t)] = integrate_chunk(level, begin, end); });
    }

    for (std::thread& worker : workers) {
        worker.join();
    }

    long double odd_sum = 0.0L;
    long double even_sum = 0.0L;
    for (const PartialSums& p : partials) {
        odd_sum += p.odd_sum;
        even_sum += p.even_sum;
    }

    const long double h = 1.0L / static_cast<long double>(denom);
    const long double f0 = 0.0L;
    const long double fn = overlap_height(0.5L, 0.5L);
    return h * (f0 + fn + 4.0L * odd_sum + 2.0L * even_sum) / 3.0L;
}

long double solve(const long double eps, const int min_level, const int max_level, const int threads) {
    long double previous = simpson_level(min_level, threads);
    for (int level = min_level + 1; level <= max_level; ++level) {
        const long double current = simpson_level(level, threads);
        if (std::fabsl(current - previous) < eps) {
            return current;
        }
        previous = current;
    }
    return previous;
}

bool run_checkpoints() {
    if (prefix_popcount(8) != 12 || prefix_popcount(16) != 32) {
        std::cerr << "Checkpoint failed for prefix popcount" << '\n';
        return false;
    }
    if (std::fabsl(takagi_dyadic(1, 1) - 0.5L) > 1e-15L || std::fabsl(takagi_dyadic(1, 2) - 0.5L) > 1e-15L ||
        std::fabsl(takagi_dyadic(1, 3) - 0.375L) > 1e-15L || std::fabsl(takagi_dyadic(3, 3) - 0.625L) > 1e-15L) {
        std::cerr << "Checkpoint failed for dyadic Takagi values" << '\n';
        return false;
    }

    const long double coarse = simpson_level(20, 1);
    const long double fine = simpson_level(21, 1);
    if (coarse <= 0.0L || fine <= 0.0L || fine >= 0.2L || std::fabsl(fine - coarse) > 1e-4L) {
        std::cerr << "Checkpoint failed for area refinement" << '\n';
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }

    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }

    const long double answer = solve(options.eps, options.min_level, options.max_level, options.threads);
    std::cout << std::fixed << std::setprecision(8) << static_cast<double>(answer) << '\n';
    return 0;
}

Python

import math
import multiprocessing

def prefix_popcount(n):
    total = 0
    for bit in range(63):
        half = 1 << bit
        if half >= n:
            break
        cycle = half << 1
        full = n // cycle
        rem = n % cycle
        total += full * half
        if rem > half:
            total += rem - half
    return total

def overlap_height(x, curve_y):
    dx = x - 0.25
    inside = 0.0625 - dx * dx
    if inside <= 0.0:
        return 0.0
    
    dy = math.sqrt(inside)
    circle_low = 0.5 - dy
    circle_high = 0.5 + dy
    lo = max(0.0, circle_low)
    hi = min(circle_high, curve_y)
    return hi - lo if hi > lo else 0.0

def integrate_chunk(args):
    level, begin, end = args
    if begin >= end:
        return 0.0, 0.0
        
    inv_denom = 1.0 / (1 << level)
    numer = level * begin - 2 * prefix_popcount(begin)
    
    odd_sum = 0.0
    even_sum = 0.0
    
    for i in range(begin, end):
        x = i * inv_denom
        curve_y = numer * inv_denom
        fx = overlap_height(x, curve_y)
        
        if i & 1:
            odd_sum += fx
        else:
            even_sum += fx
            
        numer += level - 2 * i.bit_count()
        
    return odd_sum, even_sum

def simpson_level(level, threads):
    denom = 1 << level
    steps = denom >> 1
    if steps < 2: return 0.0
    
    interior = steps - 1
    if interior == 0: return 0.0
    
    threads = min(threads, interior)
    chunk = (interior + threads - 1) // threads
    
    tasks = []
    for t in range(threads):
        begin = 1 + t * chunk
        end = min(steps, begin + chunk)
        if begin < end:
            tasks.append((level, begin, end))
            
    odd_sum = 0.0
    even_sum = 0.0
    
    if threads > 1 and len(tasks) > 1:
        with multiprocessing.Pool(threads) as pool:
            results = pool.map(integrate_chunk, tasks)
            for os, es in results:
                odd_sum += os
                even_sum += es
    else:
        for t in tasks:
            os, es = integrate_chunk(t)
            odd_sum += os
            even_sum += es
            
    h = 1.0 / denom
    f0 = 0.0
    fn = overlap_height(0.5, 0.5)
    return h * (f0 + fn + 4.0 * odd_sum + 2.0 * even_sum) / 3.0

def solve(eps=2e-9, min_level=22, max_level=28):
    threads = multiprocessing.cpu_count() or 1
    
    previous = simpson_level(min_level, threads)
    for level in range(min_level + 1, max_level + 1):
        current = simpson_level(level, threads)
        if abs(current - previous) < eps:
            return f"{current:.8f}"
        previous = current
    return f"{previous:.8f}"

if __name__ == '__main__':
    print(solve())

Java

import java.util.ArrayList;
import java.util.List;
import java.util.concurrent.*;

public class Euler226 {
    static long prefixPopcount(long n) {
        long total = 0;
        for (int bit = 0; bit < 63; ++bit) {
            long half = 1L << bit;
            if (half >= n)
                break;
            long cycle = half << 1;
            long full = n / cycle;
            long rem = n % cycle;
            total += full * half;
            if (rem > half) {
                total += rem - half;
            }
        }
        return total;
    }

    static double overlapHeight(double x, double curveY) {
        double dx = x - 0.25;
        double inside = 0.0625 - dx * dx;
        if (inside <= 0.0)
            return 0.0;

        double dy = Math.sqrt(inside);
        double circleLow = 0.5 - dy;
        double circleHigh = 0.5 + dy;
        double lo = Math.max(0.0, circleLow);
        double hi = Math.min(circleHigh, curveY);
        return hi > lo ? hi - lo : 0.0;
    }

    static class PartialSums {
        double oddSum = 0.0;
        double evenSum = 0.0;
    }

    public static String solve() {
        double eps = 2e-9;
        int minLevel = 22;
        int maxLevel = 28;
        int threads = Math.max(1, Runtime.getRuntime().availableProcessors());

        double previous = simpsonLevel(minLevel, threads);
        for (int level = minLevel + 1; level <= maxLevel; ++level) {
            double current = simpsonLevel(level, threads);
            if (Math.abs(current - previous) < eps) {
                return String.format(java.util.Locale.US, "%.8f", current);
            }
            previous = current;
        }
        return String.format(java.util.Locale.US, "%.8f", previous);
    }

    static double simpsonLevel(int level, int numThreads) {
        long denom = 1L << level;
        long steps = denom >> 1;
        if (steps < 2)
            return 0.0;

        long interior = steps - 1;
        if (interior == 0)
            return 0.0;

        int threads = (int) Math.min(numThreads, interior);
        long chunk = (interior + threads - 1) / threads;

        ExecutorService executor = Executors.newFixedThreadPool(threads);
        List<Future<PartialSums>> futures = new ArrayList<>();

        for (int t = 0; t < threads; ++t) {
            final long begin = 1 + t * chunk;
            final long end = Math.min(steps, begin + chunk);
            if (begin >= end)
                continue;

            futures.add(executor.submit(() -> {
                PartialSums p = new PartialSums();
                double invDenom = 1.0 / denom;
                double numer = (double) level * begin - 2.0 * prefixPopcount(begin);

                for (long i = begin; i < end; ++i) {
                    double x = i * invDenom;
                    double curveY = numer * invDenom;
                    double fx = overlapHeight(x, curveY);

                    if ((i & 1) != 0) {
                        p.oddSum += fx;
                    } else {
                        p.evenSum += fx;
                    }

                    numer += level - 2 * Long.bitCount(i);
                }
                return p;
            }));
        }

        double oddSum = 0.0;
        double evenSum = 0.0;

        for (Future<PartialSums> f : futures) {
            try {
                PartialSums p = f.get();
                oddSum += p.oddSum;
                evenSum += p.evenSum;
            } catch (Exception e) {
            }
        }
        executor.shutdown();

        double h = 1.0 / denom;
        double f0 = 0.0;
        double fn = overlapHeight(0.5, 0.5);
        return h * (f0 + fn + 4.0 * oddSum + 2.0 * evenSum) / 3.0;
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}