Problem Summary

Problem 906 studies two independent random preference orders on the same set of \(n\) alternatives and asks for a probability \(P(n)\); the final evaluation uses \(n=20000\). The implementations show that the whole computation can be organized around one distinguished alternative \(x\). If \(x\) is in position \(a+1\) in the first order and in position \(b+1\) in the second, then among the other \(N=n-1\) alternatives there are \(a\) above \(x\) in the first order and \(b\) above \(x\) in the second.

The remaining alternatives split into three blocks: those above \(x\) only in the first order, those above \(x\) only in the second order, and those below \(x\) in both orders. The event counted by the implementations requires the two upper blocks to be disjoint, and then it adds an auxiliary contribution coming from the common lower block. Averaging that contribution over the \(n^2\) equally likely rank pairs of \(x\) yields the desired probability.

Mathematical Approach

Let \(N=n-1\). For a fixed alternative \(x\), the two permutations induce a pair of rank coordinates. The solution turns the original probability question into a clean double sum over those coordinates.

Fix one alternative and its two rank coordinates

Write \(y\succ_1 x\) when the first preference order places \(y\) above \(x\), and \(y\succ_2 x\) for the second order. Define

$$A=\{y:y\succ_1 x\},\qquad B=\{y:y\succ_2 x\}.$$

If \(|A|=a\) and \(|B|=b\), then \(x\) has rank \(a+1\) in the first order and rank \(b+1\) in the second. Because the two orders are independent and each rank of \(x\) is uniform on \(\{1,\dots,n\}\), every pair \((a,b)\in\{0,\dots,N\}^2\) occurs with probability \(1/n^2\).

Why the upper sets must be disjoint

The implemented formula contributes only when no alternative lies above \(x\) in both orders. In set language that condition is

$$A\cap B=\varnothing.$$

For fixed sizes \(a\) and \(b\), the set \(B\) is a uniformly random \(b\)-subset of the \(N\) alternatives other than \(x\). Exactly \(\binom{N-a}{b}\) such subsets avoid the \(a\) elements already placed in \(A\). Therefore

$$q(a,b)=\Pr(A\cap B=\varnothing\mid |A|=a,\ |B|=b)=\frac{\binom{N-a}{b}}{\binom{N}{b}}.$$

This immediately explains the triangular summation region in the code: if \(a+b>N\), then disjoint sets of those sizes cannot exist, so the contribution is zero. It also explains the multiplicative update used inside the inner loop:

$$q(a,0)=1,\qquad q(a,b+1)=q(a,b)\frac{N-a-b}{N-b}.$$

The unanimous lower block

Once \(A\) and \(B\) are disjoint, every remaining alternative is below \(x\) in both orders. Denote that common lower block by

$$C=\{y:x\succ_1 y\ \text{and}\ x\succ_2 y\},\qquad |C|=m=N-a-b.$$

The implementations attach to this block the auxiliary quantity

$$R_m=\sum_{c=0}^{m}\prod_{t=0}^{c-1}\frac{m-t}{N-t}=\sum_{c=0}^{m}\frac{(m)_c}{(N)_c}=\sum_{c=0}^{m}\frac{\binom{m}{c}}{\binom{N}{c}}.$$

The term for a fixed \(c\) has a direct counting interpretation: \(\binom{m}{c}/\binom{N}{c}\) is the probability that a uniformly chosen \(c\)-subset of the \(N\) alternatives other than \(x\) lies entirely inside the unanimous lower block \(C\). The total \(R_m\) adds those probabilities for every feasible subset size \(c\).

Collapsing the auxiliary sum

The direct sum for \(R_m\) is already good enough for computation, but the mathematics becomes clearer after simplifying it. Start from

$$\frac{\binom{m}{c}}{\binom{N}{c}}=\frac{1}{\binom{N}{m}}\binom{N-c}{m-c}.$$

Substituting this into the previous formula gives

$$R_m=\frac{1}{\binom{N}{m}}\sum_{c=0}^{m}\binom{N-c}{m-c}.$$

Now set \(j=m-c\). The sum becomes

$$\sum_{j=0}^{m}\binom{N-m+j}{j},$$

which is a standard hockey-stick sum, so

$$R_m=\frac{\binom{N+1}{m}}{\binom{N}{m}}=\frac{N+1}{N+1-m}.$$

Since \(m=N-a-b\), this can be rewritten in the especially convenient form

$$R_{N-a-b}=\frac{n}{a+b+1}.$$

The final probability formula

Putting the two pieces together, the probability computed by the C++, Python, and Java implementations is

$$\boxed{P(n)=\frac{1}{n^2}\sum_{a=0}^{N}\sum_{b=0}^{N-a} q(a,b)\,R_{N-a-b}=\frac{1}{n^2}\sum_{a=0}^{N}\sum_{b=0}^{N-a}\frac{\binom{N-a}{b}}{\binom{N}{b}}\cdot\frac{n}{a+b+1}.}$$

Everything in the program is a fast way to evaluate this triangular double sum without ever forming large factorials or binomial coefficients explicitly.

Worked Example: \(n=3\)

Take \(n=3\), so \(N=2\). Then

$$R_0=1,\qquad R_1=1+\frac12=\frac32,\qquad R_2=1+1+1=3.$$

The admissible pairs \((a,b)\) satisfy \(a+b\le 2\). Their contributions are

$$\begin{aligned} (0,0)&:\ q=1,\ m=2,\ qR_m=3,\\ (0,1)&:\ q=1,\ m=1,\ qR_m=\frac32,\\ (0,2)&:\ q=1,\ m=0,\ qR_m=1,\\ (1,0)&:\ q=1,\ m=1,\ qR_m=\frac32,\\ (1,1)&:\ q=\frac12,\ m=0,\ qR_m=\frac12,\\ (2,0)&:\ q=1,\ m=0,\ qR_m=1. \end{aligned}$$

The total is \(17/2\), and dividing by \(n^2=9\) gives

$$P(3)=\frac{17/2}{9}=\frac{17}{18}.$$

This small case shows all the moving parts of the full solution: the disjointness factor, the lower-block sum, and the final average over rank pairs.

How the Code Works

Precomputation of reciprocal factors and lower-block values

The implementations first build a table of reciprocals \(1,2,\dots,N\). That allows every ratio in the double sum to be updated by multiplication instead of by repeated division. They also precompute all values of \(R_m\) for \(m=0,\dots,N\). For a fixed \(m\), the term with \(c=0\) is \(1\), and the next term is obtained from the previous one by multiplying by \((m-c+1)/(N-c+1)\). Summing those terms produces the required lower-block table.

The triangular sweep over \((a,b)\)

After the \(R_m\) table is available, the main work is the triangular loop over all \(a\) and \(b\) with \(0\le a\le N\) and \(0\le b\le N-a\). For each fixed \(a\), the implementation starts from \(q(a,0)=1\) and updates the disjointness factor with

$$q(a,b+1)=q(a,b)\frac{N-a-b}{N-b}.$$

At each lattice point \((a,b)\), it looks up \(R_{N-a-b}\), adds \(q(a,b)R_{N-a-b}\) to the running total, and finally divides the accumulated sum by \(n^2\).

Numerical accuracy and compiled-language optimizations

The calculation is purely floating-point. The compiled implementations therefore use compensated summation to reduce accumulation error, and the C++ version also keeps the main total in extended precision. The compiled versions split the outer \(a\)-range into blocks and process those blocks in parallel, while the Python implementation follows the same mathematics serially.

Complexity Analysis

Let \(N=n-1\). Building the reciprocal table costs \(O(N)\) time and \(O(N)\) memory. Precomputing all \(R_m\) values with the direct inner products used in the implementations costs \(O(N^2)\) arithmetic operations. The main triangular sum over \((a,b)\) is also \(O(N^2)\).

Therefore the overall running time is \(O(N^2)\) and the memory usage is \(O(N)\). Parallel execution reduces wall-clock time for the compiled versions, but it does not change the asymptotic order.

Footnotes and References

  1. Problem page: Project Euler 906
  2. Total order: Wikipedia - Total order
  3. Binomial coefficient: Wikipedia - Binomial coefficient
  4. Falling and rising factorials: Wikipedia - Falling and rising factorials
  5. Hypergeometric distribution: Wikipedia - Hypergeometric distribution
  6. Hockey-stick identity: Wikipedia - Hockey-stick identity
  7. Kahan summation algorithm: Wikipedia - Kahan summation algorithm

Problem Summary / Problemzusammenfassung

Problem 906 fragt nach einer Wahrscheinlichkeit \(P(n)\); in der endgültigen Auswertung wird \(n=20000\) verwendet. Aus den Implementierungen erkennt man, dass man die gesamte Rechnung um eine ausgezeichnete Alternative \(x\) herum organisieren kann. Steht \(x\) in der ersten Präferenzordnung auf Position \(a+1\) und in der zweiten auf Position \(b+1\), dann gibt es unter den übrigen \(N=n-1\) Alternativen genau \(a\) Stück über \(x\) in der ersten Ordnung und \(b\) Stück über \(x\) in der zweiten.

Die restlichen Alternativen zerfallen dann in drei Blöcke: nur in der ersten Ordnung über \(x\), nur in der zweiten Ordnung über \(x\) und in beiden Ordnungen unter \(x\). Das von den Implementierungen ausgewertete Ereignis verlangt, dass die beiden oberen Blöcke disjunkt sind; danach kommt noch ein Hilfsbeitrag aus dem gemeinsamen unteren Block hinzu. Mittelt man diesen Beitrag über die \(n^2\) gleich wahrscheinlichen Rangpaare von \(x\), erhält man die gesuchte Wahrscheinlichkeit.

Mathematical Approach / Mathematischer Ansatz

Setze \(N=n-1\). Für eine feste Alternative \(x\) erzeugen die beiden Permutationen ein Paar von Rangkoordinaten. Die Lösung verwandelt die ursprüngliche Wahrscheinlichkeitsfrage in eine saubere doppelte Summe über diese Koordinaten.

Eine Alternative fixieren und ihre beiden Rangkoordinaten betrachten

Schreibe \(y\succ_1 x\), wenn die erste Präferenzordnung \(y\) über \(x\) stellt, und \(y\succ_2 x\) für die zweite Ordnung. Definiere

$$A=\{y:y\succ_1 x\},\qquad B=\{y:y\succ_2 x\}.$$

Gilt \(|A|=a\) und \(|B|=b\), so hat \(x\) Rang \(a+1\) in der ersten und Rang \(b+1\) in der zweiten Ordnung. Weil die beiden Ordnungen unabhängig sind und jeder Rang von \(x\) gleichverteilt auf \(\{1,\dots,n\}\) ist, tritt jedes Paar \((a,b)\in\{0,\dots,N\}^2\) mit Wahrscheinlichkeit \(1/n^2\) auf.

Warum die oberen Mengen disjunkt sein muessen

Die implementierte Formel liefert nur dann einen Beitrag, wenn keine Alternative \(x\) in beiden Ordnungen zugleich schlaegt. In Mengenschreibweise lautet diese Bedingung

$$A\cap B=\varnothing.$$

Bei festen Groessen \(a\) und \(b\) ist \(B\) eine gleichverteilte \(b\)-Teilmenge der \(N\) Alternativen ausser \(x\). Genau \(\binom{N-a}{b}\) dieser Teilmengen vermeiden die \(a\) Elemente, die bereits in \(A\) liegen. Also gilt

$$q(a,b)=\Pr(A\cap B=\varnothing\mid |A|=a,\ |B|=b)=\frac{\binom{N-a}{b}}{\binom{N}{b}}.$$

Damit wird auch die dreieckige Summationsregion im Code sofort klar: Falls \(a+b>N\), koennen disjunkte Mengen dieser Groessen nicht existieren, also ist der Beitrag null. Ebenso erklaert sich die multiplikative Aktualisierung in der inneren Schleife:

$$q(a,0)=1,\qquad q(a,b+1)=q(a,b)\frac{N-a-b}{N-b}.$$

Der gemeinsame untere Block

Sobald \(A\) und \(B\) disjunkt sind, liegt jede verbleibende Alternative in beiden Ordnungen unter \(x\). Nenne diesen gemeinsamen unteren Block

$$C=\{y:x\succ_1 y\ \text{and}\ x\succ_2 y\},\qquad |C|=m=N-a-b.$$

Die Implementierungen verknuepfen mit diesem Block die Hilfsgroesse

$$R_m=\sum_{c=0}^{m}\prod_{t=0}^{c-1}\frac{m-t}{N-t}=\sum_{c=0}^{m}\frac{(m)_c}{(N)_c}=\sum_{c=0}^{m}\frac{\binom{m}{c}}{\binom{N}{c}}.$$

Der Summand fuer festes \(c\) hat eine direkte Bedeutung: \(\binom{m}{c}/\binom{N}{c}\) ist die Wahrscheinlichkeit, dass eine gleichverteilte \(c\)-Teilmenge der \(N\) Alternativen ausser \(x\) vollstaendig im gemeinsamen unteren Block \(C\) liegt. \(R_m\) addiert diese Wahrscheinlichkeiten ueber alle moeglichen Teilmengengroessen \(c\).

Die Hilfssumme auf eine geschlossene Form bringen

Fuer die Berechnung reicht die direkte Summe bereits aus, aber mathematisch wird das Bild klarer, wenn man sie vereinfacht. Ausgangspunkt ist

$$\frac{\binom{m}{c}}{\binom{N}{c}}=\frac{1}{\binom{N}{m}}\binom{N-c}{m-c}.$$

Eingesetzt in die vorige Formel ergibt das

$$R_m=\frac{1}{\binom{N}{m}}\sum_{c=0}^{m}\binom{N-c}{m-c}.$$

Setzt man \(j=m-c\), so wird daraus

$$\sum_{j=0}^{m}\binom{N-m+j}{j},$$

und das ist eine klassische Hockey-Stick-Summe. Daher folgt

$$R_m=\frac{\binom{N+1}{m}}{\binom{N}{m}}=\frac{N+1}{N+1-m}.$$

Mit \(m=N-a-b\) erhaelt man die besonders praktische Form

$$R_{N-a-b}=\frac{n}{a+b+1}.$$

Die endgueltige Wahrscheinlichkeitsformel

Setzt man beide Bausteine zusammen, so berechnen die C++-, Python- und Java-Implementierungen

$$\boxed{P(n)=\frac{1}{n^2}\sum_{a=0}^{N}\sum_{b=0}^{N-a} q(a,b)\,R_{N-a-b}=\frac{1}{n^2}\sum_{a=0}^{N}\sum_{b=0}^{N-a}\frac{\binom{N-a}{b}}{\binom{N}{b}}\cdot\frac{n}{a+b+1}.}$$

Das gesamte Programm ist also nur eine schnelle Auswertung dieser dreieckigen Doppelsumme, ohne jemals grosse Fakultaeten oder Binomialkoeffizienten explizit aufzubauen.

Durchgerechnetes Beispiel: \(n=3\)

Fuer \(n=3\) gilt \(N=2\). Dann sind

$$R_0=1,\qquad R_1=1+\frac12=\frac32,\qquad R_2=1+1+1=3.$$

Die zulaessigen Paare \((a,b)\) erfuellen \(a+b\le 2\). Ihre Beitraege sind

$$\begin{aligned} (0,0)&:\ q=1,\ m=2,\ qR_m=3,\\ (0,1)&:\ q=1,\ m=1,\ qR_m=\frac32,\\ (0,2)&:\ q=1,\ m=0,\ qR_m=1,\\ (1,0)&:\ q=1,\ m=1,\ qR_m=\frac32,\\ (1,1)&:\ q=\frac12,\ m=0,\ qR_m=\frac12,\\ (2,0)&:\ q=1,\ m=0,\ qR_m=1. \end{aligned}$$

Die Summe betraegt \(17/2\). Nach Division durch \(n^2=9\) bleibt

$$P(3)=\frac{17/2}{9}=\frac{17}{18}.$$

Dieses kleine Beispiel zeigt bereits alle wesentlichen Bestandteile der allgemeinen Loesung: den Disjunktheitsfaktor, die Hilfssumme des unteren Blocks und den abschliessenden Mittelwert ueber alle Rangpaare.

How the Code Works / Wie der Code arbeitet

Vorberechnung von Kehrwerten und unteren Blockwerten

Zunaechst bauen die Implementierungen eine Tabelle der Kehrwerte \(1,2,\dots,N\) auf. Dadurch lassen sich alle Quotienten in der Doppelsumme durch Multiplikation statt durch wiederholte Division aktualisieren. Danach werden alle Werte \(R_m\) fuer \(m=0,\dots,N\) vorab berechnet. Fuer festes \(m\) ist der Term zu \(c=0\) gleich \(1\), und jeder folgende Term entsteht aus dem vorherigen durch Multiplikation mit \((m-c+1)/(N-c+1)\). Die Summe dieser Terme liefert die benoetigte Tabelle fuer den gemeinsamen unteren Block.

Der dreieckige Durchlauf ueber \((a,b)\)

Sobald die Tabelle \(R_m\) vorliegt, besteht die Hauptarbeit aus der dreieckigen Schleife ueber alle \(a\) und \(b\) mit \(0\le a\le N\) und \(0\le b\le N-a\). Fuer jedes feste \(a\) startet die Implementierung bei \(q(a,0)=1\) und aktualisiert den Disjunktheitsfaktor mit

$$q(a,b+1)=q(a,b)\frac{N-a-b}{N-b}.$$

An jedem Gitterpunkt \((a,b)\) wird \(R_{N-a-b}\) aus der Tabelle gelesen, \(q(a,b)R_{N-a-b}\) zur laufenden Summe addiert und am Ende alles durch \(n^2\) geteilt.

Numerische Genauigkeit und Optimierungen der kompilierten Versionen

Die Rechnung arbeitet vollstaendig mit Gleitkommazahlen. Deshalb verwenden die kompilierten Implementierungen kompensierte Summation, um Akkumulationsfehler zu verringern, und die C++-Version fuehrt die Hauptsumme sogar in erweiterter Praezision. Ausserdem teilen die kompilierten Versionen den aeusseren \(a\)-Bereich in Bloecke und verarbeiten diese parallel; die Python-Implementierung verwendet dieselbe Mathematik seriell.

Complexity Analysis / Komplexitaetsanalyse

Mit \(N=n-1\) kostet der Aufbau der Kehrwerttabelle \(O(N)\) Zeit und \(O(N)\) Speicher. Die Vorberechnung aller \(R_m\)-Werte mit den direkt implementierten inneren Produkten braucht \(O(N^2)\) arithmetische Operationen. Die dreieckige Hauptsumme ueber \((a,b)\) ist ebenfalls \(O(N^2)\).

Insgesamt ergibt sich damit Laufzeit \(O(N^2)\) bei Speicherverbrauch \(O(N)\). Parallelisierung verbessert die reale Laufzeit der kompilierten Versionen, aendert aber die asymptotische Ordnung nicht.

Footnotes and References / Fussnoten und Referenzen

  1. Problemseite: Project Euler 906
  2. Totale Ordnung: Wikipedia - Total order
  3. Binomialkoeffizient: Wikipedia - Binomial coefficient
  4. Fallende und steigende Fakultaeten: Wikipedia - Falling and rising factorials
  5. Hypergeometrische Verteilung: Wikipedia - Hypergeometric distribution
  6. Hockey-Stick-Identitaet: Wikipedia - Hockey-stick identity
  7. Kahan-Summationsalgorithmus: Wikipedia - Kahan summation algorithm

Problem Summary / Problem Özeti

Problem 906 bir olasilik degeri \(P(n)\) istiyor; son hesapta \(n=20000\) kullaniliyor. Uygulamalardan goruldugu gibi tum hesap tek bir ayri alternatif \(x\) etrafinda kurulabiliyor. \(x\), birinci tercih sirasinda \(a+1\). konumda ve ikinci tercih sirasinda \(b+1\). konumda ise, kalan \(N=n-1\) alternatifin \(a\) tanesi birinci sirada \(x\)'in ustunde, \(b\) tanesi de ikinci sirada \(x\)'in ustundedir.

Baska bir deyisle, diger alternatifler uc bloga ayrilir: yalnizca birinci sirada \(x\)'in ustunde olanlar, yalnizca ikinci sirada \(x\)'in ustunde olanlar ve iki sirada da \(x\)'in altinda olanlar. Uygulamalarin topladigi olayda bu iki ust blogun ayrik olmasi gerekir; bundan sonra ortak alt blogun urettigi ek bir katsayi gelir. Bu katkilar \(x\)'in \(n^2\) esit olasilikli sira ciftinin uzerinde ortalaninca aranan olasilik elde edilir.

Mathematical Approach / Matematiksel Yaklaşım

\(N=n-1\) olsun. Sabit bir \(x\) alternatifi icin iki permutasyon, \(x\)'e ait iki sira koordinati uretir. Cozum, asil olasilik sorusunu bu koordinatlar uzerinde tanimli duzgun bir cift toplama indirger.

Tek bir alternatifi sabitlemek ve iki sira koordinatini kaydetmek

Birinci tercih duzeni \(y\)'yi \(x\)'in ustune koyuyorsa \(y\succ_1 x\), ikinci tercih duzeni icin de \(y\succ_2 x\) yazalim. O zaman

$$A=\{y:y\succ_1 x\},\qquad B=\{y:y\succ_2 x\}.$$

\(|A|=a\) ve \(|B|=b\) ise \(x\), birinci sirada \(a+1\)., ikinci sirada \(b+1\). siradadir. Iki sira bagimsiz oldugu ve \(x\)'in her iki siradaki konumu \(\{1,\dots,n\}\) uzerinde uniform dagildigi icin, \((a,b)\in\{0,\dots,N\}^2\) ciftlerinin her biri \(1/n^2\) olasilikla ortaya cikar.

Ust kumelerin neden ayrik olmasi gerektigi

Uygulanan formulde katkı ancak hicbir alternatif \(x\)'in iki sirada birden ustunde degilse gelir. Kumeler diliyle bu kosul

$$A\cap B=\varnothing$$

seklindedir. \(a\) ve \(b\) sabitlenince \(B\), \(x\) disindaki \(N\) alternatif arasindan uniform secilen bir \(b\)-alt kumedir. Bu alt kumelerin tam \(\binom{N-a}{b}\) tanesi, zaten \(A\) icinde bulunan \(a\) elemandan kacar. Dolayisiyla

$$q(a,b)=\Pr(A\cap B=\varnothing\mid |A|=a,\ |B|=b)=\frac{\binom{N-a}{b}}{\binom{N}{b}}.$$

Bu ifade kodun neden ucgensel bir bolgede topladigini da aciklar: \(a+b>N\) ise bu boyutlarda ayrik kumeler olamaz ve katkı sifirdir. Icinci dongudeki carpimsal guncelleme de buradan gelir:

$$q(a,0)=1,\qquad q(a,b+1)=q(a,b)\frac{N-a-b}{N-b}.$$

Ortak alt blok

\(A\) ve \(B\) ayrik oldugunda kalan her alternatif iki sirada da \(x\)'in altindadir. Bu ortak alt blogu

$$C=\{y:x\succ_1 y\ \text{and}\ x\succ_2 y\},\qquad |C|=m=N-a-b$$

ile gosterelim. Uygulamalar bu blok icin su yardimci buyuklugu kullanir:

$$R_m=\sum_{c=0}^{m}\prod_{t=0}^{c-1}\frac{m-t}{N-t}=\sum_{c=0}^{m}\frac{(m)_c}{(N)_c}=\sum_{c=0}^{m}\frac{\binom{m}{c}}{\binom{N}{c}}.$$

Sabit bir \(c\) icin \(\binom{m}{c}/\binom{N}{c}\), \(x\) disindaki \(N\) alternatiften uniform secilen bir \(c\)-alt kumesinin tamamen ortak alt blok \(C\) icinde kalma olasiligidir. Yani \(R_m\), mumkun tum \(c\) boyutlari icin bu olasiliklari toplar.

Yardimci toplami kapali bicime indirmek

Hesaplama icin dogrudan toplam yeterlidir; fakat matematiksel yapiyi kapali bicim daha iyi gosterir. Baslangic kimligi

$$\frac{\binom{m}{c}}{\binom{N}{c}}=\frac{1}{\binom{N}{m}}\binom{N-c}{m-c}$$

olsun. Bunu yukaridaki ifadeye koyunca

$$R_m=\frac{1}{\binom{N}{m}}\sum_{c=0}^{m}\binom{N-c}{m-c}$$

elde edilir. \(j=m-c\) degiskenini kullanirsak toplam

$$\sum_{j=0}^{m}\binom{N-m+j}{j}$$

olur; bu klasik bir hockey-stick toplamidir. Sonuc

$$R_m=\frac{\binom{N+1}{m}}{\binom{N}{m}}=\frac{N+1}{N+1-m}$$

seklindedir. \(m=N-a-b\) yazinca daha kullanisli bicim

$$R_{N-a-b}=\frac{n}{a+b+1}$$

olur.

Son olasilik formulu

Iki parcayi birlestirince C++, Python ve Java uygulamalarinin hesapladigi buyukluk

$$\boxed{P(n)=\frac{1}{n^2}\sum_{a=0}^{N}\sum_{b=0}^{N-a} q(a,b)\,R_{N-a-b}=\frac{1}{n^2}\sum_{a=0}^{N}\sum_{b=0}^{N-a}\frac{\binom{N-a}{b}}{\binom{N}{b}}\cdot\frac{n}{a+b+1}.}$$

Dolayisiyla programin yaptigi sey, buyuk faktoriyelleri ya da binom katsayilarini acikca kurmadan bu ucgensel cift toplami hizli bicimde degerlendirmektir.

Calisilmis ornek: \(n=3\)

\(n=3\) alalim; o zaman \(N=2\). Bu durumda

$$R_0=1,\qquad R_1=1+\frac12=\frac32,\qquad R_2=1+1+1=3.$$

Uygun \((a,b)\) ciftleri \(a+b\le 2\) kosulunu saglar ve katkilar sunlardir:

$$\begin{aligned} (0,0)&:\ q=1,\ m=2,\ qR_m=3,\\ (0,1)&:\ q=1,\ m=1,\ qR_m=\frac32,\\ (0,2)&:\ q=1,\ m=0,\ qR_m=1,\\ (1,0)&:\ q=1,\ m=1,\ qR_m=\frac32,\\ (1,1)&:\ q=\frac12,\ m=0,\ qR_m=\frac12,\\ (2,0)&:\ q=1,\ m=0,\ qR_m=1. \end{aligned}$$

Toplam \(17/2\) eder. Sonra \(n^2=9\)'a bolunce

$$P(3)=\frac{17/2}{9}=\frac{17}{18}$$

elde edilir. Bu kucuk ornek, genel cozumu olusturan tum yapilari bir arada gosterir: ayriklik katsayisi, ortak alt blok toplami ve sira ciftleri uzerindeki son ortalama.

How the Code Works / Kod Nasil Calisir

Tersler ve alt blok degerleri icin on hesaplama

Uygulamalar once \(1,2,\dots,N\) icin bir tersler tablosu kurar. Boylece cift toplam icindeki oranlar tekrar tekrar bolme yapmadan carpma ile guncellenebilir. Ardindan \(m=0,\dots,N\) icin tum \(R_m\) degerleri hesaplanir. Sabit bir \(m\) icin \(c=0\) terimi \(1\)'dir; sonraki her terim oncekinin \((m-c+1)/(N-c+1)\) ile carpilmasiyla elde edilir. Bu terimlerin toplami ortak alt blok tablosunu verir.

\((a,b)\) ucgeni uzerinde ana tarama

\(R_m\) tablosu hazir olduktan sonra ana is, \(0\le a\le N\) ve \(0\le b\le N-a\) kosullarini saglayan tum \((a,b)\) noktalari uzerinden gecmektir. Her sabit \(a\) icin uygulama \(q(a,0)=1\) ile baslar ve ayriklik katsayisini

$$q(a,b+1)=q(a,b)\frac{N-a-b}{N-b}$$

kuralıyla gunceller. Her kafes noktasinda \(R_{N-a-b}\) tablodan okunur, \(q(a,b)R_{N-a-b}\) toplama eklenir ve sonunda her sey \(n^2\)'ye bolunur.

Sayisal dogruluk ve derlenmis surumlerin iyilestirmeleri

Hesaplama tamamen kayan noktalidir. Bu nedenle derlenmis uygulamalar birikim hatalarini azaltmak icin telafili toplama kullanir; C++ surumu ana toplami ek olarak genisletilmis hassasiyette tutar. Derlenmis surumler distaki \(a\) araligini bloklara bolup bunlari paralel isler; Python uygulamasi ayni matematigi seri olarak uygular.

Complexity Analysis / Karmaşıklık Analizi

\(N=n-1\) olmak uzere, tersler tablosunu kurmak \(O(N)\) zaman ve \(O(N)\) bellek ister. Uygulamalardaki dogrudan ic carpimlarla tum \(R_m\) degerlerini hesaplamak \(O(N^2)\) aritmetik islem gerektirir. \((a,b)\) uzerindeki ana ucgensel toplam da yine \(O(N^2)\)'dir.

Dolayisiyla toplam calisma suresi \(O(N^2)\), bellek kullanimi ise \(O(N)\)'dir. Paralellesme derlenmis surumlerde gercek calisma suresini dusurur ama asimptotik mertebeyi degistirmez.

Footnotes and References / Dipnotlar ve Kaynaklar

  1. Problem sayfasi: Project Euler 906
  2. Tam sira: Wikipedia - Total order
  3. Binom katsayisi: Wikipedia - Binomial coefficient
  4. Dusen ve yukselen faktoriyeller: Wikipedia - Falling and rising factorials
  5. Hipergeometrik dagilim: Wikipedia - Hypergeometric distribution
  6. Hockey-stick ozdesligi: Wikipedia - Hockey-stick identity
  7. Kahan toplama algoritmasi: Wikipedia - Kahan summation algorithm

Problem Summary / Resumen del Problema

El problema 906 pide una probabilidad \(P(n)\), y la evaluacion final usa \(n=20000\). Las implementaciones muestran que todo el calculo puede organizarse alrededor de una alternativa distinguida \(x\). Si \(x\) ocupa la posicion \(a+1\) en el primer orden de preferencia y la posicion \(b+1\) en el segundo, entonces entre las otras \(N=n-1\) alternativas hay \(a\) por encima de \(x\) en el primer orden y \(b\) por encima de \(x\) en el segundo.

Las demas alternativas se reparten en tres bloques: las que estan por encima de \(x\) solo en el primer orden, las que estan por encima de \(x\) solo en el segundo y las que quedan por debajo de \(x\) en ambos ordenes. El evento evaluado por las implementaciones exige que los dos bloques superiores sean disjuntos; despues aparece una contribucion auxiliar procedente del bloque inferior comun. Al promediar esa contribucion sobre los \(n^2\) pares de rangos equiprobables de \(x\), se obtiene la probabilidad buscada.

Mathematical Approach / Enfoque Matematico

Sea \(N=n-1\). Para una alternativa fija \(x\), las dos permutaciones generan un par de coordenadas de rango. La solucion convierte la pregunta probabilistica original en una suma doble limpia sobre esas coordenadas.

Fijar una alternativa y registrar sus dos rangos

Escribamos \(y\succ_1 x\) cuando el primer orden de preferencia coloca a \(y\) por encima de \(x\), y \(y\succ_2 x\) para el segundo orden. Definimos

$$A=\{y:y\succ_1 x\},\qquad B=\{y:y\succ_2 x\}.$$

Si \(|A|=a\) y \(|B|=b\), entonces \(x\) tiene rango \(a+1\) en el primer orden y rango \(b+1\) en el segundo. Como los dos ordenes son independientes y cada rango de \(x\) es uniforme en \(\{1,\dots,n\}\), cada par \((a,b)\in\{0,\dots,N\}^2\) aparece con probabilidad \(1/n^2\).

Por que los conjuntos superiores deben ser disjuntos

La formula implementada solo aporta cuando ninguna alternativa queda por encima de \(x\) en ambos ordenes a la vez. En lenguaje de conjuntos, la condicion es

$$A\cap B=\varnothing.$$

Para tamaños fijos \(a\) y \(b\), el conjunto \(B\) es un subconjunto aleatorio uniforme de tamaño \(b\) entre las \(N\) alternativas distintas de \(x\). Exactamente \(\binom{N-a}{b}\) de esos subconjuntos evitan los \(a\) elementos ya contenidos en \(A\). Por tanto

$$q(a,b)=\Pr(A\cap B=\varnothing\mid |A|=a,\ |B|=b)=\frac{\binom{N-a}{b}}{\binom{N}{b}}.$$

Esto explica de inmediato la region triangular de sumacion en el codigo: si \(a+b>N\), no pueden existir conjuntos disjuntos con esos tamaños, de modo que la contribucion es nula. Tambien explica la actualizacion multiplicativa del bucle interior:

$$q(a,0)=1,\qquad q(a,b+1)=q(a,b)\frac{N-a-b}{N-b}.$$

El bloque inferior unanime

Una vez que \(A\) y \(B\) son disjuntos, toda alternativa restante queda por debajo de \(x\) en ambos ordenes. Llamemos a ese bloque inferior comun

$$C=\{y:x\succ_1 y\ \text{and}\ x\succ_2 y\},\qquad |C|=m=N-a-b.$$

Las implementaciones asocian a este bloque la cantidad auxiliar

$$R_m=\sum_{c=0}^{m}\prod_{t=0}^{c-1}\frac{m-t}{N-t}=\sum_{c=0}^{m}\frac{(m)_c}{(N)_c}=\sum_{c=0}^{m}\frac{\binom{m}{c}}{\binom{N}{c}}.$$

El termino con un \(c\) fijo tiene una interpretacion directa: \(\binom{m}{c}/\binom{N}{c}\) es la probabilidad de que un \(c\)-subconjunto uniforme de las \(N\) alternativas distintas de \(x\) quede completamente contenido en el bloque inferior comun \(C\). El valor \(R_m\) suma esas probabilidades para todos los tamaños \(c\) posibles.

Reducir la suma auxiliar a una forma cerrada

La suma directa ya sirve para el calculo, pero la estructura matematica se ve mejor si se simplifica. Partimos de

$$\frac{\binom{m}{c}}{\binom{N}{c}}=\frac{1}{\binom{N}{m}}\binom{N-c}{m-c}.$$

Al sustituir en la formula anterior se obtiene

$$R_m=\frac{1}{\binom{N}{m}}\sum_{c=0}^{m}\binom{N-c}{m-c}.$$

Si escribimos \(j=m-c\), la suma pasa a ser

$$\sum_{j=0}^{m}\binom{N-m+j}{j},$$

que es una suma clasica de tipo hockey-stick. Entonces

$$R_m=\frac{\binom{N+1}{m}}{\binom{N}{m}}=\frac{N+1}{N+1-m}.$$

Como \(m=N-a-b\), esto tambien se escribe como

$$R_{N-a-b}=\frac{n}{a+b+1}.$$

La formula final para la probabilidad

Al unir ambas piezas, la cantidad que calculan las implementaciones en C++, Python y Java es

$$\boxed{P(n)=\frac{1}{n^2}\sum_{a=0}^{N}\sum_{b=0}^{N-a} q(a,b)\,R_{N-a-b}=\frac{1}{n^2}\sum_{a=0}^{N}\sum_{b=0}^{N-a}\frac{\binom{N-a}{b}}{\binom{N}{b}}\cdot\frac{n}{a+b+1}.}$$

En otras palabras, el programa no hace mas que evaluar con rapidez esta suma doble triangular sin construir de forma explicita factoriales enormes ni coeficientes binomiales.

Ejemplo trabajado: \(n=3\)

Tomemos \(n=3\), de modo que \(N=2\). Entonces

$$R_0=1,\qquad R_1=1+\frac12=\frac32,\qquad R_2=1+1+1=3.$$

Los pares admisibles \((a,b)\) cumplen \(a+b\le 2\), y sus contribuciones son

$$\begin{aligned} (0,0)&:\ q=1,\ m=2,\ qR_m=3,\\ (0,1)&:\ q=1,\ m=1,\ qR_m=\frac32,\\ (0,2)&:\ q=1,\ m=0,\ qR_m=1,\\ (1,0)&:\ q=1,\ m=1,\ qR_m=\frac32,\\ (1,1)&:\ q=\frac12,\ m=0,\ qR_m=\frac12,\\ (2,0)&:\ q=1,\ m=0,\ qR_m=1. \end{aligned}$$

La suma total es \(17/2\). Al dividir por \(n^2=9\), resulta

$$P(3)=\frac{17/2}{9}=\frac{17}{18}.$$

Este caso pequeno ya muestra todas las piezas de la solucion general: el factor de disjuncion, la suma del bloque inferior y el promedio final sobre los pares de rangos.

How the Code Works / Como Funciona el Codigo

Precalculo de recíprocos y de los valores del bloque inferior

Las implementaciones construyen primero una tabla de recíprocos \(1,2,\dots,N\). Asi, cada cociente de la suma doble puede actualizarse por multiplicacion en lugar de repetir divisiones. Despues precalculan todos los valores \(R_m\) para \(m=0,\dots,N\). Para un \(m\) fijo, el termino con \(c=0\) vale \(1\), y cada termino siguiente se obtiene multiplicando el anterior por \((m-c+1)/(N-c+1)\). La suma de esos terminos produce la tabla del bloque inferior comun.

El barrido triangular sobre \((a,b)\)

Una vez disponible la tabla \(R_m\), el trabajo principal es el doble bucle sobre todos los \(a\) y \(b\) con \(0\le a\le N\) y \(0\le b\le N-a\). Para cada \(a\) fijo, la implementacion empieza con \(q(a,0)=1\) y actualiza el factor de disjuncion mediante

$$q(a,b+1)=q(a,b)\frac{N-a-b}{N-b}.$$

En cada punto \((a,b)\), toma \(R_{N-a-b}\), suma \(q(a,b)R_{N-a-b}\) al acumulado y al final divide todo por \(n^2\).

Precision numerica y optimizaciones en los lenguajes compilados

Todo el calculo es de punto flotante. Por eso las implementaciones compiladas usan suma compensada para reducir el error de acumulacion, y la version en C++ ademas mantiene el total principal en precision extendida. Las versiones compiladas dividen el rango exterior de \(a\) en bloques y procesan esos bloques en paralelo, mientras que la implementacion en Python sigue exactamente la misma matematica de forma secuencial.

Complexity Analysis / Analisis de Complejidad

Con \(N=n-1\), construir la tabla de recíprocos cuesta \(O(N)\) tiempo y \(O(N)\) memoria. Precalcular todos los valores \(R_m\) con los productos internos directos usados en las implementaciones requiere \(O(N^2)\) operaciones aritmeticas. La suma triangular principal sobre \((a,b)\) tambien es \(O(N^2)\).

Por tanto, el tiempo total es \(O(N^2)\) y el consumo de memoria es \(O(N)\). La paralelizacion mejora el tiempo real en las versiones compiladas, pero no cambia el orden asintotico.

Footnotes and References / Notas y Referencias

  1. Pagina del problema: Project Euler 906
  2. Orden total: Wikipedia - Total order
  3. Coeficiente binomial: Wikipedia - Binomial coefficient
  4. Factoriales descendentes y ascendentes: Wikipedia - Falling and rising factorials
  5. Distribucion hipergeometrica: Wikipedia - Hypergeometric distribution
  6. Identidad hockey-stick: Wikipedia - Hockey-stick identity
  7. Algoritmo de suma de Kahan: Wikipedia - Kahan summation algorithm

Problem Summary / 问题概述

第 906 题要求计算一个概率 \(P(n)\),最终取值是 \(n=20000\)。从实现可以看出,整个计算都围绕某个固定候选项 \(x\) 展开。如果 \(x\) 在第一份偏好排序中的位置是第 \(a+1\) 名,在第二份偏好排序中的位置是第 \(b+1\) 名,那么在其余 \(N=n-1\) 个候选项里,就有 \(a\) 个在第一份排序中位于 \(x\) 之上,有 \(b\) 个在第二份排序中位于 \(x\) 之上。

于是其余候选项自然分成三个部分:只在第一份排序里高于 \(x\) 的、只在第二份排序里高于 \(x\) 的,以及在两份排序里都低于 \(x\) 的。实现所计算的事件要求前两个“上方集合”互不相交,然后再乘上一个来自共同下方块的辅助贡献。最后再对 \(x\) 的 \(n^2\) 个等概率秩对取平均,就得到题目需要的概率。

Mathematical Approach / 数学方法

记 \(N=n-1\)。对一个固定的候选项 \(x\) 来说,两份随机排序给出一对秩坐标。解法的核心,就是把原来的概率问题改写成这对坐标上的一个三角形双重求和。

固定一个候选项,并记录它在两份排序中的位置

若第一份排序把 \(y\) 放在 \(x\) 前面,记作 \(y\succ_1 x\);第二份排序中类似地记作 \(y\succ_2 x\)。定义

$$A=\{y:y\succ_1 x\},\qquad B=\{y:y\succ_2 x\}.$$

当 \(|A|=a\)、\(|B|=b\) 时,\(x\) 在第一份排序中的名次是 \(a+1\),在第二份排序中的名次是 \(b+1\)。由于两份排序彼此独立,而且 \(x\) 在每一份排序中的名次都在 \(\{1,\dots,n\}\) 上均匀分布,因此每个 \((a,b)\in\{0,\dots,N\}^2\) 都以 \(1/n^2\) 的概率出现。

为什么上方集合必须互不相交

实现中的贡献只在一种情况下非零:不能存在某个候选项在两份排序里都排在 \(x\) 前面。用集合语言写就是

$$A\cap B=\varnothing.$$

在 \(a\) 和 \(b\) 固定时,集合 \(B\) 是从除 \(x\) 之外的 \(N\) 个候选项中等概率选出的一个大小为 \(b\) 的子集。其中恰好有 \(\binom{N-a}{b}\) 个子集能够避开已经落入 \(A\) 的那 \(a\) 个元素。因此

$$q(a,b)=\Pr(A\cap B=\varnothing\mid |A|=a,\ |B|=b)=\frac{\binom{N-a}{b}}{\binom{N}{b}}.$$

这也直接解释了代码中的三角形求和区域:如果 \(a+b>N\),那么这样大小的两个不相交集合根本不可能存在,所以贡献必为零。同时,这个公式还给出内层循环中的乘法递推:

$$q(a,0)=1,\qquad q(a,b+1)=q(a,b)\frac{N-a-b}{N-b}.$$

共同下方块

一旦 \(A\) 与 \(B\) 不相交,剩下的所有候选项就在两份排序里都排在 \(x\) 后面。把这个共同下方块记为

$$C=\{y:x\succ_1 y\ \text{and}\ x\succ_2 y\},\qquad |C|=m=N-a-b.$$

实现对这个块配套使用了一个辅助量

$$R_m=\sum_{c=0}^{m}\prod_{t=0}^{c-1}\frac{m-t}{N-t}=\sum_{c=0}^{m}\frac{(m)_c}{(N)_c}=\sum_{c=0}^{m}\frac{\binom{m}{c}}{\binom{N}{c}}.$$

其中固定 \(c\) 的那一项有很直接的意义:\(\binom{m}{c}/\binom{N}{c}\) 就是从除 \(x\) 之外的 \(N\) 个候选项里等概率选出一个 \(c\) 元子集,而这个子集完全落在共同下方块 \(C\) 中的概率。于是 \(R_m\) 就是把所有可能的 \(c\) 都累加起来。

把辅助和式化简成闭式

直接求和已经足够用来计算,但把它化简之后,数学结构会更清楚。先从恒等式

$$\frac{\binom{m}{c}}{\binom{N}{c}}=\frac{1}{\binom{N}{m}}\binom{N-c}{m-c}$$

开始。代入上式得到

$$R_m=\frac{1}{\binom{N}{m}}\sum_{c=0}^{m}\binom{N-c}{m-c}.$$

令 \(j=m-c\),则求和变成

$$\sum_{j=0}^{m}\binom{N-m+j}{j},$$

这正是一个标准的 hockey-stick 恒等式,因此

$$R_m=\frac{\binom{N+1}{m}}{\binom{N}{m}}=\frac{N+1}{N+1-m}.$$

再把 \(m=N-a-b\) 代回去,就得到更适合最终公式的写法

$$R_{N-a-b}=\frac{n}{a+b+1}.$$

最终概率公式

把上方集合的概率因子与下方块的辅助和结合起来,C++、Python 与 Java 实现计算的就是

$$\boxed{P(n)=\frac{1}{n^2}\sum_{a=0}^{N}\sum_{b=0}^{N-a} q(a,b)\,R_{N-a-b}=\frac{1}{n^2}\sum_{a=0}^{N}\sum_{b=0}^{N-a}\frac{\binom{N-a}{b}}{\binom{N}{b}}\cdot\frac{n}{a+b+1}.}$$

因此程序真正做的事,只是高效地计算这个三角形双重求和,而不去显式构造巨大的阶乘或二项式系数。

算例:\(n=3\)

取 \(n=3\),则 \(N=2\)。这时

$$R_0=1,\qquad R_1=1+\frac12=\frac32,\qquad R_2=1+1+1=3.$$

可行的 \((a,b)\) 满足 \(a+b\le 2\),对应贡献为

$$\begin{aligned} (0,0)&:\ q=1,\ m=2,\ qR_m=3,\\ (0,1)&:\ q=1,\ m=1,\ qR_m=\frac32,\\ (0,2)&:\ q=1,\ m=0,\ qR_m=1,\\ (1,0)&:\ q=1,\ m=1,\ qR_m=\frac32,\\ (1,1)&:\ q=\frac12,\ m=0,\ qR_m=\frac12,\\ (2,0)&:\ q=1,\ m=0,\ qR_m=1. \end{aligned}$$

这些项加起来是 \(17/2\)。再除以 \(n^2=9\),得到

$$P(3)=\frac{17/2}{9}=\frac{17}{18}.$$

这个小例子已经完整展示了通解的全部部件:上方集合的“不相交”因子、共同下方块的辅助和,以及最后对秩对的平均。

How the Code Works / 代码如何工作

预处理倒数表和下方块表

实现首先建立 \(1,2,\dots,N\) 的倒数表。这样一来,双重求和里出现的各种比值都可以通过连乘来更新,而不需要在循环中反复做除法。随后程序预先计算所有 \(R_m\) 的值。对固定的 \(m\),当 \(c=0\) 时项为 \(1\),之后每一项都由前一项乘上 \((m-c+1)/(N-c+1)\) 得到。把这些项累加起来,就得到共同下方块所需的整张表。

在 \((a,b)\) 三角区域上进行主循环

有了 \(R_m\) 之后,主计算就是遍历所有满足 \(0\le a\le N\) 且 \(0\le b\le N-a\) 的 \((a,b)\)。对每个固定的 \(a\),实现从 \(q(a,0)=1\) 开始,并按照

$$q(a,b+1)=q(a,b)\frac{N-a-b}{N-b}$$

递推更新“不相交”因子。每到一个格点,就查表取出 \(R_{N-a-b}\),把 \(q(a,b)R_{N-a-b}\) 加到总和中,最后再统一除以 \(n^2\)。

数值精度与编译型实现的优化

整个计算都是浮点运算。因此,编译型实现使用补偿求和来减小累计误差,而 C++ 版本还把主累加保存在扩展精度中。编译型版本还会把外层的 \(a\) 区间切成若干块并行处理;Python 实现则按完全相同的数学公式串行执行。

Complexity Analysis / 复杂度分析

记 \(N=n-1\)。建立倒数表需要 \(O(N)\) 时间和 \(O(N)\) 空间。用实现中的直接内层乘积来预计算全部 \(R_m\) 需要 \(O(N^2)\) 次算术运算。主三角形求和对 \((a,b)\) 的遍历同样是 \(O(N^2)\)。

因此,总时间复杂度是 \(O(N^2)\),空间复杂度是 \(O(N)\)。并行化只会改善编译型版本的实际运行时间,不会改变渐近量级。

Footnotes and References / 注释与参考资料

  1. 题目页面:Project Euler 906
  2. 全序:Wikipedia - Total order
  3. 二项式系数:Wikipedia - Binomial coefficient
  4. 下降阶乘与上升阶乘:Wikipedia - Falling and rising factorials
  5. 超几何分布:Wikipedia - Hypergeometric distribution
  6. Hockey-stick 恒等式:Wikipedia - Hockey-stick identity
  7. Kahan 求和算法:Wikipedia - Kahan summation algorithm

Problem Summary / Краткое описание задачи

В задаче 906 требуется найти вероятность \(P(n)\); в финальном вычислении используется \(n=20000\). Из реализаций видно, что всю схему удобно строить вокруг одной фиксированной альтернативы \(x\). Если \(x\) стоит на позиции \(a+1\) в первом порядке предпочтения и на позиции \(b+1\) во втором, то среди остальных \(N=n-1\) альтернатив ровно \(a\) находятся выше \(x\) в первом порядке и ровно \(b\) находятся выше \(x\) во втором.

Остальные альтернативы распадаются на три блока: выше \(x\) только в первом порядке, выше \(x\) только во втором и ниже \(x\) сразу в обоих порядках. Событие, которое учитывают реализации, требует, чтобы два верхних блока не пересекались, после чего появляется дополнительный вклад от общего нижнего блока. Усреднение этого вклада по всем \(n^2\) равновероятным парам рангов \(x\) и дает искомую вероятность.

Mathematical Approach / Математический подход

Положим \(N=n-1\). Для фиксированной альтернативы \(x\) две случайные перестановки задают пару ранговых координат. Решение сводит исходную вероятностную задачу к аккуратной двойной сумме по этим координатам.

Фиксируем одну альтернативу и записываем ее два ранга

Обозначим через \(y\succ_1 x\) факт, что первый порядок ставит \(y\) выше \(x\), а через \(y\succ_2 x\) аналогичный факт для второго порядка. Тогда

$$A=\{y:y\succ_1 x\},\qquad B=\{y:y\succ_2 x\}.$$

Если \(|A|=a\) и \(|B|=b\), то \(x\) имеет ранг \(a+1\) в первом порядке и \(b+1\) во втором. Поскольку два порядка независимы и каждый ранг \(x\) равномерно распределен по \(\{1,\dots,n\}\), любая пара \((a,b)\in\{0,\dots,N\}^2\) возникает с вероятностью \(1/n^2\).

Почему верхние множества должны быть непересекающимися

Реализованная формула дает вклад только тогда, когда не существует альтернативы, стоящей выше \(x\) сразу в обоих порядках. На языке множеств это условие записывается так:

$$A\cap B=\varnothing.$$

При фиксированных размерах \(a\) и \(b\) множество \(B\) является равновероятным \(b\)-подмножеством среди \(N\) альтернатив, отличных от \(x\). Ровно \(\binom{N-a}{b}\) таких подмножеств избегают \(a\) элементов, уже лежащих в \(A\). Следовательно,

$$q(a,b)=\Pr(A\cap B=\varnothing\mid |A|=a,\ |B|=b)=\frac{\binom{N-a}{b}}{\binom{N}{b}}.$$

Этим сразу объясняется и треугольная область суммирования в коде: если \(a+b>N\), то непересекающихся множеств таких размеров быть не может, и вклад равен нулю. Отсюда же получается мультипликативное обновление во внутреннем цикле:

$$q(a,0)=1,\qquad q(a,b+1)=q(a,b)\frac{N-a-b}{N-b}.$$

Общий нижний блок

Если \(A\) и \(B\) не пересекаются, то каждая оставшаяся альтернатива стоит ниже \(x\) в обоих порядках. Обозначим общий нижний блок через

$$C=\{y:x\succ_1 y\ \text{and}\ x\succ_2 y\},\qquad |C|=m=N-a-b.$$

Реализации связывают с этим блоком вспомогательную величину

$$R_m=\sum_{c=0}^{m}\prod_{t=0}^{c-1}\frac{m-t}{N-t}=\sum_{c=0}^{m}\frac{(m)_c}{(N)_c}=\sum_{c=0}^{m}\frac{\binom{m}{c}}{\binom{N}{c}}.$$

Слагаемое при фиксированном \(c\) имеет прямую интерпретацию: \(\binom{m}{c}/\binom{N}{c}\) равно вероятности того, что равномерно выбранное \(c\)-подмножество из \(N\) альтернатив, отличных от \(x\), целиком лежит внутри общего нижнего блока \(C\). Значит, \(R_m\) просто суммирует эти вероятности по всем допустимым размерам \(c\).

Сведение вспомогательной суммы к замкнутой форме

Для вычислений прямой суммы уже достаточно, но математическая картина становится прозрачнее после упрощения. Начнем с тождества

$$\frac{\binom{m}{c}}{\binom{N}{c}}=\frac{1}{\binom{N}{m}}\binom{N-c}{m-c}.$$

Подстановка дает

$$R_m=\frac{1}{\binom{N}{m}}\sum_{c=0}^{m}\binom{N-c}{m-c}.$$

Если положить \(j=m-c\), то сумма принимает вид

$$\sum_{j=0}^{m}\binom{N-m+j}{j},$$

а это стандартная hockey-stick-сумма. Поэтому

$$R_m=\frac{\binom{N+1}{m}}{\binom{N}{m}}=\frac{N+1}{N+1-m}.$$

С учетом \(m=N-a-b\) получаем особенно удобную запись

$$R_{N-a-b}=\frac{n}{a+b+1}.$$

Итоговая формула вероятности

Объединяя два множителя, получаем величину, которую вычисляют реализации на C++, Python и Java:

$$\boxed{P(n)=\frac{1}{n^2}\sum_{a=0}^{N}\sum_{b=0}^{N-a} q(a,b)\,R_{N-a-b}=\frac{1}{n^2}\sum_{a=0}^{N}\sum_{b=0}^{N-a}\frac{\binom{N-a}{b}}{\binom{N}{b}}\cdot\frac{n}{a+b+1}.}$$

Тем самым программа на самом деле лишь быстро вычисляет эту треугольную двойную сумму, не строя явно большие факториалы и биномиальные коэффициенты.

Разобранный пример: \(n=3\)

Возьмем \(n=3\), тогда \(N=2\). Имеем

$$R_0=1,\qquad R_1=1+\frac12=\frac32,\qquad R_2=1+1+1=3.$$

Допустимые пары \((a,b)\) удовлетворяют условию \(a+b\le 2\), а их вклады равны

$$\begin{aligned} (0,0)&:\ q=1,\ m=2,\ qR_m=3,\\ (0,1)&:\ q=1,\ m=1,\ qR_m=\frac32,\\ (0,2)&:\ q=1,\ m=0,\ qR_m=1,\\ (1,0)&:\ q=1,\ m=1,\ qR_m=\frac32,\\ (1,1)&:\ q=\frac12,\ m=0,\ qR_m=\frac12,\\ (2,0)&:\ q=1,\ m=0,\ qR_m=1. \end{aligned}$$

Сумма равна \(17/2\). После деления на \(n^2=9\) получаем

$$P(3)=\frac{17/2}{9}=\frac{17}{18}.$$

Этот небольшой пример уже показывает все ключевые части общей схемы: фактор непересечения, вспомогательную сумму по общему нижнему блоку и финальное усреднение по ранговым парам.

How the Code Works / Как работает код

Предвычисление обратных величин и таблицы нижнего блока

Сначала реализации строят таблицу обратных чисел \(1,2,\dots,N\). Это позволяет обновлять все отношения в двойной сумме умножением, а не повторяющимся делением. Затем предвычисляются все значения \(R_m\) для \(m=0,\dots,N\). При фиксированном \(m\) член при \(c=0\) равен \(1\), а каждый следующий получается из предыдущего умножением на \((m-c+1)/(N-c+1)\). Сумма этих членов и образует нужную таблицу для общего нижнего блока.

Треугольный проход по \((a,b)\)

После подготовки таблицы \(R_m\) главная работа сводится к проходу по всем \(a\) и \(b\), удовлетворяющим \(0\le a\le N\) и \(0\le b\le N-a\). Для каждого фиксированного \(a\) реализация начинает с \(q(a,0)=1\) и обновляет фактор непересечения по формуле

$$q(a,b+1)=q(a,b)\frac{N-a-b}{N-b}.$$

В каждой точке \((a,b)\) из таблицы берется \(R_{N-a-b}\), затем \(q(a,b)R_{N-a-b}\) добавляется к накопленной сумме, а в самом конце все делится на \(n^2\).

Численная точность и оптимизации в компилируемых версиях

Все вычисления выполняются в плавающей точке. Поэтому компилируемые реализации используют компенсированное суммирование, чтобы уменьшить ошибку накопления, а версия на C++ дополнительно хранит основной итог в расширенной точности. Кроме того, компилируемые версии разбивают внешний диапазон по \(a\) на блоки и обрабатывают их параллельно, тогда как Python использует ту же математику последовательно.

Complexity Analysis / Анализ сложности

При \(N=n-1\) построение таблицы обратных чисел требует \(O(N)\) времени и \(O(N)\) памяти. Предвычисление всех значений \(R_m\) с помощью прямых внутренних произведений из реализаций стоит \(O(N^2)\) арифметических операций. Основная треугольная сумма по \((a,b)\) имеет ту же сложность \(O(N^2)\).

Следовательно, общая временная сложность равна \(O(N^2)\), а потребление памяти равно \(O(N)\). Параллелизм уменьшает фактическое время работы компилируемых версий, но не меняет асимптотику.

Footnotes and References / Сноски и ссылки

  1. Страница задачи: Project Euler 906
  2. Полный порядок: Wikipedia - Total order
  3. Биномиальный коэффициент: Wikipedia - Binomial coefficient
  4. Падающие и возрастающие факториалы: Wikipedia - Falling and rising factorials
  5. Гипергеометрическое распределение: Wikipedia - Hypergeometric distribution
  6. Тождество hockey-stick: Wikipedia - Hockey-stick identity
  7. Алгоритм суммирования Кахана: Wikipedia - Kahan summation algorithm

Problem Summary / ملخص المسألة

تطلب المسألة 906 حساب احتمال \(P(n)\)، وفي التقييم النهائي يكون \(n=20000\). وتُظهر التطبيقات أن الحساب كله يمكن تنظيمه حول بديل ثابت واحد نرمز له بـ \(x\). فإذا كان \(x\) في المرتبة \(a+1\) في ترتيب التفضيل الأول، وفي المرتبة \(b+1\) في ترتيب التفضيل الثاني، فمعنى ذلك أن بين البدائل الأخرى وعددها \(N=n-1\) يوجد \(a\) بدائل فوق \(x\) في الترتيب الأول، و\(b\) بدائل فوقه في الترتيب الثاني.

وبذلك تنقسم بقية البدائل إلى ثلاثة أقسام: بدائل فوق \(x\) في الترتيب الأول فقط، وبدائل فوقه في الترتيب الثاني فقط، وبدائل تقع تحته في الترتيبين معًا. والحدث الذي تحسبه التطبيقات يشترط أن تكون المجموعتان العلويتان غير متقاطعتين، ثم يضيف مساهمةً مساعدة آتية من الكتلة السفلية المشتركة. وبعد أخذ متوسط هذه المساهمة على جميع أزواج الرتب المتساوية الاحتمال وعددها \(n^2\) للبديل \(x\)، نحصل على الاحتمال المطلوب.

Mathematical Approach / المنهج الرياضي

لنضع \(N=n-1\). فعندما نثبّت بديلًا واحدًا \(x\)، فإن الترتيبين العشوائيين يحددان له زوجًا من إحداثيات الرتبة. وفكرة الحل هي تحويل السؤال الاحتمالي الأصلي إلى مجموع ثنائي مثلثي على هذه الإحداثيات.

تثبيت بديل واحد وتسجيل مرتبتيه

نكتب \(y\succ_1 x\) إذا كان ترتيب التفضيل الأول يضع \(y\) فوق \(x\)، ونكتب \(y\succ_2 x\) للترتيب الثاني. وعندئذ نعرّف

$$A=\{y:y\succ_1 x\},\qquad B=\{y:y\succ_2 x\}.$$

إذا كان \(|A|=a\) و\(|B|=b\)، فهذا يعني أن رتبة \(x\) هي \(a+1\) في الترتيب الأول و\(b+1\) في الترتيب الثاني. وبما أن الترتيبين مستقلان، وأن رتبة \(x\) في كل منهما موزعة بانتظام على \(\{1,\dots,n\}\)، فإن كل زوج \((a,b)\in\{0,\dots,N\}^2\) يظهر باحتمال \(1/n^2\).

لماذا يجب أن تكون المجموعتان العلويتان متباعدتين

الصيغة التي تنفذها التطبيقات لا تعطي مساهمة إلا إذا لم يوجد بديل واحد يتفوق على \(x\) في الترتيبين معًا. وبصياغة المجموعات يكون الشرط

$$A\cap B=\varnothing.$$

عند تثبيت الحجمين \(a\) و\(b\)، تكون المجموعة \(B\) مجموعةً جزئية عشوائية منتظمة حجمها \(b\) من بين \(N\) بديلًا غير \(x\). ويوجد بالضبط \(\binom{N-a}{b}\) من هذه المجموعات الجزئية يتجنب العناصر \(a\) الموجودة أصلًا في \(A\). ومن ثم

$$q(a,b)=\Pr(A\cap B=\varnothing\mid |A|=a,\ |B|=b)=\frac{\binom{N-a}{b}}{\binom{N}{b}}.$$

وهذا يفسر مباشرةً لماذا يسير الجمع في الكود داخل منطقة مثلثية: فإذا كان \(a+b>N\)، فلا يمكن أصلًا أن توجد مجموعتان غير متقاطعتين بهذين الحجمين، فيصبح الإسهام صفرًا. كما يفسر أيضًا التحديث الضربي داخل الحلقة الداخلية:

$$q(a,0)=1,\qquad q(a,b+1)=q(a,b)\frac{N-a-b}{N-b}.$$

الكتلة السفلية المشتركة

متى كان \(A\) و\(B\) غير متقاطعين، فإن كل بديل متبقٍ يقع تحت \(x\) في الترتيبين معًا. ولنرمز إلى هذه الكتلة السفلية المشتركة بـ

$$C=\{y:x\succ_1 y\ \text{and}\ x\succ_2 y\},\qquad |C|=m=N-a-b.$$

وتربط التطبيقات بهذه الكتلة الكمية المساعدة

$$R_m=\sum_{c=0}^{m}\prod_{t=0}^{c-1}\frac{m-t}{N-t}=\sum_{c=0}^{m}\frac{(m)_c}{(N)_c}=\sum_{c=0}^{m}\frac{\binom{m}{c}}{\binom{N}{c}}.$$

ولكل \(c\) ثابت، فإن الكسر \(\binom{m}{c}/\binom{N}{c}\) يساوي احتمال أن تقع مجموعة جزئية منتظمة من \(c\) عناصر، مختارة من بين \(N\) بديلًا غير \(x\)، بالكامل داخل الكتلة السفلية المشتركة \(C\). ولذلك فإن \(R_m\) يجمع هذه الاحتمالات على جميع القيم الممكنة لـ \(c\).

اختزال المجموع المساعد إلى صيغة مغلقة

المجموع المباشر يكفي للحساب، لكن الصورة الرياضية تصبح أوضح بعد تبسيطه. نبدأ من الهوية

$$\frac{\binom{m}{c}}{\binom{N}{c}}=\frac{1}{\binom{N}{m}}\binom{N-c}{m-c}.$$

وبالتعويض نحصل على

$$R_m=\frac{1}{\binom{N}{m}}\sum_{c=0}^{m}\binom{N-c}{m-c}.$$

وإذا وضعنا \(j=m-c\)، صار المجموع

$$\sum_{j=0}^{m}\binom{N-m+j}{j},$$

وهو مجموع قياسي من نوع hockey-stick، ومنه نستنتج

$$R_m=\frac{\binom{N+1}{m}}{\binom{N}{m}}=\frac{N+1}{N+1-m}.$$

وبما أن \(m=N-a-b\)، فيمكن كتابة ذلك بصورة أنسب للصيغة النهائية:

$$R_{N-a-b}=\frac{n}{a+b+1}.$$

الصيغة النهائية للاحتمال

عند جمع عامل عدم التقاطع مع مساهمة الكتلة السفلية المشتركة، فإن ما تحسبه تطبيقات C++ وPython وJava هو

$$\boxed{P(n)=\frac{1}{n^2}\sum_{a=0}^{N}\sum_{b=0}^{N-a} q(a,b)\,R_{N-a-b}=\frac{1}{n^2}\sum_{a=0}^{N}\sum_{b=0}^{N-a}\frac{\binom{N-a}{b}}{\binom{N}{b}}\cdot\frac{n}{a+b+1}.}$$

وبعبارة أخرى، فإن البرنامج ليس إلا طريقة فعالة لتقييم هذا المجموع الثنائي المثلثي من غير بناء العاملية الكبيرة أو معاملات ثنائية الحد صراحةً.

مثال محلول: \(n=3\)

لنأخذ \(n=3\)، وعندها \(N=2\). تصبح القيم

$$R_0=1,\qquad R_1=1+\frac12=\frac32,\qquad R_2=1+1+1=3.$$

والأزواج الممكنة \((a,b)\) هي التي تحقق \(a+b\le 2\)، وتكون مساهماتها

$$\begin{aligned} (0,0)&:\ q=1,\ m=2,\ qR_m=3,\\ (0,1)&:\ q=1,\ m=1,\ qR_m=\frac32,\\ (0,2)&:\ q=1,\ m=0,\ qR_m=1,\\ (1,0)&:\ q=1,\ m=1,\ qR_m=\frac32,\\ (1,1)&:\ q=\frac12,\ m=0,\ qR_m=\frac12,\\ (2,0)&:\ q=1,\ m=0,\ qR_m=1. \end{aligned}$$

مجموع هذه الحدود يساوي \(17/2\). وبعد القسمة على \(n^2=9\) نحصل على

$$P(3)=\frac{17/2}{9}=\frac{17}{18}.$$

ويعرض هذا المثال الصغير بالفعل جميع مكونات الحل العام: عامل عدم التقاطع، ومجموع الكتلة السفلية، ثم المتوسط النهائي على أزواج الرتب.

How the Code Works / كيف يعمل الكود

التهيئة المسبقة للمقلوبات وقيم الكتلة السفلية

تبني التطبيقات أولًا جدولًا للمقلوبات \(1,2,\dots,N\). وبهذا يمكن تحديث كل نسبة في المجموع الثنائي بالضرب بدلًا من تكرار القسمة داخل الحلقات. ثم تُحسب جميع قيم \(R_m\) مسبقًا من أجل \(m=0,\dots,N\). فعند تثبيت \(m\)، يكون الحد الموافق لـ \(c=0\) مساويًا لـ \(1\)، ويُستخرج كل حد تالٍ من الحد السابق بضربه في \((m-c+1)/(N-c+1)\). ومجموع هذه الحدود يعطي الجدول المطلوب للكتلة السفلية المشتركة.

المسح المثلثي على \((a,b)\)

بعد تجهيز جدول \(R_m\)، تصبح المهمة الرئيسة هي المرور على جميع الأزواج \((a,b)\) التي تحقق \(0\le a\le N\) و\(0\le b\le N-a\). ولكل \(a\) ثابت، تبدأ التطبيقات من \(q(a,0)=1\)، ثم تحدّث عامل عدم التقاطع وفق

$$q(a,b+1)=q(a,b)\frac{N-a-b}{N-b}.$$

وعند كل نقطة شبكية، يُستخرج \(R_{N-a-b}\) من الجدول، ثم تضاف الكمية \(q(a,b)R_{N-a-b}\) إلى المجموع التراكمي، وفي النهاية يُقسم كل شيء على \(n^2\).

الدقة العددية وتحسينات اللغات المترجمة

كل الحسابات هنا بنقطة عائمة. لذلك تستخدم التطبيقات المترجمة جمعًا تعويضيًا لتقليل خطأ التراكم، كما تحتفظ نسخة C++ بالمجموع الرئيسي بدقة ممتدة. وتقوم النسخ المترجمة أيضًا بتقسيم المجال الخارجي الخاص بـ \(a\) إلى كتل تُعالج على نحو متوازٍ، في حين تنفذ نسخة Python الرياضيات نفسها على نحو تسلسلي.

Complexity Analysis / تحليل التعقيد

إذا كتبنا \(N=n-1\)، فإن بناء جدول المقلوبات يكلف \(O(N)\) زمنًا و\(O(N)\) ذاكرة. أما حساب جميع قيم \(R_m\) بالطريقة المباشرة المستخدمة في التطبيقات فيحتاج إلى \(O(N^2)\) عملية حسابية. والمجموع الرئيسي المثلثي على \((a,b)\) هو أيضًا من الرتبة \(O(N^2)\).

ومن ثم فإن الزمن الكلي يساوي \(O(N^2)\)، واستهلاك الذاكرة يساوي \(O(N)\). والتنفيذ المتوازي يحسن الزمن الفعلي في النسخ المترجمة، لكنه لا يغيّر الرتبة التقاربية.

Footnotes and References / الحواشي والمراجع

  1. صفحة المسألة: Project Euler 906
  2. الترتيب الكلي: Wikipedia - Total order
  3. معامل ثنائي الحد: Wikipedia - Binomial coefficient
  4. العامليات الساقطة والصاعدة: Wikipedia - Falling and rising factorials
  5. التوزيع فوق الهندسي: Wikipedia - Hypergeometric distribution
  6. متطابقة hockey-stick: Wikipedia - Hockey-stick identity
  7. خوارزمية جمع Kahan: Wikipedia - Kahan summation algorithm