Detailed mathematical approach

Problem Summary

Fix an outer square \(O=[0,n]\times[0,n]\). A valid hollow lamina is obtained by removing an interior axis-aligned rectangle

$H=[a,a+x]\times[b,b+y],\qquad 1\le x,y\le n-2,\qquad 1\le a\le n-x-1,\qquad 1\le b\le n-y-1.$

The remaining region is \(R=O\setminus H\), whose area is \(A=n^2-xy\). If two points are chosen independently and uniformly from \(R\), their expected distance is

$E(R)=\frac{1}{A^2}\int_R\int_R \|p-q\|\,dp\,dq.$

The quantity \(S(n)\) is the sum of \(E(R)\) over all valid choices of \(x,y,a,b\). The program computes this total and rounds the final value to four decimal places.

Mathematical Approach

The implementation converts the continuous geometry into a finite collection of precomputed interaction tables. The key object is a distance kernel for two unit cells, and everything else is built from that kernel by counting how many cell pairs realize each offset.

Step 1: Build the Unit-Cell Distance Kernel

Partition every rectangle into unit cells. For nonnegative integer offsets \((d_x,d_y)\), define \(K(d_x,d_y)\) as the total distance integral for two random points drawn from two unit cells whose lower-left corners differ by \((d_x,d_y)\):

$K(d_x,d_y)=\int_0^1\int_0^1\int_0^1\int_0^1 \sqrt{(d_x+u_1-u_2)^2+(d_y+v_1-v_2)^2}\,du_1\,du_2\,dv_1\,dv_2.$

After introducing the differences \(s=u_1-u_2\) and \(t=v_1-v_2\), the density becomes triangular, so the same quantity can be written as

$K(d_x,d_y)=\int_{-1}^{1}\int_{-1}^{1}(1-|s|)(1-|t|)\sqrt{(d_x+s)^2+(d_y+t)^2}\,ds\,dt.$

The implementations exploit symmetry and evaluate this integral numerically on \([0,1]^2\) with Gauss-Legendre quadrature:

$\begin{aligned} K(d_x,d_y)=\int_0^1\int_0^1 &(1-u)(1-v)\Bigl( \sqrt{(d_x+u)^2+(d_y+v)^2} +\sqrt{(d_x+u)^2+(d_y-v)^2}\\ &+\sqrt{(d_x-u)^2+(d_y+v)^2} +\sqrt{(d_x-u)^2+(d_y-v)^2} \Bigr)\,du\,dv. \end{aligned}$

A useful anchor is the special case \(K(0,0)\), the classical mean distance between two random points in the unit square:

$K(0,0)=\frac{2+\sqrt{2}+5\ln(1+\sqrt{2})}{15}.$

Step 2: Aggregate a Whole Rectangle from Cell Offsets

Let \(F(w,h)\) denote the total double integral of distance over a full \(w\times h\) rectangle:

$F(w,h)=\int_{[0,w]\times[0,h]}\int_{[0,w]\times[0,h]} \|p-q\|\,dp\,dq.$

If two unit cells have relative offset \((d_x,d_y)\), then there are exactly \((w-|d_x|)(h-|d_y|)\) ordered pairs of cells with that offset. Therefore

$F(w,h)=\sum_{d_x=-(w-1)}^{w-1}\sum_{d_y=-(h-1)}^{h-1}(w-|d_x|)(h-|d_y|)\,K(|d_x|,|d_y|).$

This formula is the backbone of the precomputed rectangle table. It gives the full distance integral, not yet the expected distance. To get the mean distance inside a solid rectangle alone, one would divide by \((wh)^2\).

Step 3: Use Inclusion-Exclusion for the Hollow Region

For any two regions \(A\) and \(B\), write

$F_{AB}=\int_A\int_B \|p-q\|\,dp\,dq.$

Since the lamina is \(R=O\setminus H\), the indicator identity \(1_R=1_O-1_H\) implies

$F_{RR}=F_{OO}-F_{OH}-F_{HO}+F_{HH}.$

Distance is symmetric, so \(F_{OH}=F_{HO}\). Hence the usable formula is

$F_{RR}=F_{OO}-2F_{OH}+F_{HH}.$

Here \(F_{OO}=F(n,n)\) comes from the outer square table, \(F_{HH}=F(x,y)\) comes from the hole-size table, and only the cross term \(F_{OH}\) depends on the placement \((a,b)\). Once \(F_{RR}\) is known, the expected distance for that lamina is

$E(R)=\frac{F_{RR}}{(n^2-xy)^2}.$

Step 4: Compute the Outer-Hole Cross Term with Prefix Sums

For each unit cell \((i,j)\) inside the \(n\times n\) outer square, define its interaction with the whole outer square by

$W_n(i,j)=\sum_{u=0}^{n-1}\sum_{v=0}^{n-1} K(|u-i|,|v-j|).$

If the hole occupies the block of cells

$a\le i\le a+x-1,\qquad b\le j\le b+y-1,$

then the cross term is simply the sum of \(W_n(i,j)\) over that block:

$F_{OH}=\sum_{i=a}^{a+x-1}\sum_{j=b}^{b+y-1} W_n(i,j).$

A two-dimensional prefix sum over the \(W_n\) table reduces each such rectangular query to \(O(1)\) time, which is essential because the same outer square must be combined with many hole sizes and many placements.

Step 5: Sum Over All Valid Holes

Combining the previous steps yields

$S(n)=\sum_{x=1}^{n-2}\sum_{y=1}^{n-2}\sum_{a=1}^{n-x-1}\sum_{b=1}^{n-y-1} \frac{F(n,n)-2F_{OH}(n;x,y,a,b)+F(x,y)}{(n^2-xy)^2}.$

This formula is exactly what the implementations evaluate. The only numerical approximation occurs in the kernel table \(K(d_x,d_y)\), and the rest is deterministic table lookup, prefix-sum extraction, and accumulation.

Worked Example: \(n=3\) and \(n=4\)

When \(n=3\), there is only one valid hole: \(x=y=1\) placed at \(a=b=1\). The lamina area is \(3^2-1=8\), so

$S(3)=\frac{F(3,3)-2F_{OH}+F(1,1)}{8^2}\approx 1.6514.$

This matches the checkpoint used by the implementations.

For \(n=4\), the number of admissible hollow laminae is

$\sum_{x=1}^{2}\sum_{y=1}^{2}(4-x-1)(4-y-1)=9.$

That count is another checkpoint: there are \(4\) placements for a \(1\times 1\) hole, \(2\) for \(1\times 2\), \(2\) for \(2\times 1\), and \(1\) for \(2\times 2\).

How the Code Works

The C++, Python, and Java implementations all follow the same pipeline. First they generate Gauss-Legendre nodes and weights on \([0,1]\) from Legendre polynomial roots and derivatives. With those nodes they fill the symmetric kernel table \(K(d_x,d_y)\) for all offsets \(0\le d_x,d_y\le n-1\).

Next they build a table of rectangle self-interactions \(F(w,h)\) for every \(1\le w,h\le n\). For the requested \(n\), they then compute the outer-to-cell interaction grid \(W_n(i,j)\) and turn it into a two-dimensional prefix sum so that each placement-dependent cross term \(F_{OH}\) is recovered by one rectangle query.

Finally they iterate over every valid hole size and every legal position, combine the three terms \(F_{OO}\), \(F_{OH}\), and \(F_{HH}\) by inclusion-exclusion, divide by the squared lamina area, and add the result to the running total. The C++ implementation also distributes independent hole-size pairs across worker threads, while preserving the same mathematics as the Python and Java implementations.

Complexity Analysis

Let \(q\) be the Gauss-Legendre order. Building the kernel table requires \(O(n^2q^2)\) time and \(O(n^2)\) memory. Building the rectangle self-interaction table costs \(O(n^4)\) time, and forming the outer interaction grid together with its prefix sum also costs \(O(n^4)\) time. The final sweep over all hole sizes and placements is

$\sum_{x=1}^{n-2}\sum_{y=1}^{n-2}(n-x-1)(n-y-1)=O(n^4),$

with \(O(1)\) work per placement after preprocessing. Therefore, for fixed quadrature order, the overall algorithm runs in \(O(n^4)\) time and uses \(O(n^2)\) memory.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=547
2. Gauss-Legendre quadrature: Wikipedia — Gauss-Legendre quadrature
3. Inclusion-exclusion principle: Wikipedia — Inclusion-exclusion principle
4. Two-dimensional prefix sums: Wikipedia — Summed-area table
5. The exact value \(K(0,0)=\frac{2+\sqrt{2}+5\ln(1+\sqrt{2})}{15}\) is the classical mean distance between two random points in a unit square.