The implementations use an equivalent tracked-ball formulation of the urn process. At stage \(t\), let
$$m=n-t+2,\qquad P_t=km.$$
The active population has \(P_t\) balls, and \(X_t\) of them are the balls whose future influence must still be tracked. We draw \(2k\) balls uniformly without replacement, and \(B_t\) denotes how many tracked balls are removed in that draw. The quantity to evaluate is
$$S(n,k)=\mathbb{E}\left[\sum_{t=1}^{n} B_t^2\right].$$
The direct state distribution becomes too large when \(n\) is large, so the key idea is to evolve only the first two moments of \(X_t\). The sample checkpoint used by the implementation is \(S(2,2)=9.6\).
The full distribution of the process is unnecessary once we notice that both the step contribution and the next-state update depend only on the first two moments of the current tracked-ball count.
Before the draw at stage \(t\), there are \(X_t\) tracked balls inside a population of size \(P_t=km\), where \(m=n-t+2\). Conditional on \(X_t=x\), the number drawn from the tracked set is hypergeometric:
$$\Pr(B_t=b\mid X_t=x)=\frac{\binom{x}{b}\binom{km-x}{2k-b}}{\binom{km}{2k}}.$$
After the draw, the next stage keeps the undrawn tracked balls and adds one fresh block of \(k\) tracked balls, so for \(t<n\),
$$X_{t+1}=X_t+k-B_t,$$
with initial value
$$X_1=k.$$
This identity is the backbone of the whole method: once we know the moments of \(X_t\), we can compute both the expected contribution of stage \(t\) and the moments of \(X_{t+1}\).
For a hypergeometric random variable, the first factorial moments are standard:
$$\mathbb{E}[B_t\mid X_t=x]=\frac{2k}{km}x=\frac{2}{m}x,$$
$$\mathbb{E}[B_t(B_t-1)\mid X_t=x]=\frac{2k(2k-1)}{km(km-1)}x(x-1)=\frac{2(2k-1)}{m(km-1)}x(x-1).$$
Define
$$\alpha_m=\frac{2}{m},\qquad \gamma_m=\frac{2(2k-1)}{m(km-1)}.$$
Because \(B_t^2=B_t(B_t-1)+B_t\), we get
$$\mathbb{E}[B_t^2\mid X_t=x]=\gamma_m x^2+(\alpha_m-\gamma_m)x.$$
If we denote
$$\mu_t=\mathbb{E}[X_t],\qquad \nu_t=\mathbb{E}[X_t^2],$$
then the expected contribution of stage \(t\) is
$$\mathbb{E}[B_t^2]=\gamma_m\nu_t+(\alpha_m-\gamma_m)\mu_t.$$
Starting from \(X_{t+1}=X_t+k-B_t\), condition on \(X_t=x\):
$$\mathbb{E}[X_{t+1}\mid X_t=x]=x+k-\mathbb{E}[B_t\mid X_t=x]=x+k-\frac{2}{m}x.$$
Hence
$$\mathbb{E}[X_{t+1}\mid X_t=x]=\frac{m-2}{m}x+k.$$
Define
$$\rho_m=\frac{m-2}{m}.$$
Taking expectations gives the linear recurrence
$$\mu_{t+1}=\rho_m\mu_t+k,$$
with initial value \(\mu_1=k\). This already shows why the process is manageable: the next mean depends only on the current mean.
Now square the state transition:
$$X_{t+1}^2=(X_t+k-B_t)^2.$$
Conditioning on \(X_t=x\) and inserting the formulas from Step 2 yields
$$\mathbb{E}[X_{t+1}^2\mid X_t=x]=\mathbb{E}[(x+k-B_t)^2\mid X_t=x].$$
After collecting the coefficients of \(x^2\), \(x\), and the constant term, the result is
$$\mathbb{E}[X_{t+1}^2\mid X_t=x]=\sigma_m x^2+\left(2k\rho_m+\rho_m-\sigma_m\right)x+k^2,$$
where
$$\sigma_m=\frac{(m-2)(k(m-2)-1)}{m(km-1)}.$$
Taking expectations gives
$$\nu_{t+1}=\sigma_m\nu_t+\left(2k\rho_m+\rho_m-\sigma_m\right)\mu_t+k^2,$$
with initial value \(\nu_1=k^2\). Therefore the whole process closes on the pair \((\mu_t,\nu_t)\), and no larger state object is required.
At the first stage, \(m=3\), so the population has \(6\) balls and \(X_1=2\). We draw \(4\) balls. The hypergeometric distribution is
$$\Pr(B_1=0)=\frac{\binom{2}{0}\binom{4}{4}}{\binom{6}{4}}=\frac{1}{15},$$
$$\Pr(B_1=1)=\frac{\binom{2}{1}\binom{4}{3}}{\binom{6}{4}}=\frac{8}{15},$$
$$\Pr(B_1=2)=\frac{\binom{2}{2}\binom{4}{2}}{\binom{6}{4}}=\frac{6}{15}.$$
Therefore
$$\mathbb{E}[B_1^2]=0^2\cdot\frac{1}{15}+1^2\cdot\frac{8}{15}+2^2\cdot\frac{6}{15}=\frac{32}{15}.$$
Next, \(X_2=4-B_1\), so \(X_2\) takes values \(4,3,2\) with the same probabilities, and
$$\mathbb{E}[X_2^2]=4^2\cdot\frac{1}{15}+3^2\cdot\frac{8}{15}+2^2\cdot\frac{6}{15}=\frac{112}{15}.$$
At the second stage, \(m=2\), so the draw size \(2k\) equals the whole remaining population. Hence \(B_2=X_2\), and
$$\mathbb{E}[B_2^2]=\mathbb{E}[X_2^2]=\frac{112}{15}.$$
The final expectation is
$$S(2,2)=\frac{32}{15}+\frac{112}{15}=\frac{144}{15}=\frac{48}{5}=9.6,$$
which matches the checkpoint embedded in the implementation.
The C++, Python, and Java implementations all follow the same recurrence. They store only the current first moment \(\mu_t\), the current second moment \(\nu_t\), and the running total of \(\mathbb{E}[B_t^2]\). For each stage they compute \(m=n-t+2\), evaluate the contribution formula
$$\gamma_m\nu_t+(\alpha_m-\gamma_m)\mu_t,$$
add it to the total, and then update \(\mu_t\) and \(\nu_t\) with the linear recurrences above unless the last stage has been reached.
The Python and Java versions use high-precision decimal arithmetic, while the C++ version uses extended floating-point arithmetic. The C++ implementation also contains tiny-case self-checks that propagate the full hypergeometric distribution and compare it against the moment recurrence. After the loop, the final expectation is rounded to the nearest integer before being printed.
The production algorithm performs one constant-size update per stage, so it runs in \(O(n)\) time and uses \(O(1)\) memory. This is the decisive improvement over a full distributional dynamic program, whose state space grows with the possible tracked-ball counts. The small validation routine is much more expensive because it keeps the whole probability distribution, but it is used only on tiny test cases.
Die Implementierungen verwenden eine äquivalente Formulierung des Urnenprozesses mit verfolgten Kugeln. In Stufe \(t\) setzen wir
$$m=n-t+2,\qquad P_t=km.$$
Die aktive Population enthält also \(P_t\) Kugeln, und \(X_t\) davon sind diejenigen, deren weiterer Einfluss noch verfolgt werden muss. Es werden \(2k\) Kugeln gleichverteilt ohne Zurücklegen gezogen; \(B_t\) bezeichnet die Anzahl der dabei entfernten verfolgten Kugeln. Gesucht ist
$$S(n,k)=\mathbb{E}\left[\sum_{t=1}^{n} B_t^2\right].$$
Für große \(n\) wäre es zu teuer, die vollständige Zustandsverteilung mitzuschleppen. Deshalb speichert die Lösung nur den ersten und zweiten Moment von \(X_t\). Der Kontrollwert aus der Implementierung ist \(S(2,2)=9.6\).
Entscheidend ist, dass sowohl der Beitrag einer Stufe als auch der Übergang zur nächsten Stufe nur von den ersten beiden Momenten der aktuellen Zahl verfolgter Kugeln abhängen.
Vor der Ziehung in Stufe \(t\) liegen \(X_t\) verfolgte Kugeln in einer Population der Größe \(P_t=km\), wobei \(m=n-t+2\) ist. Unter der Bedingung \(X_t=x\) ist die Anzahl gezogener verfolgter Kugeln hypergeometrisch verteilt:
$$\Pr(B_t=b\mid X_t=x)=\frac{\binom{x}{b}\binom{km-x}{2k-b}}{\binom{km}{2k}}.$$
Nach der Ziehung bleiben die nicht gezogenen verfolgten Kugeln erhalten, und für die nächste Stufe kommt ein frischer Block von \(k\) verfolgten Kugeln hinzu. Für \(t<n\) gilt also
$$X_{t+1}=X_t+k-B_t,$$
mit Anfangswert
$$X_1=k.$$
Genau diese Identität macht die Kompression auf wenige Momente möglich.
Für eine hypergeometrische Zufallsvariable sind die ersten Fakultätsmomente bekannt:
$$\mathbb{E}[B_t\mid X_t=x]=\frac{2k}{km}x=\frac{2}{m}x,$$
$$\mathbb{E}[B_t(B_t-1)\mid X_t=x]=\frac{2k(2k-1)}{km(km-1)}x(x-1)=\frac{2(2k-1)}{m(km-1)}x(x-1).$$
Wir definieren
$$\alpha_m=\frac{2}{m},\qquad \gamma_m=\frac{2(2k-1)}{m(km-1)}.$$
Mit \(B_t^2=B_t(B_t-1)+B_t\) folgt
$$\mathbb{E}[B_t^2\mid X_t=x]=\gamma_m x^2+(\alpha_m-\gamma_m)x.$$
Schreiben wir
$$\mu_t=\mathbb{E}[X_t],\qquad \nu_t=\mathbb{E}[X_t^2],$$
dann ist der erwartete Beitrag der Stufe \(t\)
$$\mathbb{E}[B_t^2]=\gamma_m\nu_t+(\alpha_m-\gamma_m)\mu_t.$$
Aus \(X_{t+1}=X_t+k-B_t\) erhält man bedingt auf \(X_t=x\)
$$\mathbb{E}[X_{t+1}\mid X_t=x]=x+k-\mathbb{E}[B_t\mid X_t=x]=x+k-\frac{2}{m}x.$$
Also
$$\mathbb{E}[X_{t+1}\mid X_t=x]=\frac{m-2}{m}x+k.$$
Setze
$$\rho_m=\frac{m-2}{m}.$$
Nach Erwartungsbildung folgt die lineare Rekursion
$$\mu_{t+1}=\rho_m\mu_t+k,$$
mit Anfangswert \(\mu_1=k\).
Nun quadrieren wir den Zustandsübergang:
$$X_{t+1}^2=(X_t+k-B_t)^2.$$
Nach Bedingung auf \(X_t=x\) und Einsetzen der Formeln aus Schritt 2 ergibt sich
$$\mathbb{E}[X_{t+1}^2\mid X_t=x]=\mathbb{E}[(x+k-B_t)^2\mid X_t=x].$$
Durch Zusammenfassen der Koeffizienten erhält man
$$\mathbb{E}[X_{t+1}^2\mid X_t=x]=\sigma_m x^2+\left(2k\rho_m+\rho_m-\sigma_m\right)x+k^2,$$
wobei
$$\sigma_m=\frac{(m-2)(k(m-2)-1)}{m(km-1)}.$$
Damit gilt
$$\nu_{t+1}=\sigma_m\nu_t+\left(2k\rho_m+\rho_m-\sigma_m\right)\mu_t+k^2,$$
mit Anfangswert \(\nu_1=k^2\). Der gesamte Prozess schließt also auf dem Momentenpaar \((\mu_t,\nu_t)\).
In der ersten Stufe ist \(m=3\), also gibt es \(6\) Kugeln und \(X_1=2\). Es werden \(4\) Kugeln gezogen. Die Verteilung lautet
$$\Pr(B_1=0)=\frac{\binom{2}{0}\binom{4}{4}}{\binom{6}{4}}=\frac{1}{15},$$
$$\Pr(B_1=1)=\frac{\binom{2}{1}\binom{4}{3}}{\binom{6}{4}}=\frac{8}{15},$$
$$\Pr(B_1=2)=\frac{\binom{2}{2}\binom{4}{2}}{\binom{6}{4}}=\frac{6}{15}.$$
Daher
$$\mathbb{E}[B_1^2]=0^2\cdot\frac{1}{15}+1^2\cdot\frac{8}{15}+2^2\cdot\frac{6}{15}=\frac{32}{15}.$$
Außerdem ist \(X_2=4-B_1\), also nimmt \(X_2\) die Werte \(4,3,2\) mit denselben Wahrscheinlichkeiten an und
$$\mathbb{E}[X_2^2]=4^2\cdot\frac{1}{15}+3^2\cdot\frac{8}{15}+2^2\cdot\frac{6}{15}=\frac{112}{15}.$$
In der zweiten Stufe gilt \(m=2\), also entspricht die Ziehung der gesamten Restpopulation. Somit ist \(B_2=X_2\), also
$$\mathbb{E}[B_2^2]=\mathbb{E}[X_2^2]=\frac{112}{15}.$$
Am Ende erhält man
$$S(2,2)=\frac{32}{15}+\frac{112}{15}=\frac{144}{15}=\frac{48}{5}=9.6,$$
genau den Kontrollwert der Implementierung.
Die C++-, Python- und Java-Implementierungen verwenden dieselbe Rekursion. Sie speichern nur den aktuellen ersten Moment \(\mu_t\), den aktuellen zweiten Moment \(\nu_t\) und die laufende Summe der Beiträge \(\mathbb{E}[B_t^2]\). Für jede Stufe wird \(m=n-t+2\) berechnet, dann der Beitrag
$$\gamma_m\nu_t+(\alpha_m-\gamma_m)\mu_t$$
zum Ergebnis addiert, und anschließend werden \(\mu_t\) und \(\nu_t\) mit den linearen Übergangsformeln aktualisiert, sofern noch nicht die letzte Stufe erreicht ist.
Die Python- und Java-Versionen benutzen hochpräzise Dezimalarithmetik, während die C++-Version mit erweitertem Gleitkomma arbeitet. Dort gibt es zusätzlich winzige Selbsttests, die für kleine Fälle die vollständige hypergeometrische Verteilung fortpflanzen und mit der Momentenrekursion vergleichen. Nach der Schleife wird der Erwartungswert auf die nächste ganze Zahl gerundet und ausgegeben.
Der eigentliche Algorithmus führt pro Stufe nur ein Update fester Größe aus. Daher beträgt die Laufzeit \(O(n)\) und der Speicherbedarf \(O(1)\). Genau darin liegt der Gewinn gegenüber einer vollständigen Verteilungs-Dynamik, deren Zustandsraum mit den möglichen Zahlen verfolgter Kugeln wächst. Die kleine Validierungsroutine ist deutlich teurer, wird aber nur für winzige Testfälle benutzt.
Uygulamalar, urn sürecini eşdeğer bir “izlenen bilye” modeliyle ifade ediyor. \(t\). adımda
$$m=n-t+2,\qquad P_t=km$$
olsun. Aktif popülasyonda toplam \(P_t\) bilye vardır ve bunların \(X_t\) tanesi gelecek adımlar açısından takip edilmesi gereken bilyelerdir. Bu popülasyondan geri koymadan düzgün rastgele \(2k\) bilye çekilir; çekilen izlenen bilye sayısına \(B_t\) diyelim. Hesaplanmak istenen büyüklük
$$S(n,k)=\mathbb{E}\left[\sum_{t=1}^{n} B_t^2\right]$$
ifadesidir. \(n\) büyükken tüm olasılık dağılımını adım adım taşımak pahalı olur. Bu yüzden çözüm, yalnızca \(X_t\)'nin birinci ve ikinci momentlerini günceller. Uygulamadaki örnek kontrol değeri \(S(2,2)=9.6\) sonucudur.
Ana gözlem şudur: hem bir adımın beklenen katkısı hem de bir sonraki duruma geçiş, mevcut izlenen bilye sayısının yalnızca ilk iki momentine bağlıdır.
\(t\). adımdaki çekilişten önce, büyüklüğü \(P_t=km\) olan popülasyonun içinde \(X_t\) adet izlenen bilye vardır; burada \(m=n-t+2\). \(X_t=x\) koşulu altında çekilen izlenen bilye sayısı hipergeometrik dağılıma uyar:
$$\Pr(B_t=b\mid X_t=x)=\frac{\binom{x}{b}\binom{km-x}{2k-b}}{\binom{km}{2k}}.$$
Çekilişten sonra çekilmeyen izlenen bilyeler korunur ve sonraki adıma yeni bir \(k\) bilyelik izlenen blok eklenir. Bu nedenle \(t<n\) için
$$X_{t+1}=X_t+k-B_t,$$
başlangıçta da
$$X_1=k$$
olur. Tüm yöntemin omurgası bu eşitliktir.
Hipergeometrik değişken için ilk faktöriyel momentler standarttır:
$$\mathbb{E}[B_t\mid X_t=x]=\frac{2k}{km}x=\frac{2}{m}x,$$
$$\mathbb{E}[B_t(B_t-1)\mid X_t=x]=\frac{2k(2k-1)}{km(km-1)}x(x-1)=\frac{2(2k-1)}{m(km-1)}x(x-1).$$
Şunları tanımlayalım:
$$\alpha_m=\frac{2}{m},\qquad \gamma_m=\frac{2(2k-1)}{m(km-1)}.$$
\(B_t^2=B_t(B_t-1)+B_t\) olduğundan
$$\mathbb{E}[B_t^2\mid X_t=x]=\gamma_m x^2+(\alpha_m-\gamma_m)x$$
elde edilir. Şimdi
$$\mu_t=\mathbb{E}[X_t],\qquad \nu_t=\mathbb{E}[X_t^2]$$
dersek, \(t\). adımın beklenen katkısı
$$\mathbb{E}[B_t^2]=\gamma_m\nu_t+(\alpha_m-\gamma_m)\mu_t$$
şeklindedir.
Durum geçişi \(X_{t+1}=X_t+k-B_t\) olduğundan, \(X_t=x\) koşulu altında
$$\mathbb{E}[X_{t+1}\mid X_t=x]=x+k-\mathbb{E}[B_t\mid X_t=x]=x+k-\frac{2}{m}x.$$
Dolayısıyla
$$\mathbb{E}[X_{t+1}\mid X_t=x]=\frac{m-2}{m}x+k.$$
Burada
$$\rho_m=\frac{m-2}{m}$$
tanımıyla
$$\mu_{t+1}=\rho_m\mu_t+k$$
bağıntısını buluruz. Başlangıç değeri \(\mu_1=k\)'dir.
Şimdi geçiş eşitliğinin karesini alalım:
$$X_{t+1}^2=(X_t+k-B_t)^2.$$
\(X_t=x\) koşulu altında beklenen değeri alıp Adım 2'deki formülleri yerleştirince
$$\mathbb{E}[X_{t+1}^2\mid X_t=x]=\mathbb{E}[(x+k-B_t)^2\mid X_t=x]$$
ifadesi, katsayılar toplanınca
$$\mathbb{E}[X_{t+1}^2\mid X_t=x]=\sigma_m x^2+\left(2k\rho_m+\rho_m-\sigma_m\right)x+k^2$$
biçimine iner; burada
$$\sigma_m=\frac{(m-2)(k(m-2)-1)}{m(km-1)}.$$
Böylece
$$\nu_{t+1}=\sigma_m\nu_t+\left(2k\rho_m+\rho_m-\sigma_m\right)\mu_t+k^2$$
elde edilir. Başlangıçta \(\nu_1=k^2\) olduğundan süreç bütünüyle \((\mu_t,\nu_t)\) çifti üzerinde kapanır.
İlk adımda \(m=3\) olur; yani popülasyon \(6\) bilye içerir ve \(X_1=2\)'dir. \(4\) bilye çekildiği için dağılım
$$\Pr(B_1=0)=\frac{\binom{2}{0}\binom{4}{4}}{\binom{6}{4}}=\frac{1}{15},$$
$$\Pr(B_1=1)=\frac{\binom{2}{1}\binom{4}{3}}{\binom{6}{4}}=\frac{8}{15},$$
$$\Pr(B_1=2)=\frac{\binom{2}{2}\binom{4}{2}}{\binom{6}{4}}=\frac{6}{15}$$
şeklindedir. Bu yüzden
$$\mathbb{E}[B_1^2]=0^2\cdot\frac{1}{15}+1^2\cdot\frac{8}{15}+2^2\cdot\frac{6}{15}=\frac{32}{15}.$$
Sonraki durumda \(X_2=4-B_1\) olduğundan \(X_2\), aynı olasılıklarla \(4,3,2\) değerlerini alır ve
$$\mathbb{E}[X_2^2]=4^2\cdot\frac{1}{15}+3^2\cdot\frac{8}{15}+2^2\cdot\frac{6}{15}=\frac{112}{15}$$
olur. İkinci adımda \(m=2\) olduğu için çekilen \(2k=4\) bilye, kalan popülasyonun tamamıdır. Dolayısıyla \(B_2=X_2\) ve
$$\mathbb{E}[B_2^2]=\mathbb{E}[X_2^2]=\frac{112}{15}.$$
Toplam beklenti
$$S(2,2)=\frac{32}{15}+\frac{112}{15}=\frac{144}{15}=\frac{48}{5}=9.6$$
çıkar; bu da uygulamadaki doğrulama değeriyle aynıdır.
C++, Python ve Java uygulamaları aynı moment bağıntılarını kullanır. Her biri yalnızca güncel birinci moment \(\mu_t\), güncel ikinci moment \(\nu_t\) ve biriken toplam \(\mathbb{E}[B_t^2]\) değerini tutar. Her adımda önce \(m=n-t+2\) hesaplanır, sonra
$$\gamma_m\nu_t+(\alpha_m-\gamma_m)\mu_t$$
katkısı toplama eklenir; eğer son adım değilse \(\mu_t\) ve \(\nu_t\) yukarıdaki doğrusal geçişlerle güncellenir.
Python ve Java sürümleri yüksek hassasiyetli ondalık aritmetik kullanır. C++ sürümü genişletilmiş kayan nokta aritmetiği kullanır ve ayrıca çok küçük örneklerde tam hipergeometrik dağılımı ilerletip moment bağıntısıyla karşılaştıran kendi iç doğrulamalarını da içerir. Döngü sonunda beklenen değer en yakın tam sayıya yuvarlanarak yazdırılır.
Ana algoritma her adımda sabit boyutlu bir güncelleme yaptığı için çalışma süresi \(O(n)\), bellek kullanımı \(O(1)\)'dir. Bu, olası izlenen bilye sayılarının tüm dağılımını taşımaya kıyasla temel kazançtır. Küçük doğrulama yöntemi tam dağılımı tuttuğu için çok daha pahalıdır; ancak yalnızca küçük testlerde kullanılır.
Las implementaciones usan una formulación equivalente del proceso de urnas basada en bolas rastreadas. En la etapa \(t\), definimos
$$m=n-t+2,\qquad P_t=km.$$
La población activa tiene \(P_t\) bolas, y \(X_t\) de ellas son las bolas cuyo efecto futuro todavía debe seguirse. Se extraen \(2k\) bolas de manera uniforme y sin reemplazo; \(B_t\) es el número de bolas rastreadas retiradas en esa extracción. La magnitud buscada es
$$S(n,k)=\mathbb{E}\left[\sum_{t=1}^{n} B_t^2\right].$$
Cuando \(n\) es grande, propagar toda la distribución de estados es demasiado costoso. La idea central es mantener solo el primer y el segundo momento de \(X_t\). El punto de control usado por la implementación es \(S(2,2)=9.6\).
La observación decisiva es que tanto la contribución de una etapa como la transición al siguiente estado dependen únicamente de los dos primeros momentos del número actual de bolas rastreadas.
Antes del sorteo en la etapa \(t\), hay \(X_t\) bolas rastreadas dentro de una población de tamaño \(P_t=km\), con \(m=n-t+2\). Condicionado a \(X_t=x\), el número de bolas rastreadas extraídas sigue una distribución hipergeométrica:
$$\Pr(B_t=b\mid X_t=x)=\frac{\binom{x}{b}\binom{km-x}{2k-b}}{\binom{km}{2k}}.$$
Después del sorteo, sobreviven las bolas rastreadas no extraídas y para la siguiente etapa entra un bloque nuevo de \(k\) bolas rastreadas. Por tanto, para \(t<n\),
$$X_{t+1}=X_t+k-B_t,$$
con valor inicial
$$X_1=k.$$
Esta identidad es la razón por la que se puede comprimir el proceso a unos pocos momentos.
Para una variable hipergeométrica, los primeros momentos factoriales son conocidos:
$$\mathbb{E}[B_t\mid X_t=x]=\frac{2k}{km}x=\frac{2}{m}x,$$
$$\mathbb{E}[B_t(B_t-1)\mid X_t=x]=\frac{2k(2k-1)}{km(km-1)}x(x-1)=\frac{2(2k-1)}{m(km-1)}x(x-1).$$
Definimos
$$\alpha_m=\frac{2}{m},\qquad \gamma_m=\frac{2(2k-1)}{m(km-1)}.$$
Como \(B_t^2=B_t(B_t-1)+B_t\), obtenemos
$$\mathbb{E}[B_t^2\mid X_t=x]=\gamma_m x^2+(\alpha_m-\gamma_m)x.$$
Si escribimos
$$\mu_t=\mathbb{E}[X_t],\qquad \nu_t=\mathbb{E}[X_t^2],$$
entonces la contribución esperada de la etapa \(t\) es
$$\mathbb{E}[B_t^2]=\gamma_m\nu_t+(\alpha_m-\gamma_m)\mu_t.$$
Partiendo de \(X_{t+1}=X_t+k-B_t\), y condicionando sobre \(X_t=x\), se tiene
$$\mathbb{E}[X_{t+1}\mid X_t=x]=x+k-\mathbb{E}[B_t\mid X_t=x]=x+k-\frac{2}{m}x.$$
Por tanto,
$$\mathbb{E}[X_{t+1}\mid X_t=x]=\frac{m-2}{m}x+k.$$
Definimos
$$\rho_m=\frac{m-2}{m}.$$
Al tomar esperanza, resulta
$$\mu_{t+1}=\rho_m\mu_t+k,$$
con condición inicial \(\mu_1=k\).
Ahora elevamos al cuadrado la transición de estado:
$$X_{t+1}^2=(X_t+k-B_t)^2.$$
Condicionando en \(X_t=x\) y sustituyendo las fórmulas del Paso 2, se obtiene
$$\mathbb{E}[X_{t+1}^2\mid X_t=x]=\mathbb{E}[(x+k-B_t)^2\mid X_t=x].$$
Tras reagrupar términos, aparece
$$\mathbb{E}[X_{t+1}^2\mid X_t=x]=\sigma_m x^2+\left(2k\rho_m+\rho_m-\sigma_m\right)x+k^2,$$
donde
$$\sigma_m=\frac{(m-2)(k(m-2)-1)}{m(km-1)}.$$
Por ello,
$$\nu_{t+1}=\sigma_m\nu_t+\left(2k\rho_m+\rho_m-\sigma_m\right)\mu_t+k^2,$$
con condición inicial \(\nu_1=k^2\). El proceso completo queda así cerrado sobre el par \((\mu_t,\nu_t)\).
En la primera etapa, \(m=3\), así que la población tiene \(6\) bolas y \(X_1=2\). Se extraen \(4\) bolas. La distribución es
$$\Pr(B_1=0)=\frac{\binom{2}{0}\binom{4}{4}}{\binom{6}{4}}=\frac{1}{15},$$
$$\Pr(B_1=1)=\frac{\binom{2}{1}\binom{4}{3}}{\binom{6}{4}}=\frac{8}{15},$$
$$\Pr(B_1=2)=\frac{\binom{2}{2}\binom{4}{2}}{\binom{6}{4}}=\frac{6}{15}.$$
Entonces
$$\mathbb{E}[B_1^2]=0^2\cdot\frac{1}{15}+1^2\cdot\frac{8}{15}+2^2\cdot\frac{6}{15}=\frac{32}{15}.$$
Además, \(X_2=4-B_1\), de modo que \(X_2\) toma los valores \(4,3,2\) con las mismas probabilidades, y
$$\mathbb{E}[X_2^2]=4^2\cdot\frac{1}{15}+3^2\cdot\frac{8}{15}+2^2\cdot\frac{6}{15}=\frac{112}{15}.$$
En la segunda etapa, \(m=2\), por lo que el tamaño de la extracción \(2k\) coincide con toda la población restante. Así \(B_2=X_2\), y
$$\mathbb{E}[B_2^2]=\mathbb{E}[X_2^2]=\frac{112}{15}.$$
Finalmente,
$$S(2,2)=\frac{32}{15}+\frac{112}{15}=\frac{144}{15}=\frac{48}{5}=9.6,$$
que coincide exactamente con el valor de comprobación de la implementación.
Las implementaciones en C++, Python y Java siguen la misma recurrencia. Solo guardan el primer momento actual \(\mu_t\), el segundo momento actual \(\nu_t\) y la suma acumulada de \(\mathbb{E}[B_t^2]\). En cada etapa calculan \(m=n-t+2\), evalúan la contribución
$$\gamma_m\nu_t+(\alpha_m-\gamma_m)\mu_t,$$
la añaden al total y, salvo en la última etapa, actualizan \(\mu_t\) y \(\nu_t\) con las fórmulas lineales anteriores.
Las versiones de Python y Java usan aritmética decimal de alta precisión, mientras que la versión de C++ usa coma flotante extendida. Además, la implementación en C++ incluye comprobaciones sobre casos pequeños donde se propaga la distribución hipergeométrica completa y se contrasta con la recurrencia de momentos. Al final, la esperanza calculada se redondea al entero más cercano antes de imprimirse.
El algoritmo principal realiza una actualización de tamaño constante por etapa, de modo que su coste es \(O(n)\) en tiempo y \(O(1)\) en memoria. Esa es la mejora esencial frente a una programación dinámica sobre la distribución completa, cuyo espacio de estados crece con los posibles valores del conteo rastreado. La rutina de validación pequeña es mucho más cara porque mantiene toda la distribución, pero solo se usa para pruebas diminutas.
这些实现把题目的 urn 过程改写成一个等价的“被追踪球”模型。在第 \(t\) 个阶段，记
$$m=n-t+2,\qquad P_t=km.$$
此时活动总体里共有 \(P_t\) 个球，其中有 \(X_t\) 个球仍然需要被继续追踪其后续影响。我们从这 \(P_t\) 个球中等概率地、不放回地抽出 \(2k\) 个球，并把其中被抽走的追踪球数量记为 \(B_t\)。目标量是
$$S(n,k)=\mathbb{E}\left[\sum_{t=1}^{n} B_t^2\right].$$
当 \(n\) 很大时，如果逐层维护完整的概率分布，代价会迅速失控。真正的关键在于：整个过程其实只需要维护 \(X_t\) 的一阶矩和二阶矩。实现里用到的一个样例校验值是 \(S(2,2)=9.6\)。
核心观察是：每一层对答案的贡献，以及下一层状态的更新，都只依赖当前追踪球数量的前两个矩，而不需要完整分布。
在第 \(t\) 层抽球之前，大小为 \(P_t=km\) 的总体中有 \(X_t\) 个追踪球，其中 \(m=n-t+2\)。在条件 \(X_t=x\) 下，抽到的追踪球数量 \(B_t\) 服从超几何分布：
$$\Pr(B_t=b\mid X_t=x)=\frac{\binom{x}{b}\binom{km-x}{2k-b}}{\binom{km}{2k}}.$$
抽取结束后，没有被抽走的追踪球会保留下来，并且下一层还会额外加入一个新的、大小为 \(k\) 的追踪球块。因此对 \(t<n\)，有精确关系
$$X_{t+1}=X_t+k-B_t,$$
而初始条件是
$$X_1=k.$$
这一条状态方程就是整道题的基础，因为一旦它成立，我们只要控制 \(X_t\) 的矩，就能同时控制当前层的贡献和下一层的状态。
对超几何随机变量，前两个阶乘矩有标准公式：
$$\mathbb{E}[B_t\mid X_t=x]=\frac{2k}{km}x=\frac{2}{m}x,$$
$$\mathbb{E}[B_t(B_t-1)\mid X_t=x]=\frac{2k(2k-1)}{km(km-1)}x(x-1)=\frac{2(2k-1)}{m(km-1)}x(x-1).$$
定义两个与阶段规模有关的系数
$$\alpha_m=\frac{2}{m},\qquad \gamma_m=\frac{2(2k-1)}{m(km-1)}.$$
由于
$$B_t^2=B_t(B_t-1)+B_t,$$
所以可得
$$\mathbb{E}[B_t^2\mid X_t=x]=\gamma_m x^2+(\alpha_m-\gamma_m)x.$$
再记
$$\mu_t=\mathbb{E}[X_t],\qquad \nu_t=\mathbb{E}[X_t^2],$$
那么第 \(t\) 层对总答案的期望贡献就是
$$\mathbb{E}[B_t^2]=\gamma_m\nu_t+(\alpha_m-\gamma_m)\mu_t.$$
也就是说，只要知道 \(\mu_t\) 和 \(\nu_t\)，这一层的贡献就能直接算出。
由状态方程
$$X_{t+1}=X_t+k-B_t$$
出发，在条件 \(X_t=x\) 下有
$$\mathbb{E}[X_{t+1}\mid X_t=x]=x+k-\mathbb{E}[B_t\mid X_t=x]=x+k-\frac{2}{m}x.$$
整理后得到
$$\mathbb{E}[X_{t+1}\mid X_t=x]=\frac{m-2}{m}x+k.$$
定义
$$\rho_m=\frac{m-2}{m}.$$
对上式再取期望，就得到线性递推
$$\mu_{t+1}=\rho_m\mu_t+k,$$
并且初始值为 \(\mu_1=k\)。这说明下一层的平均追踪球数量只依赖当前平均值。
接下来把状态方程平方：
$$X_{t+1}^2=(X_t+k-B_t)^2.$$
在条件 \(X_t=x\) 下取期望，再把步骤 2 里的条件矩代入，得到
$$\mathbb{E}[X_{t+1}^2\mid X_t=x]=\mathbb{E}[(x+k-B_t)^2\mid X_t=x].$$
把其中关于 \(x^2\)、\(x\) 和常数项的系数整理出来，可化为
$$\mathbb{E}[X_{t+1}^2\mid X_t=x]=\sigma_m x^2+\left(2k\rho_m+\rho_m-\sigma_m\right)x+k^2,$$
其中
$$\sigma_m=\frac{(m-2)(k(m-2)-1)}{m(km-1)}.$$
于是二阶矩满足
$$\nu_{t+1}=\sigma_m\nu_t+\left(2k\rho_m+\rho_m-\sigma_m\right)\mu_t+k^2,$$
初始值则是 \(\nu_1=k^2\)。至此可见，整个随机过程对 \((\mu_t,\nu_t)\) 这一对矩闭合，不再需要更大的状态描述。
第一层时 \(m=3\)，因此总体大小为 \(6\)，并且 \(X_1=2\)。此时抽取 \(4\) 个球。超几何分布写成
$$\Pr(B_1=0)=\frac{\binom{2}{0}\binom{4}{4}}{\binom{6}{4}}=\frac{1}{15},$$
$$\Pr(B_1=1)=\frac{\binom{2}{1}\binom{4}{3}}{\binom{6}{4}}=\frac{8}{15},$$
$$\Pr(B_1=2)=\frac{\binom{2}{2}\binom{4}{2}}{\binom{6}{4}}=\frac{6}{15}.$$
因此
$$\mathbb{E}[B_1^2]=0^2\cdot\frac{1}{15}+1^2\cdot\frac{8}{15}+2^2\cdot\frac{6}{15}=\frac{32}{15}.$$
下一层状态满足 \(X_2=4-B_1\)，所以 \(X_2\) 分别以相同概率取 \(4,3,2\)，从而
$$\mathbb{E}[X_2^2]=4^2\cdot\frac{1}{15}+3^2\cdot\frac{8}{15}+2^2\cdot\frac{6}{15}=\frac{112}{15}.$$
第二层时 \(m=2\)，于是抽取规模 \(2k=4\) 正好等于剩余总体的全部大小，所以
$$B_2=X_2,$$
从而
$$\mathbb{E}[B_2^2]=\mathbb{E}[X_2^2]=\frac{112}{15}.$$
最后得到
$$S(2,2)=\frac{32}{15}+\frac{112}{15}=\frac{144}{15}=\frac{48}{5}=9.6,$$
与实现中的校验值完全一致。
C++、Python 和 Java 三个实现都遵循同一套递推。它们只保存当前的一阶矩 \(\mu_t\)、当前的二阶矩 \(\nu_t\)，以及累计答案 \(\sum \mathbb{E}[B_t^2]\)。每一层先计算 \(m=n-t+2\)，然后用
$$\gamma_m\nu_t+(\alpha_m-\gamma_m)\mu_t$$
更新该层贡献；如果还没有到最后一层，再用上面的线性公式更新 \(\mu_t\) 和 \(\nu_t\)。
Python 与 Java 版本使用高精度十进制运算；C++ 版本使用扩展浮点精度。此外，C++ 版本还包含了很小规模的自检：它直接传播完整的超几何分布，并与矩递推的结果逐项比较。循环结束后，最终期望值会被四舍五入到最接近的整数再输出。
正式求解算法在每一层只做常数次运算，因此总时间复杂度为 \(O(n)\)，空间复杂度为 \(O(1)\)。这正是它优于完整分布动态规划的地方，因为完整分布的方法需要维护所有可能的追踪球计数。小规模验证例程由于要保存整条分布，代价高得多，但只用于极小输入的自检。
Реализации используют эквивалентную модель процесса с урнами через число отслеживаемых шаров. На шаге \(t\) положим
$$m=n-t+2,\qquad P_t=km.$$
В активной совокупности находится \(P_t\) шаров, и \(X_t\) из них нужно продолжать отслеживать в дальнейших шагах. Из этой совокупности равновероятно и без возвращения выбираются \(2k\) шаров; через \(B_t\) обозначается число удаленных отслеживаемых шаров. Нужно вычислить
$$S(n,k)=\mathbb{E}\left[\sum_{t=1}^{n} B_t^2\right].$$
Если хранить полное распределение по состояниям, то при больших \(n\) задача становится слишком дорогой. Поэтому решение сводит всё к эволюции только первого и второго моментов величины \(X_t\). Контрольное значение из реализации: \(S(2,2)=9.6\).
Главное наблюдение состоит в том, что и вклад текущего шага в ответ, и переход к следующему шагу зависят только от первых двух моментов текущего числа отслеживаемых шаров.
Перед выборкой на шаге \(t\) внутри популяции размера \(P_t=km\), где \(m=n-t+2\), находится \(X_t\) отслеживаемых шаров. При условии \(X_t=x\) число выбранных отслеживаемых шаров имеет гипергеометрическое распределение:
$$\Pr(B_t=b\mid X_t=x)=\frac{\binom{x}{b}\binom{km-x}{2k-b}}{\binom{km}{2k}}.$$
После выборки невыбранные отслеживаемые шары сохраняются, а на следующий шаг добавляется новый блок из \(k\) отслеживаемых шаров. Поэтому для \(t<n\)
$$X_{t+1}=X_t+k-B_t,$$
а начальное значение равно
$$X_1=k.$$
Именно эта формула позволяет заменить полное распределение коротким набором моментов.
Для гипергеометрической случайной величины стандартны следующие факториальные моменты:
$$\mathbb{E}[B_t\mid X_t=x]=\frac{2k}{km}x=\frac{2}{m}x,$$
$$\mathbb{E}[B_t(B_t-1)\mid X_t=x]=\frac{2k(2k-1)}{km(km-1)}x(x-1)=\frac{2(2k-1)}{m(km-1)}x(x-1).$$
Введем обозначения
$$\alpha_m=\frac{2}{m},\qquad \gamma_m=\frac{2(2k-1)}{m(km-1)}.$$
Так как
$$B_t^2=B_t(B_t-1)+B_t,$$
получаем
$$\mathbb{E}[B_t^2\mid X_t=x]=\gamma_m x^2+(\alpha_m-\gamma_m)x.$$
Если обозначить
$$\mu_t=\mathbb{E}[X_t],\qquad \nu_t=\mathbb{E}[X_t^2],$$
то вклад шага \(t\) в среднем равен
$$\mathbb{E}[B_t^2]=\gamma_m\nu_t+(\alpha_m-\gamma_m)\mu_t.$$
Из равенства \(X_{t+1}=X_t+k-B_t\) при условии \(X_t=x\) следует
$$\mathbb{E}[X_{t+1}\mid X_t=x]=x+k-\mathbb{E}[B_t\mid X_t=x]=x+k-\frac{2}{m}x.$$
Следовательно,
$$\mathbb{E}[X_{t+1}\mid X_t=x]=\frac{m-2}{m}x+k.$$
Положим
$$\rho_m=\frac{m-2}{m}.$$
После взятия полного ожидания получаем линейную рекурсию
$$\mu_{t+1}=\rho_m\mu_t+k,$$
с начальным значением \(\mu_1=k\).
Теперь возведем переход в квадрат:
$$X_{t+1}^2=(X_t+k-B_t)^2.$$
Условное ожидание при \(X_t=x\) и подстановка формул из шага 2 дают
$$\mathbb{E}[X_{t+1}^2\mid X_t=x]=\mathbb{E}[(x+k-B_t)^2\mid X_t=x].$$
После группировки коэффициентов получается
$$\mathbb{E}[X_{t+1}^2\mid X_t=x]=\sigma_m x^2+\left(2k\rho_m+\rho_m-\sigma_m\right)x+k^2,$$
где
$$\sigma_m=\frac{(m-2)(k(m-2)-1)}{m(km-1)}.$$
Отсюда следует рекурсия
$$\nu_{t+1}=\sigma_m\nu_t+\left(2k\rho_m+\rho_m-\sigma_m\right)\mu_t+k^2,$$
с начальным значением \(\nu_1=k^2\). Значит, процесс полностью замыкается на паре \((\mu_t,\nu_t)\).
На первом шаге \(m=3\), так что всего имеется \(6\) шаров и \(X_1=2\). Выбираются \(4\) шара. Тогда
$$\Pr(B_1=0)=\frac{\binom{2}{0}\binom{4}{4}}{\binom{6}{4}}=\frac{1}{15},$$
$$\Pr(B_1=1)=\frac{\binom{2}{1}\binom{4}{3}}{\binom{6}{4}}=\frac{8}{15},$$
$$\Pr(B_1=2)=\frac{\binom{2}{2}\binom{4}{2}}{\binom{6}{4}}=\frac{6}{15}.$$
Следовательно,
$$\mathbb{E}[B_1^2]=0^2\cdot\frac{1}{15}+1^2\cdot\frac{8}{15}+2^2\cdot\frac{6}{15}=\frac{32}{15}.$$
Далее \(X_2=4-B_1\), поэтому \(X_2\) принимает значения \(4,3,2\) с теми же вероятностями, а значит
$$\mathbb{E}[X_2^2]=4^2\cdot\frac{1}{15}+3^2\cdot\frac{8}{15}+2^2\cdot\frac{6}{15}=\frac{112}{15}.$$
На втором шаге \(m=2\), и размер выборки \(2k=4\) совпадает со всей оставшейся популяцией. Поэтому \(B_2=X_2\), и
$$\mathbb{E}[B_2^2]=\mathbb{E}[X_2^2]=\frac{112}{15}.$$
Итак,
$$S(2,2)=\frac{32}{15}+\frac{112}{15}=\frac{144}{15}=\frac{48}{5}=9.6,$$
что полностью совпадает с проверочным значением в реализации.
Реализации на C++, Python и Java используют одну и ту же схему. Они хранят только текущий первый момент \(\mu_t\), текущий второй момент \(\nu_t\) и накапливаемую сумму \(\mathbb{E}[B_t^2]\). На каждом шаге вычисляется \(m=n-t+2\), затем вклад
$$\gamma_m\nu_t+(\alpha_m-\gamma_m)\mu_t$$
добавляется к ответу, а если шаг не последний, то \(\mu_t\) и \(\nu_t\) обновляются по линейным формулам выше.
Версии на Python и Java используют высокоточную десятичную арифметику, а версия на C++ работает с расширенной плавающей точкой. Кроме того, реализация на C++ содержит самопроверки на маленьких случаях: там напрямую распространяется полное гипергеометрическое распределение и сравнивается с ответом, полученным через моменты. После завершения цикла итоговое ожидание округляется до ближайшего целого и выводится.
Основной алгоритм выполняет на каждом шаге лишь константное число операций, поэтому работает за \(O(n)\) времени и использует \(O(1)\) памяти. Именно это делает его практичным по сравнению с полным динамическим учетом распределения, где число состояний растет вместе с возможными значениями отслеживаемого счетчика. Маленькая проверочная процедура значительно дороже, потому что хранит всю распределенную массу, но используется только для крошечных тестов.
تستخدم التطبيقات صياغة مكافئة لعملية الأوعية تعتمد على عدد الكرات التي ما زال يلزم تتبّعها. في المرحلة \(t\) نعرّف
$$m=n-t+2,\qquad P_t=km.$$
يكون لدينا مجتمع فعّال حجمه \(P_t\) كرة، ومن بينها \(X_t\) كرة ما زال أثرها المستقبلي مهمًا في الحساب. نسحب \(2k\) كرة عشوائيًا ومن دون إرجاع، ولتكن \(B_t\) هي عدد الكرات المتتبَّعة التي أزيلت في هذا السحب. الكمية المطلوبة هي
$$S(n,k)=\mathbb{E}\left[\sum_{t=1}^{n} B_t^2\right].$$
عندما يكون \(n\) كبيرًا يصبح الاحتفاظ بالتوزيع الاحتمالي الكامل مكلفًا جدًا. لذلك تعتمد الفكرة الأساسية على متابعة العزم الأول والعزم الثاني فقط للمتغير \(X_t\). وقيمة التحقق الصغيرة المستخدمة في التطبيق هي \(S(2,2)=9.6\).
الملاحظة الحاسمة هي أن مساهمة كل مرحلة في الجواب، وكذلك انتقال الحالة إلى المرحلة التالية، يعتمدان فقط على أول عزميْن لعدد الكرات المتتبَّعة في تلك المرحلة.
قبل السحب في المرحلة \(t\)، يوجد \(X_t\) كرة متتبَّعة داخل مجتمع حجمه \(P_t=km\)، حيث \(m=n-t+2\). وبشرط \(X_t=x\)، فإن عدد الكرات المتتبَّعة المسحوبة يتبع توزيعًا فوق هندسي:
$$\Pr(B_t=b\mid X_t=x)=\frac{\binom{x}{b}\binom{km-x}{2k-b}}{\binom{km}{2k}}.$$
بعد السحب تبقى الكرات المتتبَّعة التي لم تُسحب، ثم تدخل كتلة جديدة حجمها \(k\) من الكرات المتتبَّعة إلى المرحلة التالية. لذلك عندما \(t<n\) يكون
$$X_{t+1}=X_t+k-B_t,$$
ومع الشرط الابتدائي
$$X_1=k.$$
هذه العلاقة هي أساس الحل كله، لأنها تسمح باستبدال التوزيع الكامل بتتبع عدد صغير من العزوم.
بالنسبة إلى المتغير الفوق هندسي، فإن العزوم العاملية الأولى معروفة:
$$\mathbb{E}[B_t\mid X_t=x]=\frac{2k}{km}x=\frac{2}{m}x,$$
$$\mathbb{E}[B_t(B_t-1)\mid X_t=x]=\frac{2k(2k-1)}{km(km-1)}x(x-1)=\frac{2(2k-1)}{m(km-1)}x(x-1).$$
نعرّف
$$\alpha_m=\frac{2}{m},\qquad \gamma_m=\frac{2(2k-1)}{m(km-1)}.$$
وبما أن
$$B_t^2=B_t(B_t-1)+B_t,$$
فإن
$$\mathbb{E}[B_t^2\mid X_t=x]=\gamma_m x^2+(\alpha_m-\gamma_m)x.$$
إذا كتبنا
$$\mu_t=\mathbb{E}[X_t],\qquad \nu_t=\mathbb{E}[X_t^2],$$
فإن المساهمة المتوقعة للمرحلة \(t\) تصبح
$$\mathbb{E}[B_t^2]=\gamma_m\nu_t+(\alpha_m-\gamma_m)\mu_t.$$
من العلاقة \(X_{t+1}=X_t+k-B_t\) نحصل، بشرط \(X_t=x\)، على
$$\mathbb{E}[X_{t+1}\mid X_t=x]=x+k-\mathbb{E}[B_t\mid X_t=x]=x+k-\frac{2}{m}x.$$
أي
$$\mathbb{E}[X_{t+1}\mid X_t=x]=\frac{m-2}{m}x+k.$$
لنعرّف
$$\rho_m=\frac{m-2}{m}.$$
وبأخذ التوقع الكامل نحصل على العلاقة الخطية
$$\mu_{t+1}=\rho_m\mu_t+k,$$
مع الشرط الابتدائي \(\mu_1=k\).
نربّع الآن انتقال الحالة:
$$X_{t+1}^2=(X_t+k-B_t)^2.$$
وبأخذ التوقع الشرطي عند \(X_t=x\) ثم إدخال صيغ الخطوة 2، نحصل على
$$\mathbb{E}[X_{t+1}^2\mid X_t=x]=\mathbb{E}[(x+k-B_t)^2\mid X_t=x].$$
وبعد جمع معاملات \(x^2\) و\(x\) والحد الثابت يظهر الشكل
$$\mathbb{E}[X_{t+1}^2\mid X_t=x]=\sigma_m x^2+\left(2k\rho_m+\rho_m-\sigma_m\right)x+k^2,$$
حيث
$$\sigma_m=\frac{(m-2)(k(m-2)-1)}{m(km-1)}.$$
ومن ثم
$$\nu_{t+1}=\sigma_m\nu_t+\left(2k\rho_m+\rho_m-\sigma_m\right)\mu_t+k^2,$$
مع الشرط الابتدائي \(\nu_1=k^2\). وهكذا ينغلق المسار الاحتمالي كله على الزوج \((\mu_t,\nu_t)\) فقط.
في المرحلة الأولى يكون \(m=3\)، وبالتالي حجم المجتمع \(6\) و\(X_1=2\). نسحب \(4\) كرات، وعندها
$$\Pr(B_1=0)=\frac{\binom{2}{0}\binom{4}{4}}{\binom{6}{4}}=\frac{1}{15},$$
$$\Pr(B_1=1)=\frac{\binom{2}{1}\binom{4}{3}}{\binom{6}{4}}=\frac{8}{15},$$
$$\Pr(B_1=2)=\frac{\binom{2}{2}\binom{4}{2}}{\binom{6}{4}}=\frac{6}{15}.$$
إذن
$$\mathbb{E}[B_1^2]=0^2\cdot\frac{1}{15}+1^2\cdot\frac{8}{15}+2^2\cdot\frac{6}{15}=\frac{32}{15}.$$
بعد ذلك يصبح
$$X_2=4-B_1,$$
ولذلك يأخذ \(X_2\) القيم \(4,3,2\) بالاحتمالات نفسها، ومنه
$$\mathbb{E}[X_2^2]=4^2\cdot\frac{1}{15}+3^2\cdot\frac{8}{15}+2^2\cdot\frac{6}{15}=\frac{112}{15}.$$
في المرحلة الثانية يكون \(m=2\)، فيتساوى حجم السحب \(2k=4\) مع كل المجتمع المتبقي. لذا
$$B_2=X_2,$$
وبالتالي
$$\mathbb{E}[B_2^2]=\mathbb{E}[X_2^2]=\frac{112}{15}.$$
وأخيرًا
$$S(2,2)=\frac{32}{15}+\frac{112}{15}=\frac{144}{15}=\frac{48}{5}=9.6,$$
وهو تمامًا مقدار التحقق الموجود في التنفيذ.
تتبع تطبيقات C++ وPython وJava العلاقة نفسها. فهي تحتفظ فقط بالعزم الأول الحالي \(\mu_t\)، والعزم الثاني الحالي \(\nu_t\)، والمجموع التراكمي للحدود \(\mathbb{E}[B_t^2]\). في كل مرحلة يُحسب \(m=n-t+2\)، ثم تضاف المساهمة
$$\gamma_m\nu_t+(\alpha_m-\gamma_m)\mu_t$$
إلى المجموع، وبعدها تُحدَّث \(\mu_t\) و\(\nu_t\) بالعلاقات الخطية السابقة ما لم تكن هذه هي المرحلة الأخيرة.
إصدارا Python وJava يستخدمان حسابًا عشريًا عالي الدقة، بينما يستخدم إصدار C++ دقة عائمة موسعة. كما يحتوي إصدار C++ على اختبارات صغيرة ينشر فيها التوزيع الفوق هندسي الكامل للحالات الصغيرة ويقارنه بنتيجة طريقة العزوم. بعد انتهاء الحلقة يُقرَّب التوقع النهائي إلى أقرب عدد صحيح ثم يُطبع.
الخوارزمية الأساسية تنفذ تحديثًا ثابت الحجم في كل مرحلة، ولذلك تعمل في زمن \(O(n)\) وتستخدم ذاكرة \(O(1)\). وهذا هو التفوق الجوهري على أسلوب يتابع التوزيع الكامل، لأن فضاء الحالات هناك ينمو مع عدد القيم الممكنة لعدد الكرات المتتبَّعة. أما روتين التحقق الصغير فهو أعلى تكلفة بكثير لأنه يحتفظ بالتوزيع كله، لكنه يُستخدم فقط على أمثلة صغيرة جدًا.