Detailed mathematical approach

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.

An Exact Formula Using Prefix Popcounts

Let \(s_2(r)\) be the number of 1-bits in the binary expansion of \(r\), and define the prefix sum

$P(k)=\sum_{r=0}^{k-1}s_2(r).$

For dyadic points one has the exact identity

$B\left(\frac{k}{2^m}\right)=\frac{mk-2P(k)}{2^m}.$

This is the central formula used by all three implementations. It immediately yields the recurrence

$B\left(\frac{k+1}{2^m}\right)-B\left(\frac{k}{2^m}\right)=\frac{m-2s_2(k)}{2^m},$

because \(P(k+1)=P(k)+s_2(k)\). So once one dyadic value is known, the next one is obtained by a constant-time update controlled only by the bit count of the current index.

Turning the Geometry into a Vertical Overlap Integral

The circle is

$\left(x-\frac14\right)^2+\left(y-\frac12\right)^2=\left(\frac14\right)^2.$

It occupies only the horizontal range \(0 \le x \le \tfrac12\). For such an \(x\), the corresponding vertical slice of the circle is

$y_{\pm}(x)=\frac12 \pm \sqrt{\frac1{16}-\left(x-\frac14\right)^2}.$

The blancmange region is the set

$R_B=\{(x,y): 0 \le y \le B(x)\}.$

Hence the common vertical height at position \(x\) is

$h(x)=\max\!\left(0,\ \min\!\bigl(B(x),y_+(x)\bigr)-\max\!\bigl(0,y_-(x)\bigr)\right).$

The required area is therefore

$A=\int_0^{1/2} h(x)\,dx.$

Outside \([0,\tfrac12]\) the circle contributes nothing, so the integral really is one-dimensional and compactly supported.

Worked Example: \(x=\tfrac38\)

Take \(m=3\) and \(k=3\), so \(x=3/8\). The prefix bit-count sum is

$P(3)=s_2(0)+s_2(1)+s_2(2)=0+1+1=2,$

and therefore

$B\left(\frac38\right)=\frac{3 \cdot 3-2 \cdot 2}{8}=\frac58.$

This agrees with the direct Takagi sum:

$B\left(\frac38\right)=\phi\left(\frac38\right)+\frac12\phi\left(\frac34\right)+\frac14\phi\left(\frac32\right)=\frac38+\frac18+\frac18=\frac58.$

For the circle slice, \(x-\tfrac14=\tfrac18\), so

$y_-\left(\frac38\right)=\frac12-\frac{\sqrt3}{8},\qquad y_+\left(\frac38\right)=\frac12+\frac{\sqrt3}{8}.$

Since \(\tfrac58 < y_+(\tfrac38)\), the overlap height is

$h\left(\frac38\right)=\frac58-\left(\frac12-\frac{\sqrt3}{8}\right)=\frac{1+\sqrt3}{8}.$

This is exactly what the integrand measures at each sample point: the portion of the circle's vertical segment that lies below the blancmange curve and above the \(x\)-axis.

Why Simpson Refinement Fits the Problem

At refinement level \(m\), the interval \([0,\tfrac12]\) is divided into \(N=2^{m-1}\) equal subintervals of width \(h=2^{-m}\), so every node is dyadic:

$x_i=\frac{i}{2^m},\qquad 0 \le i \le 2^{m-1}.$

That means each sampled value \(B(x_i)\) is computed exactly by the formula above. The only approximation comes from the quadrature itself. Simpson's rule gives

$A_m=\frac{h}{3}\left(h(x_0)+h(x_N)+4\sum_{\substack{1 \le i \le N-1 \\ i\text{ odd}}} h(x_i)+2\sum_{\substack{2 \le i \le N-2 \\ i\text{ even}}} h(x_i)\right).$

The implementations evaluate this on successively finer dyadic grids and stop when two consecutive levels differ by less than the requested tolerance.

How the Code Works

Counting Binary Digits in Bulk

To evaluate \(P(k)\) quickly, the implementations count 1-bits one position at a time. For bit position \(b\), the pattern repeats every \(2^{b+1}\) integers, and the number of ones among \(0,1,\dots,k-1\) is

$\left\lfloor\frac{k}{2^{b+1}}\right\rfloor 2^b+\max\!\left(0,\ k \bmod 2^{b+1}-2^b\right).$

Summing that over bit positions yields \(P(k)\) in \(O(\log k)\) arithmetic time rather than by iterating through all previous integers.

Streaming the Dyadic Blancmange Values

During one Simpson sweep at fixed level \(m\), the implementations maintain the numerator \(mk-2P(k)\) rather than recomputing it from scratch at every node. Moving from index \(k\) to \(k+1\) changes this numerator by \(m-2s_2(k)\), so the next blancmange value is obtained with one bit count and a few arithmetic operations. This invariant is the reason a very fine dyadic grid is still practical.

Accumulating and Refining the Integral

For each sample point, the code computes the circle's lower and upper \(y\)-values, clips that interval against \([0,B(x)]\), and adds the resulting overlap height to either the odd or the even Simpson accumulator. The interior sample range is split into independent chunks so partial sums can be formed separately and then combined. The C++, Python, and Java implementations all use the same mathematics; the compiled versions divide the work across threads, and the Python version can distribute chunks across worker processes when that is worthwhile. After one level is finished, the next level is computed and compared against it; refinement stops as soon as the tolerance test passes.

Complexity Analysis

At level \(m\), Simpson's rule on \([0,\tfrac12]\) uses \(2^{m-1}+1\) sample points, so one refinement level costs \(O(2^m)\) time. Because the dyadic blancmange values are streamed by the recurrence above, each additional node is only a constant amount of arithmetic plus a bit count and a square root.

The extra memory is \(O(T)\) for \(T\) worker partial sums, or \(O(1)\) in a serial evaluation. Since the code checks a short sequence of increasing dyadic levels, the total running time up to convergence is still on the order of the finest level that gets evaluated.

Footnotes and References

1. Project Euler problem page: Project Euler 226
2. Takagi function: Wikipedia - Takagi function
3. Hamming weight and binary digit sums: Wikipedia - Hamming weight
4. Simpson's rule: Wikipedia - Simpson's rule
5. Numerical integration: Wikipedia - Numerical integration