Problem 392: Enmeshed Unit Circle

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

A rectilinear grid is an orthogonal grid whose gridlines need not be evenly spaced; logarithmic graph paper is the everyday example. ^[3] Take such a grid with $N+2$ vertical and $N+2$ horizontal lines, hence $(N+1)\times(N+1)$ cells, the outermost being $x=\pm1$ and $y=\pm1$; the $N$ inner lines of each direction are free. A cell is coloured red if it overlaps the unit circle $x^2+y^2=1$ ^[2] and black otherwise, and the inner lines must be placed so that the total red area is least. The statement gives that minimum for $N=10$ as $3.3469640797$ and asks for $N=400$, again to ten decimal places. ^[1]

Overlapping means sharing area: a cell is red exactly when its interior meets the open disc $x^2+y^2\lt1$, so a cell touching the circle at a single corner is black. That is not pedantry — in the best grid many corners land exactly on the circle, and counting those cells red would give $3.6273329731$ for the same $N=10$ grid, not the stated value.

The whole solution in one sentence: the best grid is symmetric and tuned so that the staircase of red cells meets the circle exactly at its corners, which gives one equation per inner line, and the whole chain of equations follows from a single unknown, found by bisection.

Mathematical Approach

1. When is a cell red?

Write $d(I)$ for the distance from $0$ to an interval $I$ on a coordinate axis: $0$ if $I$ contains $0$, otherwise the smaller endpoint distance. A cell is a product $I\times J$ of a vertical and a horizontal interval, and the point of that rectangle nearest the origin has coordinates $\pm d(I)$ and $\pm d(J)$, so $I\times J$ is red if and only if

$$d(I)^2+d(J)^2\lt1 .$$

With inner lines at $\pm\tfrac12$ the three vertical intervals have $d=\tfrac12,0,\tfrac12$, so a corner cell has $\tfrac14+\tfrac14\lt1$ and is red; at $\pm0.8$ it has $0.64+0.64\gt1$ and turns black.

2. A lower bound for the red area

Number the vertical intervals $i$, with widths $W_i$ and distances $D_i$, and the horizontal ones $j$, with $H_j$ and $E_j$. Grouped by column, the red area is

$$A=\sum_i W_i\,T_i,\qquad T_i=\sum_{j\,:\,E_j\lt f(D_i)}H_j,$$

where $f(x)=\sqrt{1-x^2}$ is the upper half of the circle, so $D_i^2+E_j^2\lt1$ says $E_j\lt f(D_i)$. Put $c=f(D_i)$. A horizontal interval $J$ has $d(J)\lt c$ exactly when $J$ meets $(-c,c)$, so the intervals counted in $T_i$ are consecutive and their union is an interval containing $(-c,c)$; its length is at least $2c$. Hence

$$A\ \ge\ 2\sum_i W_i\,f(D_i),$$

with equality precisely when $+f(D_i)$ and $-f(D_i)$ are horizontal gridlines for every column $i$. The plan is to minimise the right-hand side, then show the minimiser meets that equality condition.

3. Splitting at the origin: two staircases

Let $0\lt b_1\lt\cdots\lt b_R=1$ be the positive vertical lines and $0\lt c_1\lt\cdots\lt c_L=1$ the absolute values of the negative ones. A vertical interval right of the origin is some $[b_k,b_{k+1}]$ with $D=b_k$, while the interval containing the origin has $D=0$, $f(0)=1$ and width $b_1+c_1$, so it splits cleanly between the sides. Hence

$$\sum_i W_i\,f(D_i)=G(b)+G(c),$$

$$G(t)=t_1+\sum_{k=1}^{R-1}(t_{k+1}-t_k)\,f(t_k).$$

The same holds with a gridline at the origin, since then $[0,b_1]$ has $D=0$ and width $b_1$. So the sum splits into two one-sided problems of the same shape. Write $g(j)$ for the least value of $G$ over chains with $j$ free lines in $(0,1)$ plus the fixed endpoint $1$. Of the $N+2$ vertical lines, $R+L=N+2$ when none is at the origin and $R+L=N+1$ when one is, and the chains have $R-1$ and $L-1$ free lines, whence

$$A\ \ge\ 2\bigl(g(j_R)+g(j_L)\bigr),\qquad j_R+j_L\le N,$$

with $j_R+j_L=N$ only if no gridline passes through the origin.

4. Why the split must be balanced

Adding one more line strictly lowers $G$: inserting $t$ inside a step $[t_k,t_{k+1}]$ replaces $(t_{k+1}-t_k)f(t_k)$ by $(t-t_k)f(t_k)+(t_{k+1}-t)f(t)$, and $f$ is strictly decreasing, so $f(t)\lt f(t_k)$ and the new pair is smaller. Hence $g$ strictly decreases: $g(0)=1$, $g(1)=\sqrt2-\tfrac12=0.9142135624$, $g(2)=0.8784148931$, $g(5)=0.8367410199$. Write $\Delta_k=g(k+1)-g(k)$ for the increments. Computing $g(m)$ for every $m\le400$ shows that it is also strictly convex: the $\Delta_k$ are negative and strictly increasing. That is what makes the balanced split best, for moving a line from a side with $b$ free lines to one with $a\le b-2$ changes the bound by $\Delta_a-\Delta_{b-1}\le0$, and repeating the move balances the split without ever raising the bound. This convexity is a computed fact, not proved here; the rest of the section is proved. For even $N$,

$$A\ \ge\ 4\,g(N/2).$$

The balanced split needs $N$ even, since $j_R=j_L$ with $R+L=N+2$ makes $N+2=2R$; for odd $N$ a symmetric grid must put a line at the origin, which costs one free line, so the bound is weaker and the construction below does not apply. The method therefore targets even $N$, and $400$ is even — though only the C++ command line enforces that.

5. Minimising one staircase

Fix $m=N/2$ free lines on one side and rename them $x_1\lt\cdots\lt x_m$, with the conventions $x_0=0$ and $x_{m+1}=1$. Because $f(0)=1$, the functional $G$ becomes

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

the area under a staircase whose step over $[x_i,x_{i+1}]$ has height $f(x_i)$. $Q$ is continuous on the compact region $0\le x_1\le\cdots\le x_m\le1$, so a minimum exists, and it cannot lie on the boundary: if two values coincide, or $x_1=0$, or $x_m=1$, the configuration uses fewer than $m$ distinct lines and the insertion argument of Section 4 lowers $Q$ strictly. The minimum is therefore interior, where every partial derivative vanishes. ^[4] Only two summands contain $x_k$, namely $(x_k-x_{k-1})f(x_{k-1})$ and $(x_{k+1}-x_k)f(x_k)$, so

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

Differentiating $f$ gives $f'(x)=-x/f(x)$, and setting the derivative to zero and solving for $x_{k+1}$ turns the $m$ equations into a recurrence:

$$x_{k+1}=x_k+\frac{\bigl(f(x_{k-1})-f(x_k)\bigr)f(x_k)}{x_k}.$$

Geometrically: at the optimum each step's loss of height $f(x_{k-1})-f(x_k)$ balances the width of the next step.

6. One unknown, found by bisection

Since $x_0=0$ is fixed, the recurrence makes the whole configuration a function of the single number $x_1$: run it for $k=1,\dots,m$ and read off $x_{m+1}$. Every difference it produces is positive while the values stay below $1$, because $f(x_{k-1})-f(x_k)\gt0$ whenever $x_{k-1}\lt x_k$, so a sequence started with $x_1\gt0$ increases by itself and can only fail by reaching $1$ too early. Call $x_1$ usable when every $x_k$ with $k\le m$ stays strictly inside $(0,1)$; the usable values form an open interval $(0,\tau_m)$, equal to $(0,1/\sqrt2)$ for $m=2$ — the endpoint excluded because it makes $x_2=1$ — and to about $(0,0.040124)$ for $m=200$. What remains is one boundary condition,

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

Reducing a boundary value problem to repeated guesses of one initial value in this way is the shooting method. ^[7] Searching by the sign of $R$ is then a bisection: the loop keeps $R\le0$ at the lower end of the bracket and $R\gt0$ at the upper, scoring unusable trials $+\infty$, which counts as positive and pushes the bracket down. ^[8] Two warnings. $R$ is not monotone: at $m=2$ it climbs to about $0.0939$ near $x_1=0.674$ before falling back. And $R$ tends to $0$ at $\tau_m$ too, where the breakpoints pile up at $1$ — a degenerate configuration the guards reject and Section 5 excluded, having $x_m=1$. For every $m\le200$ tested, $R$ is strictly negative below the solution and strictly positive above, so the sign change is unique.

7. The reflection, and why the bound is attained

The map $f$ is an involution on $[0,1]$: $f(f(x))=x$. ^[6] That lets us reflect a configuration. Given $x_1\lt\cdots\lt x_m$, put

$$\sigma x=\bigl(f(x_m),\dots,f(x_1)\bigr),$$

again increasing in $(0,1)$; the claim is $Q(\sigma x)=Q(x)$. The staircase region $\mathcal R$, the union of the rectangles $[x_i,x_{i+1}]\times[0,f(x_i)]$, has area $Q(x)$. Reflecting in the line $y=x$ maps it to the union of $[0,f(x_i)]\times[x_i,x_{i+1}]$, and a short check identifies that set: a point with first coordinate between $f(x_{i+1})$ and $f(x_i)$ belongs to it for second coordinates up to $x_{i+1}=f(f(x_{i+1}))$. It is therefore the staircase region of $\sigma x$, and reflection preserves area. For instance, with $x=(0.25,0.5,0.9)$,

$$\begin{aligned}\sigma x&=\bigl(f(0.9),f(0.5),f(0.25)\bigr)\\&=(0.435890,\ 0.866025,\ 0.968246),\end{aligned}$$

and both give $Q=0.8820606101$.

By Section 5 the minimiser is an interior stationary point, hence a usable zero of $R$, and the single sign change leaves only one such zero, so the minimiser is unique. As $\sigma$ preserves $Q$, it permutes the minimisers; with only one it must fix it:

$$f(x_k)=x_{m+1-k}\qquad (k=0,1,\dots,m+1).$$

Every $f(x_k)$ is thus itself a breakpoint, so in the symmetric grid built from them $\pm f(D_i)$ is a gridline for every column — exactly the equality condition of Section 2. The bound is attained and the minimal red area is $4\,g(N/2)$, granted the two computed facts above: the convexity of $g$ and the single sign change of $R$. The grid is interlocked with the circle at every corner of the staircase, which is what the title refers to.

8. The case $N=2$ worked out by hand

Here $m=1$ and the fixed point decides everything: $f(x_1)=x_1$ forces $1-x_1^2=x_1^2$, so $x_1=1/\sqrt2=0.7071067812$. The grid has vertical lines $-1,-1/\sqrt2,1/\sqrt2,1$ and the same horizontally: $9$ cells, with interval distances $1/\sqrt2,0,1/\sqrt2$. The central cell has $0+0\lt1$ and the four edge cells $0+\tfrac12\lt1$, so all five are red; each corner cell has $\tfrac12+\tfrac12=1$ exactly, touches the circle at one point and is black. Adding the five red areas,

$$2+4\Bigl(1-\tfrac1{\sqrt2}\Bigr)\sqrt2=4\sqrt2-2=3.6568542495,$$

agreeing with $4Q=4(\sqrt2-\tfrac12)$. This is the first checkpoint in the C++ program.

9. Two more examples: $N=4$ and $N=10$

For $m=2$ the reflection gives $x_2=f(x_1)$ and the stationarity equation at $k=1$ reads $f(x_1)=x_1+(1-f(x_1))f(x_1)/x_1$. With $x=x_1$, $y=f(x)$: multiplying by $x$ and using $y^2=1-x^2$ gives $xy=2x^2+y-1$, that is $(1-x)y=1-2x^2$, and squaring,

$$x\bigl(5x^3-2x^2-4x+2\bigr)=0 .$$

On $\tfrac12\lt x\lt\tfrac35$ both sides of $(1-x)y=1-2x^2$ are positive, so squaring adds nothing there and $x_1$ is the root of

$$5x^3-2x^2-4x+2=0$$

lying between $\tfrac12$ and $\tfrac35$, namely $x_1=0.5673070757$, with $x_2=f(x_1)=0.8235063339$. Then

$$\begin{aligned}Q&=2x_1+1-x_1^2-2x_1f(x_1)\\&=0.8784148931,\end{aligned}$$

so the red area is $3.5136595725$. The pairing counts cells exactly too: a cell off both axes, at distances $x_k$ and $x_l$, is red when $x_l\lt f(x_k)=x_{m+1-k}$, that is when $k+l\le m$, and touches the circle when $k+l=m+1$. Adding the central cell, the $4m$ on the axes and the $4\cdot m(m-1)/2$ off them, the grid has $2m^2+2m+1$ red cells and $4m$ corner-only contacts among $(2m+1)^2$: for $N=4$, $13$ and $8$ out of $25$. For $N=10$ the five positive lines are

$$0.3800728430,\ 0.5627007923,\ 0.7071067812,$$

$$0.8266606428,\ 0.9249565579,$$

and the pairing is visible in them: $f(0.3800728430)\approx0.9249565579$ and $f(0.5627007923)\approx0.8266606428$ for these rounded values, while the middle one is the exact fixed point $1/\sqrt2$. Of the $121$ cells, $61$ are red and $20$ more touch the circle at a corner and are black; the red ones total $3.3469640797$, the value in the statement.

10. How close to the circle does the staircase get?

$Q$ is a left Riemann sum of $f$ for the partition $x_0\lt\cdots\lt x_{m+1}$ of $[0,1]$, and as $f$ decreases such a sum overestimates the integral. ^[5] Since $\int_0^1f=\pi/4$, always $Q\gt\pi/4$, and so every red area exceeds $\pi$. It does converge to $\pi$, as the uniform grid shows: with cell side $h=2/(N+1)$ every point of a red cell lies within $\sqrt2\,h$ of the disc, so the minimum is at most $\pi(1+\sqrt2\,h)^2\to\pi$. The optimal values fall faster: $3.3469640797$ at $N=10$, $3.1941919969$ at $N=50$, $3.1555861880$ at $N=200$.

The rate can be predicted: a step of width $h$ near $x$ leaves about $\tfrac12|f'(x)|h^2$ of excess, and spreading $m+1$ steps with density proportional to $\sqrt{|f'|}$ makes the total least, at

$$\begin{aligned}g(m)-\frac{\pi}{4}&\approx\frac{C}{m+1},\\C&=\tfrac12\Bigl(\int_0^1\!\!\sqrt{|f'(x)|}\,dx\Bigr)^{2}.\end{aligned}$$

Substituting $t=x^2$ makes it a Beta integral, ^[9]

$$\int_0^1\!\!\sqrt{\frac{x}{\sqrt{1-x^2}}}\,dx=\tfrac12\mathrm B\bigl(\tfrac34,\tfrac34\bigr),$$

which equals $\Gamma(3/4)^2/\sqrt\pi$, so $C=\Gamma(3/4)^4/(2\pi)=0.3588850055\ldots$ The step estimate is not uniform at $x=1$, where $f'$ is unbounded, so this predicts rather than proves; but the fit is clear, $(m+1)\bigl(g(m)-\pi/4\bigr)$ being $0.3259415$ at $m=10$, $0.3471194$ at $m=40$, $0.3551605$ at $m=160$ and $0.3577722$ at $m=640$.

How the Code Works

Only the C++ file has a command line: main calls parse_arguments to fill an Options record, $N=400$ by default, with --n= to override it and --skip-checkpoints to suppress the self-checks. For inputs that fit in an int it accepts only even $N\ge2$, rejecting the rest with exit status $1$ as Section 4 requires; a failed checkpoint gives status $2$. Python and Java have no options: solve() fixes $N=400$ and runs no checkpoints. In all three, f(x) clamps the radicand with max(0, inside), so a rounded value just below zero cannot give a NaN.

build_sequence(m, x1) implements Section 6: it allocates $m+2$ slots, sets $x_0=0$, and applies the recurrence of Section 5 for $k=1,\dots,m$. It refuses three cases — a starting value outside $(0,1)$, a step with $x_k\le0$ or $f(x_k)\le0$, where the sequence has run past the end of the quarter circle, and a non-finite $x_{k+1}$ — reporting false in C++, None in Python, null in Java. shooting_residual turns that into $+\infty$, and otherwise returns $x_{m+1}-1$.

solve_breakpoints(n) brackets $x_1$ in $[10^{-18},\,0.999999999999]$ and halves the bracket $260$ times, moving hi down when the residual is positive and lo up otherwise. Both ends behave as bisection needs. At the top the residual is positive: for $m\ge2$ the second step finds $x_2\gt1$ and rejects the trial, giving $+\infty$, while for $m=1$ it survives with a finite $1.41\cdot10^{-6}$. At the bottom $f(0)-f(10^{-18})$ rounds to $0$, so no step adds anything, the sequence stays at $10^{-18}$ and the residual is $\approx-1$. The sequence is then rebuilt from the final midpoint with $x_{m+1}$ overwritten by exactly $1$, and minimal_red_area(n) sums $(x_{i+1}-x_i)f(x_i)$ for $i=0,\dots,m$ and multiplies by $4$.

run_checkpoints compares $N=2$ against $4\sqrt2-2$ within relative $10^{-12}$ and $N=10$ against the literal $3.3469640797$ within $2\cdot10^{-11}$, nearly_equal scaling by the larger magnitude. The second is tighter than it looks: that literal is the statement's ten-digit rounding, $3.4\cdot10^{-11}$ from the computed value, leaving a factor of two of slack.

The languages differ in three ways. C++ stores breakpoints and area in long double, whose layout is implementation-defined: the best-known one, on many x86 implementations, has a 64-bit significand, while elsewhere — Apple Silicon included — long double is just double with $53$ bits. ^[10] The extra precision is not needed anyway, and the file converts to double inside the loop too, for the test isfinite(static_cast<double>(x_next)). Python and Java use $53$-bit doubles throughout, Java naming the helpers buildSequence, shootingResidual, solveBreakpoints and minimalRedArea and printing with String.format("%.10f", answer).replace(',', '.'), which repairs a decimal comma but not other digit shapes a locale may use.

Two rough edges, neither affecting the result. None of the three checks the final build_sequence call: on failure the C++ would sum a partly zero vector, Python raise a TypeError, Java a NullPointerException. It did succeed for every even $N$ from $2$ to $1000$ tested, the root sitting well inside the usable interval, about $0.03999$ against $\tau_m\approx0.040124$ at $m=200$. And parse_int_after_prefix accumulates digits into an int unguarded, so an argument past the int range is undefined behaviour, not a clean rejection: a sanitiser build reports a signed overflow for --n=2147483648.

Complexity and Verification

A successful residual evaluation performs $m$ recurrence steps, each a few operations and two square roots; a rejected one stops earlier, and allocating the $m+2$ values costs $O(N)$ either way. With a fixed $260$ evaluations the program is $O(N)$ in time and memory: at most $52\,000$ steps for $N=400$, plus the $m+1$ terms of the final sum, and all three finish in a few milliseconds. The count is generous — $2^{-260}$ is far beyond floating-point resolution, and the area stops changing in its last bit after some $50$ to $60$ halvings.

Why do $53$ bits of significand suffice for ten correct decimals? Two things help. The recurrence amplifies an error in $x_1$ only by a small constant, not exponentially: at $m=200$ the derivative of $x_{m+1}$ in $x_1$ is about $5.0$ and the largest over all breakpoints about $17.3$, near index $141$; at $m=10$ they are $1.9$ and $2.7$. And $Q$ is stationary at the optimum, so displacing the interior breakpoints by $\delta$, with $x_0=0$ and $x_{m+1}=1$ held fixed as the code holds them, changes $Q$ by $O(\delta^2)$: at $m=200$ a relative perturbation of $10^{-4}$ in $x_1$ moves $Q$ by $8.7\cdot10^{-9}$ relative, and $10^{-6}$ by $8.7\cdot10^{-13}$. Neither argument bounds the rounding of the square roots and the summation, so the real evidence is a recomputation at $70$ digits: the $53$-bit result is off by $6.7\cdot10^{-17}$ at $N=10$ and $9.2\cdot10^{-16}$ at $N=400$, a sixth of a unit and about two units in the last place — four orders of magnitude inside the $5\cdot10^{-11}$ that ten decimals need.

Beyond the two packaged checkpoints, which hold, this article verified the following.

At their default $N=400$ all three print byte-identical output on this machine, whose locale uses ASCII digits and a point.
Enumerating all $(N+1)^2$ cells of the optimal grid reproduces $4Q$ and the counts of Section 9 for $N=2,4,6,10,20,40$ — provided the corner contacts are recognised as such. Testing $D_i^2+E_j^2\lt1$ in plain double arithmetic on rounded breakpoints lets rounding decide them instead, giving $73$ red cells rather than $61$ at $N=10$ and overstating the area by $0.1668127428$.
The bound of Section 2 held on $2000$ random grids, with equality for the optimal symmetric ones, and $Q(\sigma x)=Q(x)$ held to $5.6\cdot10^{-16}$ on $200\,000$ random configurations.
Nelder–Mead then Powell, from six starts per $m$, reproduces the shooting value of $g(m)$ to $2\cdot10^{-16}$ for every $m\le8$, and the numerical Hessian of $Q$ there is positive definite (smallest eigenvalue $0.571$ at $m=2$, $0.0798$ at $m=8$). This tests the candidate; it is not itself a proof of global optimality.
The computed $g$ is strictly decreasing and strictly convex over $0\le m\le400$, and among the splits $j_R+j_L=400$ the balanced one uniquely minimises $g(j_R)+g(j_L)$, $199+201$ being worse. These are values of the shooting candidate, not a proof about the exact minima.

Footnotes and References

The references supply the background facts used above. The lower bound, the reduction to two staircases, the reflection identity and the asymptotic constant are derived in this article.

Project Euler 392 — Enmeshed Unit Circle. The official statement defines a rectilinear grid with N+2 gridlines in each direction, outermost at x = ±1 and y = ±1, colours a cell red if it overlaps the unit circle, asks where to put the N inner lines so that the red area is least, gives that minimum for N = 10 as 3.3469640797, and requests the value for N = 400 to ten decimal places.
Wikipedia — Unit circle. The unit circle has radius 1 and, centred at the origin as here, is the set of points with x² + y² = 1 — the equation Section 1 tests.
Wikipedia — Logarithmic scale. On a logarithmic scale "numbers with equal distance between them such as 1, 2, 3, 4, 5 are not equally spaced", and logarithmic graph paper was once a standard tool: the statement's example of uneven spacing.
Wikipedia — Fermat's theorem (stationary points). If a differentiable function has a local extremum at an interior point, its derivative vanishes there; the condition is necessary, not sufficient. Section 5 applies it to each variable of Q.
Wikipedia — Riemann sum. Riemann sums are defined for arbitrary partitions, widths need not be equal, and "the left Riemann sum amounts to an overestimation if f is monotonically decreasing on this interval" — Section 10's lower bound.
Wikipedia — Involution (mathematics). An involution is its own inverse, f(f(x)) = x for every x in its domain, and every involution is a bijection; reflection is listed as a geometric example. Section 7 uses both facts.
Wikipedia — Shooting method. The method solves a boundary value problem by "finding solutions to the initial value problem for different initial conditions", adjusting the unknown initial data until the far condition holds — Section 6, with x₁ as the unknown.
Wikipedia — Bisection method. The method "applies to any continuous function for which one knows two values with opposite signs", such a pair bracketing a root by the intermediate value theorem, and "an interval that contains a zero of f is reduced in width by 50% at each step" — the sign condition and halving rate used here.
NIST Digital Library of Mathematical Functions — §5.12, equation 5.12.1. The integral of t^a−1(1 − t)^b−1 over [0, 1] equals B(a, b) = Γ(a)Γ(b)/Γ(a + b) for positive a and b. With a = b = 3/4 and Γ(3/2) = √π/2 this evaluates the integral in Section 10.
cppreference — Fundamental types describes double as usually IEEE-754 binary64 and long double as an implementation-dependent extended type whose best-known layout, on many x86 and x86-64 implementations, is the 80-bit x87 format. Java — Double fixes double at 64 bits with a 53-bit significand, and Python — Floating-point arithmetic notes that "almost all platforms map Python floats to IEEE 754 binary64", with the same 53 bits.

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