Problem 226: A Scoop of Blancmange

Mathematical Foundation

Theorem 1 (Convergence of the Blancmange Series). The series $\operatorname{blanc}(x) = \sum_{n=0}^{\infty} \frac{s(2^n x)}{2^n}$ converges uniformly on $\mathbb{R}$.

Proof. Since $0 \leq s(t) \leq 1/2$ for all $t$, each term satisfies $0 \leq \frac{s(2^n x)}{2^n} \leq \frac{1}{2^{n+1}}$. The series $\sum_{n=0}^{\infty} \frac{1}{2^{n+1}} = 1$ converges, so by the Weierstrass $M$-test, $\operatorname{blanc}(x)$ converges uniformly. $\square$

Theorem 2 (Continuity and Non-Differentiability). $\operatorname{blanc}(x)$ is continuous everywhere and differentiable nowhere on $[0,1]$.

Proof. Continuity follows from the uniform limit of continuous functions $s(2^n x)/2^n$. Nowhere-differentiability is the classical Takagi theorem (1903). $\square$

Lemma 1 (Truncation Error). Truncating the series at $N$ terms gives error bounded by $2^{-N}$:

Proof. Direct from $|s(t)| \leq 1/2$ and geometric series summation. $\square$

Theorem 3 (Symmetry). $\operatorname{blanc}(x) = \operatorname{blanc}(1 - x)$ for all $x$.

Proof. $s(t) = s(1 - t)$ for all $t$ (distance to nearest integer is symmetric about $1/2$ modulo 1). Hence $s(2^n x) = s(2^n(1-x))$ for all $n$ (since $2^n(1-x) = 2^n - 2^n x$ and $s$ is periodic with period 1). $\square$

Lemma 2 (Circle Arcs). The circle $C$ has lower and upper arcs on $[0, 1/2]$:

Proof. Solve $(x - 1/4)^2 + (y - 1/2)^2 = 1/16$ for $y$. $\square$

Theorem 4 (Intersection and Area). The Blancmange curve intersects the circle at $x_2 = 1/2$ (where $\operatorname{blanc}(1/2) = 1/2$) and at some $x_1 \approx 0.0789$. Between $x_1$ and $x_2$, the curve lies above the lower arc and inside $C$. The enclosed area is

Proof. At $x = 1/2$: $\operatorname{blanc}(1/2) = 1/2$ and $(1/2 - 1/4)^2 + (1/2 - 1/2)^2 = 1/16$, confirming the point lies on $C$. The left intersection $x_1$ is the unique root of $(x - 1/4)^2 + (\operatorname{blanc}(x) - 1/2)^2 = 1/16$ on $(0, 1/2)$ (verified numerically). Between $x_1$ and $x_2$, the Blancmange curve lies above $y_{\text{low}}(x)$ and below $y_{\text{up}}(x)$, so the area between the curve and the lower arc gives the enclosed region. $\square$

Editorial

We find x1 by Brent’s method on f(x) = (x-1/4)^2 + (blanc(x)-1/2)^2 - 1/16. Finally, simpson’s rule with M subdivisions.

Pseudocode

Find x1 by Brent's method on f(x) = (x-1/4)^2 + (blanc(x)-1/2)^2 - 1/16
Simpson's rule with M subdivisions

Complexity Analysis

• Time: $O(N \cdot M)$ where $N = 60$ (series terms) and $M = 2 \times 10^6$ (quadrature points). Total: $\approx 1.2 \times 10^8$ floating-point operations.
• Space: $O(1)$ (each quadrature point is computed on the fly).

Answer

#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;
}

"""
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

# Find intersection point
x1 = bisect_root(on_circle, 0.05, 0.1)
x2 = 0.5

# Compute area using Simpson's rule with fine grid
N = 100000  # 100K subdivisions for reasonable Python speed
area = simpsons(integrand, x1, x2, N)

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

# Optional visualization