Detailed mathematical approach

Problem Summary

The implementations evaluate an expected median area arising from a random rectangle construction. After normalizing two complementary horizontal gaps and two complementary vertical gaps, the geometry depends only on two ratios \(x,y\in[0,1]\). The target value becomes a two-dimensional integral of a normalized kernel, followed by a constant rescaling by \(1/4\).

Mathematical Approach

Introduce the complementary fractions

$x_0=1-x,\qquad y_0=1-y.$

The normalized geometry is completely determined by these four numbers. Each admissible relative arrangement contributes the median of three candidate areas, and the kernel is the average over six equally weighted arrangements.

Step 1: Express the Median of Three Areas

For any triple \((a,b,c)\), the median can be written without sorting as

$\operatorname{med}(a,b,c)=a+b+c-\min(a,b,c)-\max(a,b,c).$

This identity is exactly what the implementation uses. It avoids a separate ordering step and gives the middle value directly.

Step 2: Build the Six Normalized Configurations

From the normalized side fractions, the six median terms are

$\begin{aligned} T_1(x,y)&=\operatorname{med}(xy,\ x_0y_0,\ 1),\\ T_2(x,y)&=\operatorname{med}(xy,\ x_0,\ y_0),\\ T_3(x,y)&=\operatorname{med}(xy_0,\ x_0y,\ 1),\\ T_4(x,y)&=\operatorname{med}(xy_0,\ x_0,\ y),\\ T_5(x,y)&=\operatorname{med}(x,\ x_0y,\ y_0),\\ T_6(x,y)&=\operatorname{med}(x,\ x_0y_0,\ y). \end{aligned}$

The normalized kernel is their average:

$M(x,y)=\frac{T_1(x,y)+T_2(x,y)+T_3(x,y)+T_4(x,y)+T_5(x,y)+T_6(x,y)}{6}.$

This formula is the mathematical core of the solution: once \(M(x,y)\) is available, the rest of the task is numerical integration.

Step 3: Exploit the Symmetries

Replacing \(x\) with \(1-x\), replacing \(y\) with \(1-y\), or swapping \(x\) and \(y\) only permutes the six triples above. Therefore

$M(x,y)=M(1-x,y)=M(x,1-y)=M(y,x).$

These identities are important for two reasons. First, they confirm that the six-term average has the expected geometric symmetry. Second, they justify computing the double quadrature with symmetry in the node indices.

Step 4: Turn the Expectation into an Integral

The normalized mean contribution is

$J=\int_0^1\int_0^1 M(x,y)\,dx\,dy.$

After the geometric normalization, the original expected area is obtained by multiplying by the scale factor left outside the ratio variables. In the implemented derivation that factor is \(1/4\), so the final quantity is

$E=\frac{J}{4}.$

The entire computational problem is therefore reduced to evaluating \(J\) very accurately.

Step 5: Approximate the Integral with Gauss-Legendre Quadrature

Let \((\xi_i,w_i)_{i=1}^n\) be the \(n\)-point Gauss-Legendre nodes and weights on \([0,1]\). Then

$J_n=\sum_{i=1}^{n}\sum_{j=1}^{n} w_i w_j M(\xi_i,\xi_j)$

approximates \(J\). Because the kernel is symmetric in its two arguments, the implementations reorganize the sum as

$J_n=\sum_{i=1}^{n} w_i^2 M(\xi_i,\xi_i)+2\sum_{1\le i\lt j\le n} w_i w_j M(\xi_i,\xi_j).$

This keeps the same quadrature rule while cutting the off-diagonal work almost in half.

Step 6: Cancel the Leading Error with Richardson Extrapolation

The kernel is smooth inside regions separated by median-switching boundaries, so the global quadrature error behaves like a piecewise-smooth integral rather than a fully analytic one. The implementation models the leading term as

$E_n=E+\frac{C}{n^2}+O\left(\frac{1}{n^4}\right).$

Computing once with \(n\) nodes and once with \(2n\) nodes gives

$E_{2n}=E+\frac{C}{4n^2}+O\left(\frac{1}{n^4}\right).$

Eliminating the unknown constant \(C\) yields

$E_{\text{rich}}=E_{2n}+\frac{E_{2n}-E_n}{3}.$

This Richardson-refined value is the reported answer.

Worked Example: Evaluate the Kernel at \(x=0\), \(y=\frac{1}{4}\)

Here \(x_0=1\) and \(y_0=\frac{3}{4}\). The six median terms become

$\begin{aligned} T_1&=\operatorname{med}\left(0,\frac{3}{4},1\right)=\frac{3}{4},\\ T_2&=\operatorname{med}\left(0,1,\frac{3}{4}\right)=\frac{3}{4},\\ T_3&=\operatorname{med}\left(0,\frac{1}{4},1\right)=\frac{1}{4},\\ T_4&=\operatorname{med}\left(0,1,\frac{1}{4}\right)=\frac{1}{4},\\ T_5&=\operatorname{med}\left(0,\frac{1}{4},\frac{3}{4}\right)=\frac{1}{4},\\ T_6&=\operatorname{med}\left(0,\frac{3}{4},\frac{1}{4}\right)=\frac{1}{4}. \end{aligned}$

Therefore

$M\left(0,\frac{1}{4}\right)=\frac{\frac{3}{4}+\frac{3}{4}+\frac{1}{4}+\frac{1}{4}+\frac{1}{4}+\frac{1}{4}}{6}=\frac{5}{12},$

which is one of the checkpoint values verified by the implementation.

How the Code Works

The C++, Python, and Java implementations all follow the same numerical plan. They evaluate the six median expressions for any requested pair \((x,y)\), generate Gauss-Legendre nodes and weights on \([0,1]\) from Legendre roots, and accumulate the symmetric double sum for two resolutions, \(n\) and \(2n\).

Each integral approximation is divided by \(4\), then the two results are combined with Richardson extrapolation. The C++ and Java implementations parallelize the outer loop over quadrature rows, while the Python implementation delegates to the compiled numerical solver and returns the parsed final value.

Before reporting the answer, the workflow checks fixed kernel values, symmetry identities, and benchmark quadrature outputs at moderate orders. Those checks guard against both algebraic mistakes in the kernel and numerical mistakes in the quadrature generator.

Complexity Analysis

For order \(n\), generating the Gauss-Legendre nodes and weights costs \(O(n^2)\) arithmetic overall because each root evaluation uses a degree-\(n\) recurrence and there are \(O(n)\) roots. The dominant cost is the double quadrature itself, which requires \(O(n^2)\) kernel evaluations and \(O(n)\) memory for the node and weight arrays. Because the sum is split by rows, the wall-clock time benefits directly from multi-core parallelism.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=695
2. Gaussian quadrature: Wikipedia — Gaussian quadrature
3. Legendre polynomials: Wikipedia — Legendre polynomials
4. Richardson extrapolation: Wikipedia — Richardson extrapolation
5. Median: Wikipedia — Median