Problem 392: Enmeshed Unit Circle

View on Project Euler

Project Euler Problem 392 Solution

EulerSolve provides an optimized solution for Project Euler Problem 392, Enmeshed Unit Circle, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For even \(N\), the implementation models the red region as a symmetric staircase built around the unit circle. By symmetry, the four quadrants contribute the same amount, so it is enough to optimize a single quadrant and multiply by \(4\). If \(m=N/2\), the numerical task is to choose \(m\) interior breakpoints on the quarter-circle so that this staircase area is as small as possible. Mathematical Approach 1. Reduce the geometry to one quadrant In the first quadrant the circle is the graph $$f(x)=\sqrt{1-x^2}, \qquad 0 \le x \le 1.$$ Choose breakpoints $$0=x_0 \lt x_1 \lt \cdots \lt x_m \lt x_{m+1}=1,$$ and set \(y_i=f(x_i)\). On each interval \([x_i,x_{i+1}]\), the staircase uses the constant height \(y_i\). Therefore the quadrant contribution is $$Q(x_1,\dots,x_m)=\sum_{i=0}^{m}(x_{i+1}-x_i)f(x_i).$$ The full red area is then $$A_N=4Q(x_1,\dots,x_m).$$ This is exactly the quantity accumulated by the code after the optimal breakpoints have been found. 2. Derive the optimality recurrence The key point is that the solver is not using an arbitrary mesh. It chooses the breakpoints so that \(Q\) is stationary with respect to every interior variable \(x_k\). Only two summands depend on \(x_k\): $$ (x_k-x_{k-1})f(x_{k-1}) \qquad \text{and} \qquad (x_{k+1}-x_k)f(x_k)....

Detailed mathematical approach

Problem Summary

For even \(N\), the implementation models the red region as a symmetric staircase built around the unit circle. By symmetry, the four quadrants contribute the same amount, so it is enough to optimize a single quadrant and multiply by \(4\). If \(m=N/2\), the numerical task is to choose \(m\) interior breakpoints on the quarter-circle so that this staircase area is as small as possible.

Mathematical Approach

1. Reduce the geometry to one quadrant

In the first quadrant the circle is the graph

$$f(x)=\sqrt{1-x^2}, \qquad 0 \le x \le 1.$$

Choose breakpoints

$$0=x_0 \lt x_1 \lt \cdots \lt x_m \lt x_{m+1}=1,$$

and set \(y_i=f(x_i)\). On each interval \([x_i,x_{i+1}]\), the staircase uses the constant height \(y_i\). Therefore the quadrant contribution is

$$Q(x_1,\dots,x_m)=\sum_{i=0}^{m}(x_{i+1}-x_i)f(x_i).$$

The full red area is then

$$A_N=4Q(x_1,\dots,x_m).$$

This is exactly the quantity accumulated by the code after the optimal breakpoints have been found.

2. Derive the optimality recurrence

The key point is that the solver is not using an arbitrary mesh. It chooses the breakpoints so that \(Q\) is stationary with respect to every interior variable \(x_k\). Only two summands depend on \(x_k\):

$$ (x_k-x_{k-1})f(x_{k-1}) \qquad \text{and} \qquad (x_{k+1}-x_k)f(x_k). $$

Differentiating \(Q\) with respect to \(x_k\) gives

$$\frac{\partial Q}{\partial x_k}=f(x_{k-1})-f(x_k)+(x_{k+1}-x_k)f'(x_k).$$

For the unit circle,

$$f'(x)=\frac{-x}{\sqrt{1-x^2}}=-\frac{x}{f(x)}.$$

Setting \(\frac{\partial Q}{\partial x_k}=0\) and rearranging yields

$$\boxed{x_{k+1}=x_k+\frac{(f(x_{k-1})-f(x_k))f(x_k)}{x_k}}, \qquad 1 \le k \le m.$$

This boxed identity is the exact recurrence implemented in the C++, Python, and Java solutions.

3. Why one unknown is enough: the shooting formulation

Once \(x_1\) is fixed, the recurrence determines \(x_2,x_3,\dots,x_{m+1}\) one after another. That reduces the original \(m\)-variable minimization problem to a one-variable boundary condition:

$$x_{m+1}=1.$$

The implementations encode this through the residual function

$$R(x_1)=x_{m+1}(x_1)-1.$$

Then a shooting step becomes: guess \(x_1\), build the whole sequence, inspect the sign of \(R(x_1)\), and refine the guess. The code brackets \(x_1\) inside \((10^{-18},0.999999999999)\) and performs 260 bisection iterations. If a trial leaves the admissible region or produces a non-finite value, the build is rejected and the residual is treated as \(+\infty\), which safely pushes the search back toward valid values.

4. Evaluate the area after solving the boundary condition

After bisection, the sequence is rebuilt from the final \(x_1\). The code then forces \(x_{m+1}=1\) exactly to eliminate tiny floating-point boundary drift. With the final sequence in hand, it computes

$$Q=\sum_{i=0}^{m}(x_{i+1}-x_i)f(x_i),$$

and returns

$$A_N=4Q.$$

So the numerical work has two clean phases: first satisfy the optimality equations with the endpoint condition, then evaluate the resulting staircase area.

5. Checkpoints and numerical interpretation

For \(N=2\) we have \(m=1\). The optimum occurs at

$$x_1=\frac{1}{\sqrt{2}}, \qquad f(x_1)=\frac{1}{\sqrt{2}},$$

which gives

$$A_2=4\left(x_1+(1-x_1)f(x_1)\right)=4\sqrt{2}-2.$$

This is the first checkpoint used in the C++ program.

For \(N=10\), shooting produces the breakpoints

$$0,\ 0.3800728430,\ 0.5627007923,\ 0.7071067812,\ 0.8266606428,\ 0.9249565579,\ 1,$$

and the total area

$$A_{10}\approx 3.3469640797.$$

This is the second checkpoint in the C++ file. With the default input \(N=400\), the implementations output

$$A_{400}\approx 3.1486734435.$$

How the Code Works

f(x) evaluates the quarter-circle height. build_sequence stores the breakpoints and applies the recurrence step by step. shooting_residual returns \(R(x_1)\). solve_breakpoints performs the fixed-count bisection loop and then enforces the right endpoint numerically. Finally, minimal_red_area accumulates the strip areas and multiplies by \(4\). The Python and Java files are direct translations of the same method; the C++ version additionally checks \(N=2\) and \(N=10\) before printing the default \(N=400\) answer.

Complexity Analysis

Let \(m=N/2\). One call to build_sequence computes \(m+1\) new values, so it costs \(O(m)\) time and \(O(m)\) memory. Bisection performs a fixed 260 residual evaluations, hence the implemented runtime is \(O(260m)=O(N)\). If the iteration count is written symbolically as \(B\), the method is \(O(BN)\) time and \(O(N)\) memory.

References

  1. Problem page: https://projecteuler.net/problem=392
  2. Unit circle: Wikipedia — Unit circle
  3. Bisection method: Wikipedia — Bisection method
  4. Shooting method: Wikipedia — Shooting method
  5. Riemann sum: Wikipedia — Riemann sum

Problem 392 source code

C++

#include <cmath>
#include <iomanip>
#include <iostream>
#include <limits>
#include <string>
#include <vector>

namespace {

struct Options {
    int n = 400;
    bool run_checkpoints = true;
};

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 ch : tail) {
        if (ch < '0' || ch > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<int>(ch - '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_int_after_prefix(arg, "--n=", options.n)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.n >= 2 && (options.n % 2 == 0);
}

long double f(const long double x) {
    const long double inside = 1.0L - x * x;
    return std::sqrt(std::max(0.0L, inside));
}

bool build_sequence(const int m, const long double x1, std::vector<long double>& x) {
    x.assign(static_cast<std::size_t>(m + 2), 0.0L);
    x[0] = 0.0L;
    x[1] = x1;

    if (!(x1 > 0.0L && x1 < 1.0L)) {
        return false;
    }

    for (int k = 1; k <= m; ++k) {
        const long double x_prev = x[static_cast<std::size_t>(k - 1)];
        const long double x_cur = x[static_cast<std::size_t>(k)];
        const long double y_prev = f(x_prev);
        const long double y_cur = f(x_cur);

        if (x_cur <= 0.0L || y_cur <= 0.0L) {
            return false;
        }

        const long double x_next = x_cur + (y_prev - y_cur) * y_cur / x_cur;
        if (!std::isfinite(static_cast<double>(x_next))) {
            return false;
        }
        x[static_cast<std::size_t>(k + 1)] = x_next;
    }
    return true;
}

long double shooting_residual(const int m, const long double x1) {
    std::vector<long double> x;
    if (!build_sequence(m, x1, x)) {
        return std::numeric_limits<long double>::infinity();
    }
    return x[static_cast<std::size_t>(m + 1)] - 1.0L;
}

std::vector<long double> solve_breakpoints(const int n) {
    const int m = n / 2;
    long double lo = 1e-18L;
    long double hi = 0.999999999999L;

    for (int iter = 0; iter < 260; ++iter) {
        const long double mid = (lo + hi) / 2.0L;
        const long double r = shooting_residual(m, mid);
        if (r > 0.0L) {
            hi = mid;
        } else {
            lo = mid;
        }
    }

    const long double x1 = (lo + hi) / 2.0L;
    std::vector<long double> x;
    build_sequence(m, x1, x);
    x[static_cast<std::size_t>(m + 1)] = 1.0L;  // enforce boundary exactly after shooting.
    return x;
}

long double minimal_red_area(const int n) {
    const int m = n / 2;
    const std::vector<long double> x = solve_breakpoints(n);

    long double quadrant_area = 0.0L;
    for (int i = 0; i <= m; ++i) {
        const long double width = x[static_cast<std::size_t>(i + 1)] - x[static_cast<std::size_t>(i)];
        quadrant_area += width * f(x[static_cast<std::size_t>(i)]);
    }

    return 4.0L * quadrant_area;
}

bool nearly_equal(long double a, long double b, long double tol = 1e-13L) {
    const long double scale = std::max(std::fabsl(a), std::fabsl(b));
    if (scale == 0.0L) {
        return true;
    }
    return std::fabsl(a - b) <= tol * scale;
}

bool run_checkpoints() {
    const long double area_n2 = minimal_red_area(2);
    const long double expected_n2 = 4.0L * std::sqrt(2.0L) - 2.0L;
    if (!nearly_equal(area_n2, expected_n2, 1e-12L)) {
        std::cerr << "Checkpoint failed for N=2\n";
        return false;
    }

    const long double area_n10 = minimal_red_area(10);
    const long double expected_n10 = 3.3469640797L;
    if (!nearly_equal(area_n10, expected_n10, 2e-11L)) {
        std::cerr << "Checkpoint failed for N=10 (got " << std::setprecision(16)
                  << static_cast<double>(area_n10) << ")\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 = minimal_red_area(options.n);
    std::cout << std::fixed << std::setprecision(10) << static_cast<double>(answer) << '\n';
    return 0;
}

Python

import math

def f(x):
    inside = 1.0 - x * x
    return math.sqrt(max(0.0, inside))

def build_sequence(m, x1):
    x = [0.0] * (m + 2)
    x[0] = 0.0
    x[1] = x1

    if not (0.0 < x1 < 1.0):
        return None

    for k in range(1, m + 1):
        x_prev = x[k - 1]
        x_cur = x[k]
        y_prev = f(x_prev)
        y_cur = f(x_cur)

        if x_cur <= 0.0 or y_cur <= 0.0:
            return None

        x_next = x_cur + (y_prev - y_cur) * y_cur / x_cur
        if not math.isfinite(x_next):
            return None
        x[k + 1] = x_next
    return x

def shooting_residual(m, x1):
    x = build_sequence(m, x1)
    if x is None:
        return float('inf')
    return x[m + 1] - 1.0

def solve_breakpoints(n):
    m = n // 2
    lo = 1e-18
    hi = 0.999999999999

    for _ in range(260):
        mid = (lo + hi) / 2.0
        r = shooting_residual(m, mid)
        if r > 0.0:
            hi = mid
        else:
            lo = mid

    x1 = (lo + hi) / 2.0
    x = build_sequence(m, x1)
    x[m + 1] = 1.0
    return x

def minimal_red_area(n):
    m = n // 2
    x = solve_breakpoints(n)

    quadrant_area = 0.0
    for i in range(m + 1):
        width = x[i + 1] - x[i]
        quadrant_area += width * f(x[i])

    return 4.0 * quadrant_area

def solve():
    n = 400
    answer = minimal_red_area(n)
    return "{:.10f}".format(answer)

if __name__ == '__main__':
    print(solve())

Java

public class Euler392 {
    private static double f(double x) {
        double inside = 1.0 - x * x;
        return Math.sqrt(Math.max(0.0, inside));
    }

    private static double[] buildSequence(int m, double x1) {
        double[] x = new double[m + 2];
        x[0] = 0.0;
        x[1] = x1;

        if (!(x1 > 0.0 && x1 < 1.0)) {
            return null;
        }

        for (int k = 1; k <= m; ++k) {
            double xPrev = x[k - 1];
            double xCur = x[k];
            double yPrev = f(xPrev);
            double yCur = f(xCur);

            if (xCur <= 0.0 || yCur <= 0.0) {
                return null;
            }

            double xNext = xCur + (yPrev - yCur) * yCur / xCur;
            if (!Double.isFinite(xNext)) {
                return null;
            }
            x[k + 1] = xNext;
        }
        return x;
    }

    private static double shootingResidual(int m, double x1) {
        double[] x = buildSequence(m, x1);
        if (x == null) {
            return Double.POSITIVE_INFINITY;
        }
        return x[m + 1] - 1.0;
    }

    private static double[] solveBreakpoints(int n) {
        int m = n / 2;
        double lo = 1e-18;
        double hi = 0.999999999999;

        for (int iter = 0; iter < 260; ++iter) {
            double mid = (lo + hi) / 2.0;
            double r = shootingResidual(m, mid);
            if (r > 0.0) {
                hi = mid;
            } else {
                lo = mid;
            }
        }

        double x1 = (lo + hi) / 2.0;
        double[] x = buildSequence(m, x1);
        x[m + 1] = 1.0;
        return x;
    }

    private static double minimalRedArea(int n) {
        int m = n / 2;
        double[] x = solveBreakpoints(n);

        double quadrantArea = 0.0;
        for (int i = 0; i <= m; ++i) {
            double width = x[i + 1] - x[i];
            quadrantArea += width * f(x[i]);
        }

        return 4.0 * quadrantArea;
    }

    public static String solve() {
        int n = 400;
        double answer = minimalRedArea(n);
        return String.format("%.10f", answer).replace(',', '.');
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}