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.
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.
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}.$$
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\).
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}.$$
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.
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.
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\).
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.
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.
Fixiert sei das äußere Quadrat \(O=[0,n]\times[0,n]\). Eine zulässige hohle Lamina entsteht, indem man ein inneres achsenparalleles Rechteck
$$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$$
entfernt. Die verbleibende Region ist \(R=O\setminus H\) mit Fläche \(A=n^2-xy\). Wählt man zwei Punkte unabhängig und gleichverteilt in \(R\), dann ist ihr Erwartungsabstand
$$E(R)=\frac{1}{A^2}\int_R\int_R \|p-q\|\,dp\,dq.$$
Gesucht ist \(S(n)\), also die Summe von \(E(R)\) über alle zulässigen Werte von \(x,y,a,b\). Das Programm berechnet genau diese Summe und rundet das Ergebnis auf vier Nachkommastellen.
Die Lösung ersetzt das kontinuierliche Geometrieproblem durch vorab berechnete Wechselwirkungstabellen zwischen Einheitszellen. Entscheidend ist ein Distanzkern für zwei Zellen mit festem Gitterversatz; aus diesem Kern werden dann Rechtecke und schließlich die hohlen Regionen zusammengesetzt.
Zerlege jedes Rechteck in Einheitszellen. Für nichtnegative ganzzahlige Verschiebungen \((d_x,d_y)\) sei \(K(d_x,d_y)\) das doppelte Distanzintegral für zwei zufällige Punkte in zwei Zellen mit genau diesem Versatz:
$$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.$$
Mit den Differenzen \(s=u_1-u_2\) und \(t=v_1-v_2\) erhält man eine Dreiecksverteilung, also
$$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.$$
Die Implementierungen nutzen Symmetrie und werten diese Formel auf \([0,1]^2\) per Gauß-Legendre-Quadratur aus:
$$\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}$$
Der Spezialfall \(K(0,0)\) ist die bekannte mittlere Distanz zweier Zufallspunkte im Einheitsquadrat:
$$K(0,0)=\frac{2+\sqrt{2}+5\ln(1+\sqrt{2})}{15}.$$
Bezeichne mit \(F(w,h)\) das vollständige Doppelintegral der Distanz über ein volles \(w\times h\)-Rechteck:
$$F(w,h)=\int_{[0,w]\times[0,h]}\int_{[0,w]\times[0,h]} \|p-q\|\,dp\,dq.$$
Zu einem Versatz \((d_x,d_y)\) gibt es genau \((w-|d_x|)(h-|d_y|)\) geordnete Zellpaare. Daher folgt
$$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|).$$
Das ist die Grundlage der Rechtecktabelle. Gezählt wird hier das gesamte Integral, noch nicht der Erwartungswert. Für den Mittelwert innerhalb eines massiven Rechtecks müsste man anschließend durch \((wh)^2\) teilen.
Für zwei Regionen \(A\) und \(B\) schreiben wir
$$F_{AB}=\int_A\int_B \|p-q\|\,dp\,dq.$$
Da \(R=O\setminus H\) gilt und \(1_R=1_O-1_H\), erhält man
$$F_{RR}=F_{OO}-F_{OH}-F_{HO}+F_{HH}.$$
Wegen der Symmetrie der Distanz ist \(F_{OH}=F_{HO}\), also
$$F_{RR}=F_{OO}-2F_{OH}+F_{HH}.$$
Hier stammt \(F_{OO}=F(n,n)\) aus der Tabelle für das äußere Quadrat, \(F_{HH}=F(x,y)\) aus der Tabelle für die Lochgröße, und nur \(F_{OH}\) hängt von der Position \((a,b)\) ab. Ist \(F_{RR}\) bekannt, dann lautet der Erwartungsabstand
$$E(R)=\frac{F_{RR}}{(n^2-xy)^2}.$$
Für jede Einheitszelle \((i,j)\) im äußeren \(n\times n\)-Quadrat definieren wir ihre Wechselwirkung mit dem gesamten Außenquadrat durch
$$W_n(i,j)=\sum_{u=0}^{n-1}\sum_{v=0}^{n-1} K(|u-i|,|v-j|).$$
Besetzt das Loch den Zellenblock
$$a\le i\le a+x-1,\qquad b\le j\le b+y-1,$$
dann ist der Kreuzterm einfach die Summe über diesen Block:
$$F_{OH}=\sum_{i=a}^{a+x-1}\sum_{j=b}^{b+y-1} W_n(i,j).$$
Eine zweidimensionale Präfixsumme über \(W_n\) reduziert jede solche Rechteckabfrage auf \(O(1)\) Zeit. Das ist entscheidend, weil für dasselbe Außenquadrat sehr viele Lochgrößen und Positionen geprüft werden.
Zusammengesetzt ergibt sich
$$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}.$$
Genau diese Formel wird ausgewertet. Der einzige numerische Approximationsschritt steckt im Kern \(K(d_x,d_y)\); der Rest besteht aus Tabellennachschlagen, Präfixsummen und Summation.
Für \(n=3\) gibt es nur ein einziges zulässiges Loch: \(x=y=1\) an der Position \(a=b=1\). Die Lamina hat Fläche \(3^2-1=8\), also
$$S(3)=\frac{F(3,3)-2F_{OH}+F(1,1)}{8^2}\approx 1.6514.$$
Das stimmt mit dem Kontrollwert der Implementierungen überein.
Für \(n=4\) ist die Anzahl zulässiger Laminae
$$\sum_{x=1}^{2}\sum_{y=1}^{2}(4-x-1)(4-y-1)=9.$$
Diese \(9\) setzen sich zusammen aus \(4\) Positionen für ein \(1\times 1\)-Loch, \(2\) für \(1\times 2\), \(2\) für \(2\times 1\) und \(1\) für \(2\times 2\).
Die C++, Python- und Java-Implementierungen folgen derselben Pipeline. Zuerst erzeugen sie Gauß-Legendre-Knoten und -Gewichte auf \([0,1]\) aus Nullstellen der Legendre-Polynome und deren Ableitungen. Damit füllen sie die symmetrische Kerntabelle \(K(d_x,d_y)\) für alle Versätze \(0\le d_x,d_y\le n-1\).
Danach wird für jedes \(1\le w,h\le n\) die Tabelle der Rechteck-Selbstwechselwirkungen \(F(w,h)\) aufgebaut. Für das konkrete \(n\) berechnen die Programme anschließend das Außen-zu-Zelle-Gitter \(W_n(i,j)\) und daraus eine zweidimensionale Präfixsumme, sodass jeder positionsabhängige Kreuzterm \(F_{OH}\) mit einer einzigen Rechteckabfrage gewonnen wird.
Zum Schluss iterieren sie über alle zulässigen Lochgrößen und Positionen, setzen \(F_{OO}\), \(F_{OH}\) und \(F_{HH}\) per Inklusion-Exklusion zusammen, teilen durch die quadrierte Laminafläche und addieren alles zur Gesamtsumme. Die C++-Implementierung verteilt dabei unabhängige Lochgrößenpaare zusätzlich auf mehrere Threads, ohne die Mathematik zu verändern.
Sei \(q\) die Ordnung der Gauß-Legendre-Quadratur. Der Aufbau der Kerntabelle benötigt \(O(n^2q^2)\) Zeit und \(O(n^2)\) Speicher. Die Rechtecktabelle kostet \(O(n^4)\) Zeit, und auch das Außen-Interaktionsgitter mit Präfixsummen kostet \(O(n^4)\) Zeit. Der abschließende Durchlauf über alle Lochgrößen und Positionen hat
$$\sum_{x=1}^{n-2}\sum_{y=1}^{n-2}(n-x-1)(n-y-1)=O(n^4)$$
Terme und nach der Vorverarbeitung nur \(O(1)\) Arbeit pro Term. Für feste Quadraturordnung liegt die Gesamtlaufzeit also bei \(O(n^4)\), der Speicherverbrauch bei \(O(n^2)\).
Dış kareyi \(O=[0,n]\times[0,n]\) olarak sabitleyelim. Geçerli bir boş lamina, içeriden şu eksenlere paralel dikdörtgen çıkarılarak elde edilir:
$$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.$$
Geriye kalan bölge \(R=O\setminus H\) olur ve alanı \(A=n^2-xy\)'dir. \(R\) içinden bağımsız ve düzgün dağılımla seçilen iki noktanın beklenen uzaklığı
$$E(R)=\frac{1}{A^2}\int_R\int_R \|p-q\|\,dp\,dq$$
şeklindedir. \(S(n)\), bu değerin tüm geçerli \(x,y,a,b\) seçimleri üzerindeki toplamıdır. Program tam olarak bu toplamı hesaplayıp sonucu dört ondalık basamağa yuvarlar.
Çözüm, sürekli geometri problemini sonlu sayıda tabloya indirger. Temel nesne, iki birim hücre arasındaki uzaklık çekirdeğidir; daha sonra bütün dikdörtgenler ve laminalar bu çekirdekten türetilir.
Her dikdörtgeni birim hücrelere ayıralım. Negatif olmayan tamsayı ofsetler \((d_x,d_y)\) için, alt-sol köşeleri bu kadar kaymış iki hücredeki rasgele iki noktanın toplam uzaklık integralini \(K(d_x,d_y)\) ile gösterelim:
$$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.$$
\(s=u_1-u_2\) ve \(t=v_1-v_2\) farkları alındığında yoğunluk üçgensel hale gelir ve
$$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$$
yazılır. Uygulamalar simetriyi kullanıp bu ifadeyi \([0,1]^2\) üzerinde Gauss-Legendre ile sayısal olarak hesaplar:
$$\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}$$
Özel durum olan \(K(0,0)\), birim karede iki rastgele nokta arasındaki klasik ortalama uzaklıktır:
$$K(0,0)=\frac{2+\sqrt{2}+5\ln(1+\sqrt{2})}{15}.$$
\(F(w,h)\), dolu bir \(w\times h\) dikdörtgen üzerindeki toplam çift integrali göstersin:
$$F(w,h)=\int_{[0,w]\times[0,h]}\int_{[0,w]\times[0,h]} \|p-q\|\,dp\,dq.$$
Ofseti \((d_x,d_y)\) olan tam olarak \((w-|d_x|)(h-|d_y|)\) adet sıralı hücre çifti vardır. Bu yüzden
$$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|)$$
olur. Bu formül önceden hesaplanan dikdörtgen tablosunun temelidir. Burada hesaplanan şey beklenen değer değil, toplam integraldir; yalnızca dolu dikdörtgenin ortalama uzaklığı istenseydi \((wh)^2\)'ye bölmek gerekirdi.
Herhangi iki bölge \(A\) ve \(B\) için
$$F_{AB}=\int_A\int_B \|p-q\|\,dp\,dq$$
tanımını yapalım. \(R=O\setminus H\) ve \(1_R=1_O-1_H\) olduğundan
$$F_{RR}=F_{OO}-F_{OH}-F_{HO}+F_{HH}$$
elde edilir. Uzaklık simetrik olduğu için \(F_{OH}=F_{HO}\), dolayısıyla
$$F_{RR}=F_{OO}-2F_{OH}+F_{HH}.$$
Burada \(F_{OO}=F(n,n)\) dış kare tablosundan, \(F_{HH}=F(x,y)\) delik boyutu tablosundan gelir; konuma \((a,b)\) bağlı olan tek terim \(F_{OH}\)'dir. Sonra o laminanın beklenen uzaklığı
$$E(R)=\frac{F_{RR}}{(n^2-xy)^2}$$
olur.
\(n\times n\) dış karedeki her \((i,j)\) birim hücresi için, bütün dış kare ile etkileşimini
$$W_n(i,j)=\sum_{u=0}^{n-1}\sum_{v=0}^{n-1} K(|u-i|,|v-j|)$$
şeklinde tanımlayalım. Delik şu hücre bloğunu kaplıyorsa
$$a\le i\le a+x-1,\qquad b\le j\le b+y-1,$$
çapraz terim doğrudan bu blok toplamıdır:
$$F_{OH}=\sum_{i=a}^{a+x-1}\sum_{j=b}^{b+y-1} W_n(i,j).$$
\(W_n\) tablosu üzerinde iki boyutlu prefix sum kurulduğunda, her dikdörtgen sorgusu \(O(1)\) sürede alınır. Bu hızlanma kritik önemdedir, çünkü aynı dış kare için çok sayıda delik boyutu ve konumu denenir.
Önceki adımlar birleşince
$$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}$$
elde edilir. Uygulamaların değerlendirdiği ifade tam olarak budur. Sayısal yaklaşım yalnızca çekirdek tablosu \(K(d_x,d_y)\) içinde vardır; geri kalan her şey tablo erişimi, prefix sum ve toplamadan ibarettir.
\(n=3\) için yalnızca tek bir delik mümkündür: \(x=y=1\) ve \(a=b=1\). Lamina alanı \(3^2-1=8\) olduğundan
$$S(3)=\frac{F(3,3)-2F_{OH}+F(1,1)}{8^2}\approx 1.6514.$$
Bu değer uygulamaların kullandığı kontrol noktasına tam uyar.
\(n=4\) için geçerli lamina sayısı
$$\sum_{x=1}^{2}\sum_{y=1}^{2}(4-x-1)(4-y-1)=9$$
olur. Bunlar \(1\times 1\) delik için \(4\), \(1\times 2\) için \(2\), \(2\times 1\) için \(2\) ve \(2\times 2\) için \(1\) yerleşimden gelir.
C++, Python ve Java implementasyonlarının izlediği hat aynıdır. Önce Legendre polinomlarının kökleri ve türevleri kullanılarak \([0,1]\) aralığında Gauss-Legendre düğümleri ve ağırlıkları üretilir. Bunlarla \(0\le d_x,d_y\le n-1\) için simetrik çekirdek tablosu \(K(d_x,d_y)\) doldurulur.
Daha sonra her \(1\le w,h\le n\) için dikdörtgen öz-etkileşim tablosu \(F(w,h)\) kuruludur. İstenen \(n\) için ayrıca dış-kareden tek hücreye etkileşim tablosu \(W_n(i,j)\) hesaplanır ve iki boyutlu prefix sum haline getirilir; böylece yerleşime bağlı her \(F_{OH}\) değeri tek bir dikdörtgen sorgusuyla alınır.
Son adımda tüm geçerli delik boyutları ve konumları gezilir; \(F_{OO}\), \(F_{OH}\) ve \(F_{HH}\) dahil et-çıkar ile birleştirilir, kare alanının karesine bölünür ve toplama eklenir. C++ implementasyonu ayrıca bağımsız delik-boyutu çiftlerini iş parçacıklarına dağıtarak duvar saati süresini azaltır; Python ve Java aynı matematiği seri biçimde uygular.
\(q\), Gauss-Legendre mertebesi olsun. Çekirdek tablosunun kurulması \(O(n^2q^2)\) zaman ve \(O(n^2)\) bellek gerektirir. Dikdörtgen öz-etkileşim tablosu \(O(n^4)\) zamanda, dış-etkileşim tablosu ve prefix sum da yine \(O(n^4)\) zamanda kurulur. Son tarama ise
$$\sum_{x=1}^{n-2}\sum_{y=1}^{n-2}(n-x-1)(n-y-1)=O(n^4)$$
adet yerleşim içerir ve önişlemeden sonra her biri \(O(1)\) iş ister. Bu yüzden sabit quadrature mertebesinde toplam süre \(O(n^4)\), bellek kullanımı \(O(n^2)\)'dir.
Fijemos el cuadrado exterior \(O=[0,n]\times[0,n]\). Una lámina hueca válida se obtiene quitando un rectángulo interior alineado con los ejes:
$$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.$$
La región restante es \(R=O\setminus H\), con área \(A=n^2-xy\). Si se eligen dos puntos de manera independiente y uniforme dentro de \(R\), su distancia esperada es
$$E(R)=\frac{1}{A^2}\int_R\int_R \|p-q\|\,dp\,dq.$$
La cantidad \(S(n)\) es la suma de \(E(R)\) sobre todas las elecciones válidas de \(x,y,a,b\). El programa calcula exactamente esa suma y redondea el resultado final a cuatro decimales.
La solución transforma la geometría continua en un conjunto finito de tablas precomputadas. El objeto central es un núcleo de distancia entre dos celdas unitarias; a partir de él se reconstruyen rectángulos completos y, finalmente, las láminas huecas.
Partimos cada rectángulo en celdas unitarias. Para desplazamientos enteros no negativos \((d_x,d_y)\), definimos \(K(d_x,d_y)\) como la integral total de distancia para dos puntos aleatorios tomados de dos celdas separadas por ese desplazamiento:
$$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.$$
Al introducir las diferencias \(s=u_1-u_2\) y \(t=v_1-v_2\), aparece una densidad triangular, y por tanto
$$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.$$
Las implementaciones aprovechan la simetría y evalúan esta fórmula en \([0,1]^2\) mediante cuadratura de Gauss-Legendre:
$$\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}$$
El caso especial \(K(0,0)\) es la distancia media clásica entre dos puntos aleatorios del cuadrado unidad:
$$K(0,0)=\frac{2+\sqrt{2}+5\ln(1+\sqrt{2})}{15}.$$
Sea \(F(w,h)\) la integral doble total de la distancia sobre un rectángulo sólido de tamaño \(w\times h\):
$$F(w,h)=\int_{[0,w]\times[0,h]}\int_{[0,w]\times[0,h]} \|p-q\|\,dp\,dq.$$
Para un desplazamiento \((d_x,d_y)\) existen exactamente \((w-|d_x|)(h-|d_y|)\) pares ordenados de celdas con ese desplazamiento. Luego
$$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|).$$
Esta es la base de la tabla de rectángulos. La fórmula entrega la integral completa, no la esperanza. Si quisiéramos solo la distancia media en un rectángulo macizo, dividiríamos después por \((wh)^2\).
Para dos regiones \(A\) y \(B\), escribimos
$$F_{AB}=\int_A\int_B \|p-q\|\,dp\,dq.$$
Como \(R=O\setminus H\) y \(1_R=1_O-1_H\), se obtiene
$$F_{RR}=F_{OO}-F_{OH}-F_{HO}+F_{HH}.$$
La distancia es simétrica, así que \(F_{OH}=F_{HO}\), y por tanto
$$F_{RR}=F_{OO}-2F_{OH}+F_{HH}.$$
Aquí \(F_{OO}=F(n,n)\) viene de la tabla del cuadrado exterior, \(F_{HH}=F(x,y)\) de la tabla del hueco, y el único término que depende de la posición \((a,b)\) es \(F_{OH}\). Una vez conocido \(F_{RR}\), la distancia esperada de esa lámina es
$$E(R)=\frac{F_{RR}}{(n^2-xy)^2}.$$
Para cada celda \((i,j)\) del cuadrado exterior \(n\times n\), definimos su interacción con todo el cuadrado exterior como
$$W_n(i,j)=\sum_{u=0}^{n-1}\sum_{v=0}^{n-1} K(|u-i|,|v-j|).$$
Si el hueco ocupa el bloque
$$a\le i\le a+x-1,\qquad b\le j\le b+y-1,$$
entonces el término cruzado es simplemente la suma de \(W_n(i,j)\) sobre ese bloque:
$$F_{OH}=\sum_{i=a}^{a+x-1}\sum_{j=b}^{b+y-1} W_n(i,j).$$
Una suma prefija bidimensional convierte cada consulta rectangular en \(O(1)\), lo cual es esencial porque el mismo cuadrado exterior debe combinarse con muchas dimensiones y posiciones del hueco.
Uniendo todo lo anterior resulta
$$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}.$$
Esa es exactamente la expresión que evalúan las implementaciones. La única aproximación numérica está en el núcleo \(K(d_x,d_y)\); el resto del método es acceso a tablas, extracción por prefijos y acumulación.
Cuando \(n=3\), solo existe un hueco válido: \(x=y=1\) y \(a=b=1\). El área de la lámina es \(3^2-1=8\), de modo que
$$S(3)=\frac{F(3,3)-2F_{OH}+F(1,1)}{8^2}\approx 1.6514.$$
Ese valor coincide con el punto de control usado por las implementaciones.
Para \(n=4\), el número de láminas admisibles es
$$\sum_{x=1}^{2}\sum_{y=1}^{2}(4-x-1)(4-y-1)=9.$$
Se descompone en \(4\) posiciones para un hueco \(1\times 1\), \(2\) para \(1\times 2\), \(2\) para \(2\times 1\) y \(1\) para \(2\times 2\).
Las implementaciones en C++, Python y Java siguen la misma secuencia. Primero generan nodos y pesos de Gauss-Legendre sobre \([0,1]\) a partir de raíces de polinomios de Legendre y sus derivadas. Con esos nodos llenan la tabla simétrica del núcleo \(K(d_x,d_y)\) para todos los desplazamientos \(0\le d_x,d_y\le n-1\).
Después construyen la tabla de auto-interacción rectangular \(F(w,h)\) para cada \(1\le w,h\le n\). Para el valor concreto de \(n\), calculan además la rejilla \(W_n(i,j)\) de interacción entre todo el exterior y cada celda individual, y a partir de ella crean una suma prefija 2D que permite recuperar cualquier \(F_{OH}\) con una sola consulta rectangular.
Por último recorren todos los tamaños de hueco y todas las posiciones legales, combinan \(F_{OO}\), \(F_{OH}\) y \(F_{HH}\) mediante inclusión-exclusión, dividen por el cuadrado del área de la lámina y acumulan el resultado. La versión en C++ reparte además pares independientes de tamaños del hueco entre varios hilos, sin alterar la fórmula matemática.
Sea \(q\) el orden de la cuadratura de Gauss-Legendre. Construir la tabla del núcleo cuesta \(O(n^2q^2)\) tiempo y \(O(n^2)\) memoria. La tabla de auto-interacción rectangular cuesta \(O(n^4)\) tiempo, y la construcción de la rejilla exterior con su suma prefija también cuesta \(O(n^4)\) tiempo. El barrido final sobre tamaños y posiciones del hueco contiene
$$\sum_{x=1}^{n-2}\sum_{y=1}^{n-2}(n-x-1)(n-y-1)=O(n^4)$$
configuraciones, con \(O(1)\) trabajo por configuración después del preprocesamiento. Por tanto, para cuadratura de orden fijo, el algoritmo total es \(O(n^4)\) en tiempo y \(O(n^2)\) en memoria.
固定外部正方形 \(O=[0,n]\times[0,n]\)。一个合法的空心层板由删除内部一个与坐标轴平行的矩形得到：
$$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.$$
剩余区域记为 \(R=O\setminus H\)，其面积为 \(A=n^2-xy\)。如果从 \(R\) 中独立且均匀地随机选取两点，那么它们的期望距离为
$$E(R)=\frac{1}{A^2}\int_R\int_R \|p-q\|\,dp\,dq.$$
\(S(n)\) 定义为对所有合法的 \(x,y,a,b\) 把这个 \(E(R)\) 加总。程序要做的正是这个总和，并把最终结果四舍五入到小数点后四位。
解法的核心是把连续区域上的距离积分转化为“单位小格之间的相互作用”。一旦知道两个单位格在给定格点偏移下的平均距离贡献，整个矩形、外方形与内洞之间的交互、以及最后的空心区域，都可以通过组合这些局部贡献得到。
把每个矩形分成单位小格。对于非负整数偏移 \((d_x,d_y)\)，定义 \(K(d_x,d_y)\) 为两点分别落在两个单位格中时的总距离积分，其中这两个单位格的左下角相差 \((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.$$
令 \(s=u_1-u_2\)、\(t=v_1-v_2\)。两个独立均匀变量之差在 \([-1,1]\) 上服从三角形密度 \(1-|s|\)、\(1-|t|\)，因此同一个核还可以写成
$$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.$$
实现中利用对称性，把积分折叠到 \([0,1]^2\) 上，再用 Gauss-Legendre 求积来近似：
$$\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}$$
其中最特殊的 \(K(0,0)\) 其实有已知闭式，它正是单位正方形中两随机点之间的平均距离：
$$K(0,0)=\frac{2+\sqrt{2}+5\ln(1+\sqrt{2})}{15}.$$
记 \(F(w,h)\) 为一个实心 \(w\times h\) 矩形上的总双重距离积分：
$$F(w,h)=\int_{[0,w]\times[0,h]}\int_{[0,w]\times[0,h]} \|p-q\|\,dp\,dq.$$
如果两个单位格的相对偏移为 \((d_x,d_y)\)，那么具有这个偏移的有序格对恰好有 \((w-|d_x|)(h-|d_y|)\) 个。因此
$$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|).$$
这就是程序预处理“矩形自作用表”的数学来源。注意这里得到的是总积分，而不是已经归一化后的期望值；如果单独要求一个实心矩形中的平均距离，那么还要再除以 \((wh)^2\)。
对任意两个区域 \(A\) 与 \(B\)，定义
$$F_{AB}=\int_A\int_B \|p-q\|\,dp\,dq.$$
由于空心区域 \(R=O\setminus H\)，并且指标函数满足 \(1_R=1_O-1_H\)，于是有
$$F_{RR}=F_{OO}-F_{OH}-F_{HO}+F_{HH}.$$
距离函数是对称的，所以 \(F_{OH}=F_{HO}\)。因此真正使用的公式是
$$F_{RR}=F_{OO}-2F_{OH}+F_{HH}.$$
这里 \(F_{OO}=F(n,n)\) 来自外方形表，\(F_{HH}=F(x,y)\) 来自内洞尺寸表，而只有 \(F_{OH}\) 依赖于内洞的放置位置 \((a,b)\)。一旦 \(F_{RR}\) 求出，该空心层板的期望距离就是
$$E(R)=\frac{F_{RR}}{(n^2-xy)^2}.$$
对外部 \(n\times n\) 方格中的每个单位格 \((i,j)\)，定义它与整个外方形的总相互作用
$$W_n(i,j)=\sum_{u=0}^{n-1}\sum_{v=0}^{n-1} K(|u-i|,|v-j|).$$
如果内洞覆盖的单位格块为
$$a\le i\le a+x-1,\qquad b\le j\le b+y-1,$$
那么外方形与内洞之间的交互总量就是这个块上的和：
$$F_{OH}=\sum_{i=a}^{a+x-1}\sum_{j=b}^{b+y-1} W_n(i,j).$$
实现中先对 \(W_n\) 建立二维前缀和，这样每个矩形块的求和都能在 \(O(1)\) 时间内完成。因为同一个外方形要与大量不同的洞尺寸和洞位置配对，这一步的加速非常关键。
把前面各部分合并起来，可得
$$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}.$$
这就是程序直接计算的表达式。整个方法中唯一的数值近似发生在核表 \(K(d_x,d_y)\) 的求积阶段；其余部分都是确定性的表格累积与查询。
当 \(n=3\) 时，只存在一个合法内洞：\(x=y=1\)，并且必须放在 \(a=b=1\) 的中心位置。此时空心区域面积为 \(3^2-1=8\)，因此
$$S(3)=\frac{F(3,3)-2F_{OH}+F(1,1)}{8^2}\approx 1.6514.$$
这正好对应实现里使用的校验值。
当 \(n=4\) 时，合法空心层板的总数为
$$\sum_{x=1}^{2}\sum_{y=1}^{2}(4-x-1)(4-y-1)=9.$$
具体来说，\(1\times 1\) 洞有 \(4\) 种放法，\(1\times 2\) 有 \(2\) 种，\(2\times 1\) 有 \(2\) 种，\(2\times 2\) 只有 \(1\) 种。
C++、Python 和 Java 实现遵循同一条计算链。第一步是根据 Legendre 多项式的根及其导数，在 \([0,1]\) 上生成 Gauss-Legendre 节点与权重。随后用这些节点填充所有 \(0\le d_x,d_y\le n-1\) 的对称核表 \(K(d_x,d_y)\)。
接下来，程序为所有 \(1\le w,h\le n\) 建立矩形总距离积分表 \(F(w,h)\)。对于给定的 \(n\)，程序再计算外方形到单个单位格的相互作用表 \(W_n(i,j)\)，并把它转成二维前缀和，从而每个与位置相关的 \(F_{OH}\) 都能通过一次矩形查询取得。
最后，程序遍历所有合法的洞尺寸与放置位置，用容斥组合 \(F_{OO}\)、\(F_{OH}\) 和 \(F_{HH}\)，再除以空心区域面积的平方并累加到总答案里。C++ 实现还会把彼此独立的洞尺寸对分配到多个工作线程上，但计算公式与 Python、Java 完全一致。
设 Gauss-Legendre 的阶数为 \(q\)。核表的构造需要 \(O(n^2q^2)\) 时间和 \(O(n^2)\) 内存。矩形总积分表需要 \(O(n^4)\) 时间；外方形相互作用表及其前缀和也需要 \(O(n^4)\) 时间。最终对所有洞尺寸和位置的遍历共有
$$\sum_{x=1}^{n-2}\sum_{y=1}^{n-2}(n-x-1)(n-y-1)=O(n^4)$$
项，而每一项在预处理后只需 \(O(1)\) 工作。因此在固定求积阶数下，整个算法的总时间复杂度是 \(O(n^4)\)，空间复杂度是 \(O(n^2)\)。
Зафиксируем внешний квадрат \(O=[0,n]\times[0,n]\). Допустимая полая ламина получается удалением внутреннего прямоугольника, параллельного осям:
$$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.$$
Оставшаяся область равна \(R=O\setminus H\), ее площадь \(A=n^2-xy\). Если выбрать две точки независимо и равномерно внутри \(R\), то их ожидаемое расстояние равно
$$E(R)=\frac{1}{A^2}\int_R\int_R \|p-q\|\,dp\,dq.$$
Величина \(S(n)\) определяется как сумма \(E(R)\) по всем допустимым наборам \(x,y,a,b\). Программа вычисляет именно эту сумму и выводит результат, округленный до четырех знаков после запятой.
Решение сводит непрерывную геометрию к конечному набору таблиц. Главный объект здесь — ядро расстояния между двумя единичными клетками при фиксированном сдвиге. После его предвычисления можно собирать интегралы для сплошных прямоугольников, а затем и для полой области.
Разобьем каждый прямоугольник на единичные клетки. Для неотрицательных целых сдвигов \((d_x,d_y)\) обозначим через \(K(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.$$
Если перейти к разностям \(s=u_1-u_2\) и \(t=v_1-v_2\), возникает треугольная плотность, и формула становится
$$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.$$
Реализация использует симметрию, сворачивает интеграл к области \([0,1]^2\) и считает его численно квадратурой Гаусса-Лежандра:
$$\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}$$
Частный случай \(K(0,0)\) известен в замкнутом виде: это среднее расстояние между двумя случайными точками в единичном квадрате,
$$K(0,0)=\frac{2+\sqrt{2}+5\ln(1+\sqrt{2})}{15}.$$
Пусть \(F(w,h)\) — полный двойной интеграл расстояния по сплошному прямоугольнику размера \(w\times h\):
$$F(w,h)=\int_{[0,w]\times[0,h]}\int_{[0,w]\times[0,h]} \|p-q\|\,dp\,dq.$$
Для фиксированного сдвига \((d_x,d_y)\) существует ровно \((w-|d_x|)(h-|d_y|)\) упорядоченных пар клеток с таким сдвигом. Поэтому
$$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|).$$
Именно эта формула лежит в основе таблицы прямоугольников. Она дает полный интеграл, а не среднее значение; чтобы получить среднее расстояние в сплошном прямоугольнике, нужно дополнительно делить на \((wh)^2\).
Для любых двух областей \(A\) и \(B\) введем обозначение
$$F_{AB}=\int_A\int_B \|p-q\|\,dp\,dq.$$
Так как \(R=O\setminus H\) и \(1_R=1_O-1_H\), получаем
$$F_{RR}=F_{OO}-F_{OH}-F_{HO}+F_{HH}.$$
По симметрии расстояния \(F_{OH}=F_{HO}\), следовательно
$$F_{RR}=F_{OO}-2F_{OH}+F_{HH}.$$
Здесь \(F_{OO}=F(n,n)\) берется из таблицы для внешнего квадрата, \(F_{HH}=F(x,y)\) — из таблицы для размеров отверстия, а только \(F_{OH}\) зависит от положения \((a,b)\). После этого ожидаемое расстояние для данной ламины равно
$$E(R)=\frac{F_{RR}}{(n^2-xy)^2}.$$
Для каждой единичной клетки \((i,j)\) внутри внешнего квадрата \(n\times n\) определим ее взаимодействие со всем внешним квадратом:
$$W_n(i,j)=\sum_{u=0}^{n-1}\sum_{v=0}^{n-1} K(|u-i|,|v-j|).$$
Если отверстие занимает блок клеток
$$a\le i\le a+x-1,\qquad b\le j\le b+y-1,$$
то перекрестный член равен сумме \(W_n(i,j)\) по этому блоку:
$$F_{OH}=\sum_{i=a}^{a+x-1}\sum_{j=b}^{b+y-1} W_n(i,j).$$
Двумерная префиксная сумма по таблице \(W_n\) позволяет извлекать такой прямоугольный блок за \(O(1)\). Это и дает нужное ускорение, потому что один и тот же внешний квадрат комбинируется со многими размерами и положениями отверстия.
В итоге получаем
$$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}.$$
Именно это выражение и вычисляется. Единственный численный этап — приближенное вычисление таблицы \(K(d_x,d_y)\); все остальное сводится к табличным суммам и извлечению прямоугольных блоков.
При \(n=3\) существует только одно допустимое отверстие: \(x=y=1\) и \(a=b=1\). Площадь ламины равна \(3^2-1=8\), поэтому
$$S(3)=\frac{F(3,3)-2F_{OH}+F(1,1)}{8^2}\approx 1.6514.$$
Это совпадает с контрольным значением, используемым в реализациях.
При \(n=4\) число допустимых ламин равно
$$\sum_{x=1}^{2}\sum_{y=1}^{2}(4-x-1)(4-y-1)=9.$$
Оно складывается из \(4\) положений для отверстия \(1\times 1\), \(2\) для \(1\times 2\), \(2\) для \(2\times 1\) и \(1\) для \(2\times 2\).
Реализации на C++, Python и Java следуют одной и той же схеме. Сначала они строят узлы и веса Гаусса-Лежандра на \([0,1]\), используя корни полиномов Лежандра и их производные. Затем по этим узлам заполняют симметричную таблицу ядра \(K(d_x,d_y)\) для всех сдвигов \(0\le d_x,d_y\le n-1\).
После этого строится таблица полных прямоугольных интегралов \(F(w,h)\) для всех \(1\le w,h\le n\). Для заданного \(n\) программа дополнительно вычисляет сетку \(W_n(i,j)\), описывающую взаимодействие всего внешнего квадрата с каждой клеткой, и преобразует ее в двумерную префиксную сумму, чтобы любой член \(F_{OH}\) получался одной прямоугольной выборкой.
На заключительном этапе перебираются все допустимые размеры отверстия и все его позиции, значения \(F_{OO}\), \(F_{OH}\) и \(F_{HH}\) объединяются по формуле включений-исключений, результат делится на квадрат площади ламины и добавляется к общей сумме. В C++-версии независимые пары размеров отверстия дополнительно распределяются по потокам, но сама математика остается той же.
Пусть \(q\) — порядок квадратуры Гаусса-Лежандра. Построение таблицы ядра требует \(O(n^2q^2)\) времени и \(O(n^2)\) памяти. Таблица прямоугольных интегралов строится за \(O(n^4)\), и столько же стоит построение сетки внешних взаимодействий вместе с префиксными суммами. Финальный перебор размеров и положений отверстия содержит
$$\sum_{x=1}^{n-2}\sum_{y=1}^{n-2}(n-x-1)(n-y-1)=O(n^4)$$
слагаемых, причем после предвычислений каждое обрабатывается за \(O(1)\). Следовательно, при фиксированном порядке квадратуры общая сложность равна \(O(n^4)\) по времени и \(O(n^2)\) по памяти.
نثبت المربع الخارجي \(O=[0,n]\times[0,n]\). وتتكوّن الصفيحة المجوفة المسموح بها بحذف مستطيل داخلي موازٍ للمحاور:
$$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.$$
المنطقة الباقية هي \(R=O\setminus H\)، ومساحتها \(A=n^2-xy\). إذا اخترنا نقطتين بصورة مستقلة وبانتظام من \(R\)، فإن المسافة المتوقعة بينهما تساوي
$$E(R)=\frac{1}{A^2}\int_R\int_R \|p-q\|\,dp\,dq.$$
والكمية \(S(n)\) هي مجموع \(E(R)\) على جميع الاختيارات المسموح بها لـ \(x,y,a,b\). البرنامج يحسب هذا المجموع مباشرة ثم يقرّب الناتج النهائي إلى أربع منازل عشرية.
الفكرة الأساسية هي تحويل التكامل الهندسي المستمر إلى جداول مسبقة الحساب بين خلايا وحدوية. فإذا عرفنا مساهمة المسافة بين خليتين عند إزاحة شبكية محددة، أمكننا تركيب تكامل أي مستطيل ممتلئ، ثم تطبيق مبدأ الاحتواء والاستبعاد على المنطقة المجوفة.
نقسم كل مستطيل إلى خلايا وحدة. ولكل إزاحة صحيحة غير سالبة \((d_x,d_y)\)، نعرّف \(K(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.$$
بعد الانتقال إلى الفرقين \(s=u_1-u_2\) و\(t=v_1-v_2\)، تظهر كثافة مثلثية، ولذلك يمكن كتابة المقدار نفسه على الصورة
$$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.$$
وتستفيد التطبيقات من التناظر لتقليص التكامل إلى \([0,1]^2\)، ثم تحسبه عددياً بطريقة Gauss-Legendre:
$$\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}$$
أما الحالة الخاصة \(K(0,0)\) فلها صيغة مغلقة معروفة، وهي متوسط المسافة بين نقطتين عشوائيتين داخل مربع الوحدة:
$$K(0,0)=\frac{2+\sqrt{2}+5\ln(1+\sqrt{2})}{15}.$$
لنرمز بـ \(F(w,h)\) إلى التكامل المزدوج الكلي للمسافة على مستطيل ممتلئ أبعاده \(w\times h\):
$$F(w,h)=\int_{[0,w]\times[0,h]}\int_{[0,w]\times[0,h]} \|p-q\|\,dp\,dq.$$
إذا كان فرق الإزاحة بين خليتين هو \((d_x,d_y)\)، فعدد الأزواج المرتبة من الخلايا التي تحقق ذلك يساوي تماماً \((w-|d_x|)(h-|d_y|)\). ومن ثم
$$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|).$$
هذه الصيغة هي أساس جدول المستطيلات المسبق. وهي تعطي التكامل الكلي، لا متوسط المسافة بعد التطبيع؛ فإذا أردنا متوسط المسافة داخل مستطيل ممتلئ فقط، لوجب القسمة على \((wh)^2\).
لكل منطقتين \(A\) و\(B\) نكتب
$$F_{AB}=\int_A\int_B \|p-q\|\,dp\,dq.$$
وبما أن \(R=O\setminus H\) وأن \(1_R=1_O-1_H\)، فإننا نحصل على
$$F_{RR}=F_{OO}-F_{OH}-F_{HO}+F_{HH}.$$
ولأن المسافة متناظرة فلدينا \(F_{OH}=F_{HO}\)، وبالتالي
$$F_{RR}=F_{OO}-2F_{OH}+F_{HH}.$$
هنا \(F_{OO}=F(n,n)\) يأتي من جدول المربع الخارجي، و\(F_{HH}=F(x,y)\) يأتي من جدول أحجام الفتحة، أما الحد الوحيد الذي يعتمد على الموضع \((a,b)\) فهو \(F_{OH}\). وبعد معرفة \(F_{RR}\) تصبح المسافة المتوقعة لهذه الصفيحة
$$E(R)=\frac{F_{RR}}{(n^2-xy)^2}.$$
لكل خلية وحدة \((i,j)\) داخل المربع الخارجي \(n\times n\)، نعرّف تفاعلها مع المربع الخارجي كله بالكمية
$$W_n(i,j)=\sum_{u=0}^{n-1}\sum_{v=0}^{n-1} K(|u-i|,|v-j|).$$
إذا شغلت الفتحة الكتلة الخلوية
$$a\le i\le a+x-1,\qquad b\le j\le b+y-1,$$
فإن الحد المتقاطع يساوي مجموع \(W_n(i,j)\) على تلك الكتلة:
$$F_{OH}=\sum_{i=a}^{a+x-1}\sum_{j=b}^{b+y-1} W_n(i,j).$$
وعند بناء مجموع سابق ثنائي الأبعاد لجدول \(W_n\)، تصبح كل عملية جمع على مستطيل فرعي منتهية في \(O(1)\) زمن. وهذا ضروري لأن المربع الخارجي نفسه يُقرن بعدد كبير من أحجام الفتحات ومواضعها.
بضم الخطوات السابقة نحصل على
$$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}.$$
وهذا هو التعبير الذي تنفذه التطبيقات حرفياً. التقريب العددي الوحيد يقع في حساب نواة المسافة \(K(d_x,d_y)\)، أما بقية العمل فهو تجميع حتمي باستخدام الجداول والمجاميع السابقة.
عندما \(n=3\)، لا توجد إلا فتحة واحدة ممكنة: \(x=y=1\) عند الموضع \(a=b=1\). مساحة الصفيحة هي \(3^2-1=8\)، ومن ثم
$$S(3)=\frac{F(3,3)-2F_{OH}+F(1,1)}{8^2}\approx 1.6514.$$
وهذا يطابق قيمة التحقق المستخدمة في التطبيقات.
أما عند \(n=4\)، فإن عدد الصفائح المجوفة المسموح بها يساوي
$$\sum_{x=1}^{2}\sum_{y=1}^{2}(4-x-1)(4-y-1)=9.$$
وهذه \(9\) تأتي من \(4\) مواضع لفتحة \(1\times 1\)، و\(2\) لفتحة \(1\times 2\)، و\(2\) لفتحة \(2\times 1\)، وموضع واحد فقط لفتحة \(2\times 2\).
تتبع تطبيقات C++ وPython وJava المسار الحسابي نفسه. أولاً تُولَّد عقد وأوزان Gauss-Legendre على الفترة \([0,1]\) انطلاقاً من جذور كثيرات حدود Legendre ومشتقاتها. وبعد ذلك تُملأ الجدول المتناظر للنواة \(K(d_x,d_y)\) لجميع الإزاحات \(0\le d_x,d_y\le n-1\).
ثم يُبنى جدول تكاملات المستطيلات الكاملة \(F(w,h)\) لكل \(1\le w,h\le n\). وللقيمة المطلوبة من \(n\)، يُحسب أيضاً جدول \(W_n(i,j)\) الذي يصف تفاعل المربع الخارجي كله مع كل خلية منفردة، ثم يُحوَّل إلى مجموع سابق ثنائي الأبعاد بحيث يمكن استخراج أي \(F_{OH}\) باستعلام مستطيلي واحد.
وفي النهاية تُفحَص جميع أحجام الفتحات المسموح بها وجميع مواضعها القانونية، وتُدمج القيم \(F_{OO}\) و\(F_{OH}\) و\(F_{HH}\) بصيغة الاحتواء والاستبعاد، ثم تُقسَم على مربع مساحة الصفيحة وتُضاف إلى المجموع الكلي. كما أن تنفيذ C++ يوزع أزواج أحجام الفتحات المستقلة على عدة خيوط لتحسين الزمن الفعلي، من غير أي تغيير في الصيغة الرياضية.
إذا كانت رتبة Gauss-Legendre هي \(q\)، فإن بناء جدول النواة يحتاج إلى \(O(n^2q^2)\) زمناً و\(O(n^2)\) ذاكرة. كما أن بناء جدول تكاملات المستطيلات يكلف \(O(n^4)\) زمناً، وبناء جدول التفاعل الخارجي مع مجموعاته السابقة يكلف أيضاً \(O(n^4)\) زمناً. أما المرور النهائي على جميع أحجام الفتحات ومواضعها فيحتوي على
$$\sum_{x=1}^{n-2}\sum_{y=1}^{n-2}(n-x-1)(n-y-1)=O(n^4)$$
حالات، مع \(O(1)\) عمل لكل حالة بعد مرحلة التمهيد. لذلك، عند تثبيت رتبة التقريب العددي، تكون الكلفة الكلية \(O(n^4)\) في الزمن و\(O(n^2)\) في الذاكرة.