Project Euler Problem 226

Solution

Answer: 0.11316017

Let

$$B(x)=\sum_{n=0}^{\infty}\frac{s(2^n x)}{2^n}, \qquad s(x)=\text{distance from }x\text{ to the nearest integer}.$$

This is the blancmange function (or Takagi function).

We seek the area enclosed simultaneously:

• under the blancmange curve $y=B(x)$,
• and inside the circle

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

The required region is the small “scoop” bounded by the lower semicircle and the blancmange curve.

---

1. Mathematical analysis

Step 1 — The circle equation

The lower half of the circle is

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

We need the intersection point $x_0$ satisfying

$$B(x_0)=y_c(x_0).$$

By symmetry of the blancmange curve,

$$B(x)=B(1-x),$$

and the relevant enclosed region runs from $x_0$ to $1/2$.

Hence the desired area is

$$A=\int_{x_0}^{1/2}\left(B(x)-y_c(x)\right),dx.$$

---

Step 2 — Efficient evaluation of $B(x)$

The definition converges very rapidly because each term is divided by $2^n$.

For numerical computation:

$$B(x)\approx \sum_{n=0}^{N-1}\frac{s(2^n x)}{2^n},$$

with $N=50$ already giving extremely high precision.

For any real $t$,

$$s(t)=|t-\operatorname{round}(t)|.$$

---

Step 3 — Finding the intersection

Define

$$f(x)=B(x)-y_c(x).$$

Using bisection on the interval $(0,0.5)$, we obtain

$$x_0 \approx 0.07890778796534.$$

---

Step 4 — Numerical integration

We compute

$$A=\int_{x_0}^{1/2}(B(x)-y_c(x)),dx$$

using Simpson’s rule with a fine subdivision.

The numerical value converges to

$$A \approx 0.11316016345\ldots$$

which rounds to eight decimal places as

$$0.11316017.$$

---

2. Python implementation

from math import sqrt

# Distance to nearest integer
def s(x):
    return abs(x - round(x))

# Blancmange function
def B(x, terms=60):
    total = 0.0
    scale = 1.0

    for n in range(terms):
        total += s(scale * x) / scale
        scale *= 2.0

    return total

# Lower semicircle
def circle_lower(x):
    return 0.5 - sqrt(0.0625 - (x - 0.25) ** 2)

# Difference between curves
def f(x):
    return B(x) - circle_lower(x)

# Bisection to find intersection point
a, b = 0.0, 0.5

for _ in range(100):
    m = (a + b) / 2

    if f(m) > 0:
        b = m
    else:
        a = m

x0 = (a + b) / 2

# Simpson integration
def simpson(g, a, b, n=20000):
    if n % 2 == 1:
        n += 1

    h = (b - a) / n
    total = g(a) + g(b)

    for i in range(1, n):
        coeff = 4 if i % 2 == 1 else 2
        total += coeff * g(a + i * h)

    return total * h / 3

area = simpson(
    lambda x: B(x) - circle_lower(x),
    x0,
    0.5,
    20000
)

print(f"{area:.8f}")

---

3. Code walkthrough

Blancmange function

def B(x, terms=60):

Computes the truncated infinite series.

Because terms decay like $2^{-n}$, 60 terms is far more than enough.

---

Circle

def circle_lower(x):

Implements the lower branch of the circle equation.

---

Root finding

for _ in range(100):

Performs bisection to locate the intersection point $x_0$.

---

Integration

simpson(...)

Uses Simpson’s rule for accurate numerical integration of

$$B(x)-y_c(x).$$

---

4. Final verification

The answer should be:

• positive,
• much smaller than the total blancmange area $1/2$,
• roughly the size of the visible “scoop”.

Our computed value

$$0.113160163\ldots$$

fits these expectations and is numerically stable under increased precision.

---

Answer: 0.11316017