Project Euler

A Scoop of Blancmange

The Blancmange curve is blanc(x)=sum_(n=0)^infinity(s(2^n x))/(2^n), s(x)=min_(kinmathbb Z)|x-k|. Consider the circle (x-frac14)^2+(y-frac12)^2=frac116. Find the area of the region that lies both b...

Source sync May 21, 2026

Problem #0226

Level Level 10

Solved By 2,047

Languages C++, Python

Answer 0.11316017

Length 254 words

geometrymodular_arithmeticalgebra

The blancmange curve is the set of points $(x, y)$ such that $0 \le x \le 1$ and $\displaystyle y = \sum \limits_{n = 0}^{\infty} {\dfrac{s(2^n x)}{2^n}}$, where $s(x)$ is the distance from $x$ to the nearest integer.

The area under the blancmange curve is equal to Â½, shown in pink in the diagram below.

Problem illustration

Let $C$ be the circle with centre $\left ( \frac{1}{4}, \frac{1}{2} \right )$ and radius $\frac{1}{4}$, shown in black in the diagram.

What area under the blancmange curve is enclosed by $C$?

Give your answer rounded to eight decimal places in the form $0.abcdefgh$.

Problem 226: A Scoop of Blancmange

Mathematical Development

Since $0 \le s(t) \le \tfrac12$ , the tail of the Blancmange series is bounded by

\sum_{n=N}^{\infty}\frac{1}{2^{n+1}}=2^{-N}.

So truncating after $60$ terms leaves an error below $10^{-18}$ , which is far smaller than the required precision.

The lower arc of the circle is

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

The relevant intersection points satisfy

\left(x-\frac14\right)^2+\left(\operatorname{blanc}(x)-\frac12\right)^2=\frac1{16}.

One of them is obvious:

x=\frac12, \qquad \operatorname{blanc}\!\left(\frac12\right)=\frac12.

The other lies near $x \approx 0.0789$ and is found numerically. Once that left intersection $x_1$ is known, the desired area is

\int_{x_1}^{1/2}\bigl(\operatorname{blanc}(x)-y_{\mathrm{low}}(x)\bigr)\,dx.

So the whole problem reduces to two numerical tasks:

find the left intersection accurately;
integrate the difference between the curve and the lower arc.

Editorial

The function itself is not the obstacle. Because the Blancmange series has a geometric tail, $60$ terms already give much more accuracy than the final answer needs. That lets us treat $\operatorname{blanc}(x)$ as an ordinary smooth-enough numeric function everywhere in the relevant interval.

From there, the geometry is clean. The enclosed region starts at the nontrivial intersection of the curve with the circle, ends at $x=\tfrac12$ , and its height is simply

\operatorname{blanc}(x)-y_{\mathrm{low}}(x).

So the program first locates the left intersection with a robust one-dimensional root finder, then integrates that height with Simpson’s rule on a fine grid. There is no symbolic trick after that; the important part is keeping the truncation and quadrature errors comfortably below $10^{-8}$ .

Pseudocode

Approximate blanc(x) by summing the first 60 terms of the series.

Define the circle equation
    f(x) = (x - 1/4)^2 + (blanc(x) - 1/2)^2 - 1/16.

Use a bracketing root finder on the interval near 0.08
to locate the left intersection x1.

Define the lower circle arc y_low(x).
Define the integrand blanc(x) - y_low(x).

Apply Simpson's rule on [x1, 1/2] with a sufficiently fine subdivision.
Print the area to 8 decimal places.

Complexity Analysis

Time: $O(T \cdot M)$ where $T=60$ Blancmange terms are evaluated at each of the $M$ quadrature points.
Space: $O(1)$ .

Answer

\boxed{0.11316017}

C++ project_euler/problem_226/solution.cpp

#include <bits/stdc++.h>
using namespace std;

// Triangle wave: distance to nearest integer
double s(double x) {
    return abs(x - round(x));
}

// Blancmange curve value at x (60 terms suffice for ~18 digit precision)
double blanc(double x) {
    double result = 0.0;
    double pow2 = 1.0;
    for (int n = 0; n < 60; n++) {
        result += s(pow2 * x) / pow2;
        pow2 *= 2.0;
    }
    return result;
}

// Circle: (x - 1/4)^2 + (y - 1/2)^2 = 1/16
// Lower arc: y = 1/2 - sqrt(1/16 - (x - 1/4)^2)
double circle_lower(double x) {
    double val = 1.0 / 16.0 - (x - 0.25) * (x - 0.25);
    if (val < 0) return 0.5;
    return 0.5 - sqrt(val);
}

// Function whose root gives the intersection point
double on_circle(double x) {
    double b = blanc(x);
    return (x - 0.25) * (x - 0.25) + (b - 0.5) * (b - 0.5) - 1.0 / 16.0;
}

// Brent's method to find root in [a, b]
double brent(double a, double b, double tol = 1e-15) {
    double fa = on_circle(a), fb = on_circle(b);
    if (fa * fb > 0) return a;
    double c = a, fc = fa, d = b - a, e = d;
    for (int i = 0; i < 200; i++) {
        if (fb * fc > 0) { c = a; fc = fa; d = e = b - a; }
        if (fabs(fc) < fabs(fb)) { a = b; b = c; c = a; fa = fb; fb = fc; fc = fa; }
        double tol1 = 2e-16 * fabs(b) + 0.5 * tol;
        double m = 0.5 * (c - b);
        if (fabs(m) <= tol1 || fb == 0) return b;
        if (fabs(e) >= tol1 && fabs(fa) > fabs(fb)) {
            double s_val;
            if (a == c) {
                s_val = fb / fa; double p = 2 * m * s_val; double q = 1 - s_val;
                if (p > 0) q = -q; else p = -p;
                if (2*p < 3*m*q - fabs(tol1*q) && 2*p < fabs(e*q)) { e = d; d = p/q; }
                else { d = m; e = m; }
            } else { d = m; e = m; }
        } else { d = m; e = m; }
        a = b; fa = fb;
        if (fabs(d) > tol1) b += d;
        else b += (m > 0 ? tol1 : -tol1);
        fb = on_circle(b);
    }
    return b;
}

int main() {
    // Find intersection point x1
    double x1 = brent(0.05, 0.1);
    double x2 = 0.5;

    // Integrand: blanc(x) - circle_lower(x)
    auto integrand = [](double x) -> double {
        return blanc(x) - (0.5 - sqrt(max(0.0, 1.0 / 16.0 - (x - 0.25) * (x - 0.25))));
    };

    // Simpson's rule with 2M subdivisions
    int N = 2000000;
    double h = (x2 - x1) / N;
    double sum = integrand(x1) + integrand(x2);

    for (int i = 1; i < N; i++) {
        double x = x1 + i * h;
        sum += (i % 2 == 0 ? 2.0 : 4.0) * integrand(x);
    }

    double area = sum * h / 3.0;
    printf("%.8f\n", area);
    return 0;
}

Python project_euler/problem_226/solution.py

"""
Problem 226: A Scoop of Blancmange

Find the area under the Blancmange curve enclosed by the circle
(x - 1/4)^2 + (y - 1/2)^2 = 1/16, to 8 decimal places.
"""

import math

def s(x):
    """Triangle wave: distance to nearest integer."""
    return abs(x - round(x))

def blanc(x, terms=60):
    """Blancmange (Takagi) curve value at x."""
    result = 0.0
    pow2 = 1.0
    for _ in range(terms):
        result += s(pow2 * x) / pow2
        pow2 *= 2.0
    return result

def on_circle(x):
    """Returns (x-1/4)^2 + (blanc(x)-1/2)^2 - 1/16. Zero on circle."""
    b = blanc(x)
    return (x - 0.25)**2 + (b - 0.5)**2 - 1.0/16

def circle_lower(x):
    """Lower arc of circle: y = 1/2 - sqrt(1/16 - (x-1/4)^2)."""
    val = 1.0/16 - (x - 0.25)**2
    if val < 0:
        return 0.5
    return 0.5 - math.sqrt(val)

def bisect_root(f, a, b, tol=1e-14):
    """Find root of f in [a, b] using bisection."""
    fa, fb = f(a), f(b)
    if fa * fb > 0:
        return a
    for _ in range(100):
        mid = (a + b) / 2
        fm = f(mid)
        if abs(b - a) < tol:
            return mid
        if fa * fm <= 0:
            b, fb = mid, fm
        else:
            a, fa = mid, fm
    return (a + b) / 2

def integrand(x):
    """blanc(x) - circle_lower(x)"""
    return blanc(x) - circle_lower(x)

def simpsons(f, a, b, n):
    """Simpson's rule with n subdivisions (n must be even)."""
    h = (b - a) / n
    result = f(a) + f(b)
    for i in range(1, n):
        x = a + i * h
        result += (2 if i % 2 == 0 else 4) * f(x)
    return result * h / 3

def solve():
    x1 = bisect_root(on_circle, 0.05, 0.1)
    x2 = 0.5
    subdivisions = 100000
    area = simpsons(integrand, x1, x2, subdivisions)
    print(f"{area:.8f}")


if __name__ == "__main__":
    solve()