Solution

Answer: 11730879.0023

Let

$$E(R)=\frac1{\operatorname{area}(R)^2}\iint_{R\times R} |P-Q|,dP,dQ$$

be the expected distance between two uniformly random points in a region $R$.

For a hollow square lamina, the region is

$$L = O\setminus H,$$

where

$O$ is the outer $n\times n$ square,
$H$ is the removed $x\times y$ rectangle.

We must compute

$$S(n)=\sum_{\text{all laminae}} E(L).$$

For $n=40$, brute force geometric integration is impossible unless we derive a closed-form rectangle integral.

1. Mathematical analysis

Step 1 — Rectangle pair integral

Define

$$I(A,B)=\iint_{A\times B}\sqrt{(x_1-x_2)^2+(y_1-y_2)^2},dA,dB.$$

Then

$$E(R)=\frac{I(R,R)}{\operatorname{area}(R)^2}.$$

So for a lamina $L=O\setminus H$,

$$I(L,L) = I(O,O)-2I(O,H)+I(H,H).$$

Thus everything reduces to evaluating $I$ for axis-aligned rectangles.

Step 2 — Antiderivative

Let

$$r=\sqrt{x^2+y^2}.$$

A useful identity is

$$\frac{\partial^2}{\partial x,\partial y} \left( \frac{xyr}{3} +\frac{x^3}{6}\operatorname{asinh}\frac{y}{x} +\frac{y^3}{6}\operatorname{asinh}\frac{x}{y} \right) =r.$$

Integrating once more in both variables gives a fourth-order primitive:

$$Q(x,y)= \frac{x^4y}{24}\operatorname{asinh}\frac{y}{x} +\frac{xy^4}{24}\operatorname{asinh}\frac{x}{y} -\frac{x^4r}{60} +\frac{x^2y^2r}{20} -\frac{y^4r}{60}.$$

Then by repeated inclusion–exclusion,

for rectangles

$$R_1=[a,b]\times[c,d], \qquad R_2=[e,f]\times[g,h],$$

we obtain

$$I(R_1,R_2) = \sum (-1)^{i+j+k+l} Q(x_i-u_j,; y_k-v_l),$$

where

$$x_i\in{a,b}, \quad u_j\in{e,f}, \quad y_k\in{c,d}, \quad v_l\in{g,h}.$$

This gives an exact closed form with only 16 evaluations.

Step 3 — Enumerating all laminae

For fixed $n$:

hole width $x$: $1\le x\le n-2$,
hole height $y$: $1\le y\le n-2$,
left offset $a$: $1\le a\le n-x-1$,
bottom offset $b$: $1\le b\le n-y-1$.

Each placement counts separately, including rotations/reflections exactly as stated.

For each lamina:

$$E(L) = \frac{ I(O,O)-2I(O,H)+I(H,H) }{ (n^2-xy)^2 }.$$

Summing all these values yields $S(n)$.

2. Python implementation

import math

# Fourth-order primitive
def Q(x, y):
    x = abs(x)
    y = abs(y)

    r = math.hypot(x, y)

    val = 0.0

    if x > 0:
        val += x**4 * y * math.asinh(y / x) / 24.0

    if y > 0:
        val += x * y**4 * math.asinh(x / y) / 24.0

    val += (
        -x**4 * r / 60.0
        + x**2 * y**2 * r / 20.0
        - y**4 * r / 60.0
    )

    return val

# Rectangle pair integral
# rectangle format: (x1, x2, y1, y2)
def rect_pair_integral(R1, R2):

    a, b, c, d = R1
    e, f, g, h = R2

    xs = [a, b]
    us = [e, f]
    ys = [c, d]
    vs = [g, h]

    total = 0.0

    for i, X in enumerate(xs):
        for j, U in enumerate(us):
            for k, Y in enumerate(ys):
                for l, V in enumerate(vs):

                    total += (
                        (-1) ** (i + j + k + l)
                        * Q(X - U, Y - V)
                    )

    return total

def S(n):

    outer = (0, n, 0, n)

    # Cache repeated rectangle computations
    cache = {}

    def RP(A, B):
        key = (A, B)

        if key not in cache:
            cache[key] = rect_pair_integral(A, B)

        return cache[key]

    outer_outer = RP(outer, outer)

    total = 0.0

    for x in range(1, n - 1):
        for y in range(1, n - 1):

            # self-integral of hole depends only on x,y
            HH = RP((0, x, 0, y), (0, x, 0, y))

            area = n * n - x * y

            for a in range(1, n - x):
                if a + x > n - 1:
                    continue

                for b in range(1, n - y):
                    if b + y > n - 1:
                        continue

                    hole = (a, a + x, b, b + y)

                    I = (
                        outer_outer
                        - 2 * RP(outer, hole)
                        + HH
                    )

                    total += I / (area * area)

    return total

print(round(S(3), 4))
print(round(S(4), 4))
print(round(S(40), 4))

3. Code walkthrough

`Q(x,y)`

Implements the closed-form fourth-order primitive derived above.

`rect_pair_integral`

Uses the 16-corner inclusion–exclusion formula:

$$\sum (-1)^{i+j+k+l}Q(\cdots).$$

This exactly evaluates

$$I(R_1,R_2).$$

`S(n)`

Enumerates every valid inner rectangle placement.

For each lamina:

compute $I(O,O)$,
subtract $2I(O,H)$,
add $I(H,H)$,
divide by area squared.

Then sum over all laminae.

4. Verification

The program reproduces the values given in the problem:

$$S(3)=1.6514$$

and

$$S(4)=19.6564.$$

So the implementation matches the statement exactly.

Finally:

$$S(40)=11730879.002273293\ldots$$

Rounded to four decimal places:

$$11730879.0023$$

Answer: 11730879.0023

Project Euler Problem 547