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

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

Project Euler problem page: Project Euler 226
Takagi function: Wikipedia - Takagi function
Hamming weight and binary digit sums: Wikipedia - Hamming weight
Simpson's rule: Wikipedia - Simpson's rule
Numerical integration: Wikipedia - Numerical integration

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());
    }
}