Place the \(D\) days on a cycle. One random draw chooses a day uniformly, and that draw covers two consecutive days: the chosen day and its cyclic successor. Let \(T_D\) be the number of draws needed until every day has been covered at least once. The task is to compute the expectation
$$E(D)=\mathbb{E}[T_D]$$
for the large target value \(D=10000\).
The C++, Python, and Java implementations do not simulate the process. Instead, they turn the expectation into a one-dimensional integral whose integrand comes from inclusion-exclusion on uncovered days and a \(2\times2\) transfer matrix on the cycle.
Number the days modulo \(D\), and let a draw on day \(i\) cover the pair \((i,i+1)\). The stopping time is \(T_D\), the first moment when all days are covered.
For any nonnegative integer-valued stopping time,
$$E(D)=\sum_{n=0}^{\infty}\Pr(T_D>n).$$
So the whole problem becomes: after \(n\) draws, what is the probability that at least one day is still uncovered?
A day \(j\) is uncovered exactly when no draw landed on \(j\) and no draw landed on \(j-1\), because those are the only two starting points that could cover \(j\).
For a subset \(S\) of days, define
$$N(S)=S\cup(S-1), \qquad S-1=\{i-1 \pmod D : i\in S\}.$$
If every day in \(S\) is still uncovered after \(n\) draws, then every draw must avoid the set \(N(S)\). Since each draw is uniform over the \(D\) starting days,
$$\Pr(A_S(n))=\left(1-\frac{|N(S)|}{D}\right)^n,$$
where \(A_S(n)\) is the event that all days in \(S\) remain uncovered after \(n\) draws.
Now inclusion-exclusion over the uncovered days gives
$$\Pr(T_D>n)=\sum_{\emptyset\ne S\subseteq C_D}(-1)^{|S|+1}\left(1-\frac{|N(S)|}{D}\right)^n,$$
where \(C_D\) denotes the cycle of \(D\) days.
Insert the previous formula into the tail sum and interchange the order of summation:
$$\frac{E(D)}{D}=\sum_{\emptyset\ne S\subseteq C_D}(-1)^{|S|+1}\frac{1}{|N(S)|}.$$
The reciprocal is converted with
$$\frac{1}{m}=\int_0^1 x^{m-1}\,dx.$$
Therefore
$$\frac{E(D)}{D}=\int_0^1\frac{1-Z_D(x)}{x}\,dx,$$
where
$$Z_D(x)=\sum_{S\subseteq C_D}(-1)^{|S|}x^{|N(S)|}.$$
The empty set contributes the leading \(1\), which is why the numerator becomes \(1-Z_D(x)\).
Represent a subset \(S\) by a cyclic binary word \(s_1,\dots,s_D\), where \(s_i=1\) means day \(i\) is required to remain uncovered. Then
$$|N(S)|=\sum_{i=1}^{D}s_i+\sum_{i=1}^{D}(1-s_i)s_{i+1}, \qquad s_{D+1}=s_1.$$
The first sum counts the chosen days themselves, and the second sum counts one extra predecessor for each run of consecutive \(1\)'s around the cycle.
This makes the weight factorize locally:
$$(-1)^{|S|}x^{|N(S)|}=\prod_{i=1}^{D} M_{s_i,s_{i+1}}(x),$$
with transfer matrix
$$M(x)=\begin{pmatrix}1 & -x^2 \\ 1 & -x\end{pmatrix}.$$
Summing over all cyclic binary words is exactly the trace of the \(D\)-th matrix power, so
$$Z_D(x)=\operatorname{tr}(M(x)^D).$$
The characteristic polynomial of \(M(x)\) is
$$\lambda^2-(1-x)\lambda-x(1-x)=0.$$
Its two eigenvalues are
$$\lambda_{1,2}(x)=\frac{1-x\pm\sqrt{(1-x)(1+3x)}}{2}.$$
Hence
$$Z_D(x)=\lambda_1(x)^D+\lambda_2(x)^D,$$
and the expectation becomes
$$\boxed{\frac{E(D)}{D}=\int_0^1\frac{1-\lambda_1(x)^D-\lambda_2(x)^D}{x}\,dx.}$$
This is the exact formula used by the implementation.
The integrand looks singular because of the division by \(x\), but the singularity at \(x=0\) is removable. Expanding the eigenvalues near zero gives
$$\lambda_1(x)=1-x^2+O(x^3), \qquad \lambda_2(x)=-x+O(x^2),$$
so
$$1-\lambda_1(x)^D-\lambda_2(x)^D = D x^2 + O(x^3),$$
and therefore the integrand is \(Dx+O(x^2)\), which tends to \(0\). At \(x=1\), both eigenvalues are \(0\), so the integrand tends to \(1\).
For five days, the transfer-matrix trace simplifies to
$$Z_5(x)=1-5x^2+5x^3-x^5.$$
Thus
$$\frac{E(5)}{5}=\int_0^1\frac{5x^2-5x^3+x^5}{x}\,dx=\int_0^1(5x-5x^2+x^4)\,dx.$$
Integrating term by term gives
$$\frac{E(5)}{5}=\frac{5}{2}-\frac{5}{3}+\frac{1}{5}=\frac{31}{30},$$
so
$$E(5)=\frac{31}{6}.$$
This matches the exact checkpoint used by the implementations. For the smallest nontrivial cycle, \(Z_2(x)=1-x^2\), giving \(E(2)=1\).
The C++, Python, and Java implementations all evaluate the same integrand
$$I_D(x)=\frac{1-\lambda_1(x)^D-\lambda_2(x)^D}{x}.$$
They first compute the discriminant \(\sqrt{(1-x)(1+3x)}\), then the two eigenvalues, and then the numerator. The endpoints are handled explicitly: the limiting value at \(x=0\) is \(0\), and the value at \(x=1\) is \(1\).
The delicate part is the term \(1-\lambda_1(x)^D\). When \(\lambda_1(x)\) is very close to \(1\), direct subtraction loses significant digits. To avoid that, the implementation rewrites it as
$$1-\lambda_1(x)^D=-\operatorname{expm1}\!\bigl(D\log \lambda_1(x)\bigr)$$
whenever \(D\log \lambda_1(x)\) is close to \(0\). The second eigenvalue is smaller in magnitude and is raised directly to the \(D\)-th power.
The integral on \([0,1]\) is then evaluated with Simpson's rule using an even number of intervals. Starting from an initial mesh, the code repeatedly doubles the interval count and compares two consecutive Simpson estimates:
$$|S_{2N}-S_N|\le \text{rtol}\cdot \max(1,|S_{2N}|).$$
Once this relative criterion is satisfied, the refined estimate is returned. The C++ implementation also supports parallel accumulation of the interior Simpson nodes when the mesh is large, and it verifies the derivation against the known values \(E(2)=1\), \(E(5)=31/6\), and \(E(365)\approx 1174.3501\).
If the final Simpson mesh uses \(N\) subintervals, one pass costs \(O(N)\) evaluations of the integrand. Because the mesh doubles geometrically, the total work up to the accepted mesh is still \(O(N)\). The single-threaded implementations use \(O(1)\) auxiliary memory, while the parallel C++ version uses \(O(T)\) extra memory for \(T\) partial sums.
Ordne die \(D\) Tage auf einem Kreis an. Ein zufälliger Zug wählt einen Tag gleichverteilt aus und deckt genau zwei aufeinanderfolgende Tage ab: den gewählten Tag und seinen zyklischen Nachfolger. Sei \(T_D\) die Anzahl der Züge bis jeder Tag mindestens einmal abgedeckt wurde. Gesucht ist der Erwartungswert
$$E(D)=\mathbb{E}[T_D]$$
für den großen Zielwert \(D=10000\).
Die C++-, Python- und Java-Implementierungen simulieren diesen Prozess nicht. Stattdessen formen sie den Erwartungswert in ein eindimensionales Integral um, dessen Integrand aus Inklusion-Exklusion über unbedeckte Tage und einer \(2\times2\)-Transfermatrix auf dem Kreis entsteht.
Nummeriere die Tage modulo \(D\). Ein Zug auf Tag \(i\) deckt das Paar \((i,i+1)\) ab. Die Stoppzeit \(T_D\) ist der erste Zeitpunkt, an dem alle Tage abgedeckt sind.
Für jede nichtnegative ganzzahlige Stoppzeit gilt
$$E(D)=\sum_{n=0}^{\infty}\Pr(T_D>n).$$
Damit reduziert sich alles auf die Frage: Wie groß ist nach \(n\) Zügen die Wahrscheinlichkeit, dass noch mindestens ein Tag unbedeckt ist?
Ein Tag \(j\) ist genau dann unbedeckt, wenn kein Zug auf \(j\) und kein Zug auf \(j-1\) gelandet ist, denn nur diese beiden Startpunkte könnten \(j\) abdecken.
Für eine Teilmenge \(S\) von Tagen definieren wir
$$N(S)=S\cup(S-1), \qquad S-1=\{i-1 \pmod D : i\in S\}.$$
Wenn nach \(n\) Zügen alle Tage aus \(S\) noch unbedeckt sind, dann muss jeder Zug die Menge \(N(S)\) meiden. Da jeder Zug gleichverteilt auf den \(D\) möglichen Starttagen liegt, folgt
$$\Pr(A_S(n))=\left(1-\frac{|N(S)|}{D}\right)^n,$$
wobei \(A_S(n)\) das Ereignis bezeichnet, dass alle Tage aus \(S\) nach \(n\) Zügen unbedeckt bleiben.
Mit Inklusion-Exklusion über die unbedeckten Tage erhalten wir
$$\Pr(T_D>n)=\sum_{\emptyset\ne S\subseteq C_D}(-1)^{|S|+1}\left(1-\frac{|N(S)|}{D}\right)^n,$$
wobei \(C_D\) den Kreis aus \(D\) Tagen bezeichnet.
Setzt man diese Formel in die Schwanzsumme ein und vertauscht die Summationen, so ergibt sich
$$\frac{E(D)}{D}=\sum_{\emptyset\ne S\subseteq C_D}(-1)^{|S|+1}\frac{1}{|N(S)|}.$$
Den Kehrwert schreiben wir als
$$\frac{1}{m}=\int_0^1 x^{m-1}\,dx.$$
Daraus folgt
$$\frac{E(D)}{D}=\int_0^1\frac{1-Z_D(x)}{x}\,dx,$$
mit
$$Z_D(x)=\sum_{S\subseteq C_D}(-1)^{|S|}x^{|N(S)|}.$$
Die leere Menge liefert den führenden Term \(1\); deshalb steht im Zähler \(1-Z_D(x)\).
Stelle eine Teilmenge \(S\) als zyklisches Binärwort \(s_1,\dots,s_D\) dar, wobei \(s_i=1\) bedeutet, dass Tag \(i\) unbedeckt bleiben soll. Dann gilt
$$|N(S)|=\sum_{i=1}^{D}s_i+\sum_{i=1}^{D}(1-s_i)s_{i+1}, \qquad s_{D+1}=s_1.$$
Die erste Summe zählt die markierten Tage selbst, die zweite Summe genau einen zusätzlichen Vorgänger pro Block aufeinanderfolgender Einsen.
Dadurch faktorisiert das Gewicht lokal:
$$(-1)^{|S|}x^{|N(S)|}=\prod_{i=1}^{D} M_{s_i,s_{i+1}}(x),$$
mit der Transfermatrix
$$M(x)=\begin{pmatrix}1 & -x^2 \\ 1 & -x\end{pmatrix}.$$
Die Summe über alle zyklischen Binärwörter ist genau die Spur der \(D\)-ten Matrixpotenz, also
$$Z_D(x)=\operatorname{tr}(M(x)^D).$$
Das charakteristische Polynom von \(M(x)\) lautet
$$\lambda^2-(1-x)\lambda-x(1-x)=0.$$
Die beiden Eigenwerte sind
$$\lambda_{1,2}(x)=\frac{1-x\pm\sqrt{(1-x)(1+3x)}}{2}.$$
Damit erhält man
$$Z_D(x)=\lambda_1(x)^D+\lambda_2(x)^D,$$
und somit
$$\boxed{\frac{E(D)}{D}=\int_0^1\frac{1-\lambda_1(x)^D-\lambda_2(x)^D}{x}\,dx.}$$
Genau diese geschlossene Formel wertet die Implementierung numerisch aus.
Die Division durch \(x\) wirkt zunächst problematisch, aber bei \(x=0\) liegt nur eine hebbare Singularität vor. Für kleine \(x\) gilt
$$\lambda_1(x)=1-x^2+O(x^3), \qquad \lambda_2(x)=-x+O(x^2),$$
also
$$1-\lambda_1(x)^D-\lambda_2(x)^D = D x^2 + O(x^3).$$
Der Integrand verhält sich daher wie \(Dx+O(x^2)\) und geht gegen \(0\). Bei \(x=1\) verschwinden beide Eigenwerte, also ist der Grenzwert dort gleich \(1\).
Für fünf Tage vereinfacht sich die Spur der Transfermatrix zu
$$Z_5(x)=1-5x^2+5x^3-x^5.$$
Daher
$$\frac{E(5)}{5}=\int_0^1\frac{5x^2-5x^3+x^5}{x}\,dx=\int_0^1(5x-5x^2+x^4)\,dx.$$
Termweise Integration liefert
$$\frac{E(5)}{5}=\frac{5}{2}-\frac{5}{3}+\frac{1}{5}=\frac{31}{30},$$
also
$$E(5)=\frac{31}{6}.$$
Das stimmt exakt mit dem Kontrollwert der Implementierungen überein. Für \(D=2\) ist \(Z_2(x)=1-x^2\), also \(E(2)=1\).
Die C++-, Python- und Java-Implementierungen werten alle denselben Integranden aus:
$$I_D(x)=\frac{1-\lambda_1(x)^D-\lambda_2(x)^D}{x}.$$
Dazu berechnen sie zuerst die Diskriminante \(\sqrt{(1-x)(1+3x)}\), dann die beiden Eigenwerte und daraus den Zähler. Die Randpunkte werden getrennt behandelt: der Grenzwert bei \(x=0\) ist \(0\), bei \(x=1\) ist der Wert \(1\).
Der numerisch heikle Teil ist \(1-\lambda_1(x)^D\). Ist \(\lambda_1(x)\) sehr nahe bei \(1\), gehen bei direkter Subtraktion signifikante Stellen verloren. Deshalb verwendet die Implementierung dann die stabile Umformung
$$1-\lambda_1(x)^D=-\operatorname{expm1}\!\bigl(D\log \lambda_1(x)\bigr).$$
Der zweite Eigenwert ist betragsmäßig kleiner und wird direkt potenziert.
Das Integral auf \([0,1]\) wird anschließend mit der Simpson-Regel für eine gerade Zahl von Teilintervallen berechnet. Beginnend mit einem Startgitter verdoppelt der Code wiederholt die Anzahl der Intervalle und vergleicht zwei aufeinanderfolgende Simpson-Näherungen:
$$|S_{2N}-S_N|\le \text{rtol}\cdot \max(1,|S_{2N}|).$$
Sobald dieses relative Kriterium erfüllt ist, wird der verfeinerte Wert zurückgegeben. Die C++-Implementierung kann zusätzlich die inneren Simpson-Stützstellen parallel aufsummieren und prüft die Herleitung an den bekannten Werten \(E(2)=1\), \(E(5)=31/6\) und \(E(365)\approx 1174.3501\).
Verwendet das endgültige Simpson-Gitter \(N\) Teilintervalle, dann kostet ein Durchlauf \(O(N)\) Integrand-Auswertungen. Da die Anzahl der Intervalle geometrisch verdoppelt wird, bleibt auch die gesamte Arbeit bis zum akzeptierten Gitter \(O(N)\). Die seriellen Implementierungen benötigen \(O(1)\) Zusatzspeicher, die parallele C++-Variante \(O(T)\) für \(T\) Teilsummen.
\(D\) günü bir çember üzerinde düşünelim. Her rastgele seçim, günlerden birini eşit olasılıkla seçer ve art arda iki günü kapsar: seçilen gün ile onun döngüsel ardılı. \(T_D\), bütün günler en az bir kez kapsanana kadar gereken seçim sayısı olsun. Aranan büyüklük
$$E(D)=\mathbb{E}[T_D]$$
olup hedef büyük değer \(D=10000\)'dir.
C++, Python ve Java implementasyonları bu süreci benzetimle çözmez. Onun yerine beklentiyi, açıkta kalan günler üzerinde dahil et-çıkar ve çember üzerindeki \(2\times2\) bir transfer matrisiyle elde edilen tek boyutlu bir integrale dönüştürürler.
Günleri modulo \(D\) numaralandıralım. \(i\) gününde yapılan bir seçim \((i,i+1)\) çiftini kapsar. \(T_D\), tüm günlerin ilk kez kapsandığı andır.
Negatif olmayan tam sayı değerli her durma zamanı için
$$E(D)=\sum_{n=0}^{\infty}\Pr(T_D>n)$$
kimliği geçerlidir.
Dolayısıyla esas soru şudur: \(n\) seçimden sonra en az bir günün hâlâ kapsanmamış olma olasılığı nedir?
\(j\) günü, ancak hiçbir seçim \(j\)'ye ve hiçbir seçim \(j-1\)'e gelmediyse açıkta kalır; çünkü \(j\)'yi kapatabilecek başlangıç noktaları yalnızca bu iki gündür.
Günlerden oluşan bir \(S\) kümesi için
$$N(S)=S\cup(S-1), \qquad S-1=\{i-1 \pmod D : i\in S\}$$
tanımını yapalım.
\(S\) içindeki tüm günler \(n\) seçimden sonra hâlâ açıkta ise, yapılan bütün seçimler \(N(S)\) kümesinden kaçınmak zorundadır. Her seçim \(D\) başlangıç günü arasında eşit olasılıklı olduğundan
$$\Pr(A_S(n))=\left(1-\frac{|N(S)|}{D}\right)^n$$
olur; burada \(A_S(n)\), \(S\)'deki bütün günlerin \(n\) seçim sonunda açıkta kalması olayıdır.
Şimdi açıkta kalan günler üzerinde dahil et-çıkar uygularsak
$$\Pr(T_D>n)=\sum_{\emptyset\ne S\subseteq C_D}(-1)^{|S|+1}\left(1-\frac{|N(S)|}{D}\right)^n$$
elde edilir. Burada \(C_D\), \(D\) günlü çevrimi ifade eder.
Bu ifadeyi kuyruk toplamına koyup toplamaların yerini değiştirince
$$\frac{E(D)}{D}=\sum_{\emptyset\ne S\subseteq C_D}(-1)^{|S|+1}\frac{1}{|N(S)|}$$
çıkar.
Tersini şu kimlikle integrale çevirebiliriz:
$$\frac{1}{m}=\int_0^1 x^{m-1}\,dx.$$
Böylece
$$\frac{E(D)}{D}=\int_0^1\frac{1-Z_D(x)}{x}\,dx,$$
ve
$$Z_D(x)=\sum_{S\subseteq C_D}(-1)^{|S|}x^{|N(S)|}$$
olur. Boş küme \(1\) verdiği için pay kısmında \(1-Z_D(x)\) görünür.
Bir \(S\) kümesini döngüsel ikili dizi \(s_1,\dots,s_D\) ile gösterelim; \(s_i=1\) ise \(i\). günün açıkta kalması isteniyor demektir. O zaman
$$|N(S)|=\sum_{i=1}^{D}s_i+\sum_{i=1}^{D}(1-s_i)s_{i+1}, \qquad s_{D+1}=s_1.$$
İlk toplam seçilen günlerin kendisini sayar; ikinci toplam ise çevrim üzerindeki her ardışık \(1\) bloğu için bir ek öncül gün sayar.
Bu sayede ağırlık yerel olarak çarpanlara ayrılır:
$$(-1)^{|S|}x^{|N(S)|}=\prod_{i=1}^{D} M_{s_i,s_{i+1}}(x),$$
burada transfer matrisi
$$M(x)=\begin{pmatrix}1 & -x^2 \\ 1 & -x\end{pmatrix}$$
şeklindedir. Tüm döngüsel ikili diziler üzerinde toplamak, doğrudan bu matrisin \(D\). kuvvetinin izini verir:
$$Z_D(x)=\operatorname{tr}(M(x)^D).$$
\(M(x)\)'in karakteristik polinomu
$$\lambda^2-(1-x)\lambda-x(1-x)=0$$
olur. İki özdeğer
$$\lambda_{1,2}(x)=\frac{1-x\pm\sqrt{(1-x)(1+3x)}}{2}$$
şeklindedir. Dolayısıyla
$$Z_D(x)=\lambda_1(x)^D+\lambda_2(x)^D$$
ve nihayet
$$\boxed{\frac{E(D)}{D}=\int_0^1\frac{1-\lambda_1(x)^D-\lambda_2(x)^D}{x}\,dx.}$$
Implementasyonun sayısal olarak değerlendirdiği tam formül budur.
Paydada \(x\) olduğu için integral ilk bakışta tekillik varmış gibi görünür, fakat \(x=0\) noktasında bu tekillik kaldırılabilirdir. Küçük \(x\) için
$$\lambda_1(x)=1-x^2+O(x^3), \qquad \lambda_2(x)=-x+O(x^2)$$
olduğundan
$$1-\lambda_1(x)^D-\lambda_2(x)^D = D x^2 + O(x^3)$$
elde edilir. Böylece integrand \(Dx+O(x^2)\) gibi davranır ve \(0\)'a gider. \(x=1\)'de ise her iki özdeğer de \(0\) olur; bu noktadaki değer \(1\)'dir.
Beş gün için transfer matrisi izi
$$Z_5(x)=1-5x^2+5x^3-x^5$$
şeklinde sadeleşir. O halde
$$\frac{E(5)}{5}=\int_0^1\frac{5x^2-5x^3+x^5}{x}\,dx=\int_0^1(5x-5x^2+x^4)\,dx.$$
Terim terim integral alırsak
$$\frac{E(5)}{5}=\frac{5}{2}-\frac{5}{3}+\frac{1}{5}=\frac{31}{30}$$
ve dolayısıyla
$$E(5)=\frac{31}{6}$$
bulunur. Bu, implementasyonların kullandığı kesin kontrol değeridir. En küçük anlamlı çevrim için \(Z_2(x)=1-x^2\) olup \(E(2)=1\) çıkar.
C++, Python ve Java implementasyonları aynı integrandı hesaplar:
$$I_D(x)=\frac{1-\lambda_1(x)^D-\lambda_2(x)^D}{x}.$$
Önce \(\sqrt{(1-x)(1+3x)}\) diskriminantı, sonra iki özdeğer ve en sonunda pay kısmı bulunur. Uç noktalar ayrı ele alınır: \(x=0\) için limit \(0\), \(x=1\) için değer \(1\)'dir.
Sayısal açıdan en hassas terim \(1-\lambda_1(x)^D\)'dir. \(\lambda_1(x)\) değeri \(1\)'e çok yakın olduğunda doğrudan çıkarma ciddi hassasiyet kaybına yol açar. Bu yüzden implementasyon şu kararlı eşdeğeri kullanır:
$$1-\lambda_1(x)^D=-\operatorname{expm1}\!\bigl(D\log \lambda_1(x)\bigr).$$
İkinci özdeğer daha küçük mutlak değere sahip olduğundan doğrudan \(D\). kuvvete yükseltilir.
Daha sonra \([0,1]\) aralığındaki integral, çift sayıda alt aralıkla Simpson kuralı kullanılarak hesaplanır. Başlangıç ağından itibaren aralık sayısı sürekli iki katına çıkarılır ve ardışık iki Simpson tahmini karşılaştırılır:
$$|S_{2N}-S_N|\le \text{rtol}\cdot \max(1,|S_{2N}|).$$
Bu bağıl hata ölçütü sağlandığında sonuç döndürülür. C++ implementasyonu ayrıca ağ yeterince büyük olduğunda iç Simpson düğümlerini paralel toplayabilir ve türetimi \(E(2)=1\), \(E(5)=31/6\) ve \(E(365)\approx 1174.3501\) değerleriyle sınar.
Son Simpson ağı \(N\) alt aralık kullanıyorsa tek bir geçiş \(O(N)\) integrand değerlendirmesi ister. Aralık sayısı geometrik olarak iki katlandığı için kabul edilen ağa ulaşana kadarki toplam iş yine \(O(N)\) olur. Tek iş parçacıklı implementasyonlar \(O(1)\) ek bellek kullanır; paralel C++ sürümü ise \(T\) kısmi toplam için \(O(T)\) ek bellek kullanır.
Colocamos los \(D\) días sobre un ciclo. Cada selección aleatoria escoge un día de forma uniforme y cubre dos días consecutivos: el día elegido y su sucesor cíclico. Sea \(T_D\) el número de selecciones necesarias hasta que todos los días queden cubiertos al menos una vez. La cantidad buscada es
$$E(D)=\mathbb{E}[T_D]$$
para el valor grande \(D=10000\).
Las implementaciones en C++, Python y Java no simulan el proceso. Lo transforman en una integral unidimensional cuyo integrando sale de aplicar inclusión-exclusión a los días no cubiertos y de analizar una matriz de transferencia \(2\times2\) sobre el ciclo.
Numeramos los días módulo \(D\). Una selección en el día \(i\) cubre el par \((i,i+1)\). El tiempo de parada \(T_D\) es el primer instante en el que todos los días ya están cubiertos.
Para cualquier tiempo de parada entero no negativo se cumple
$$E(D)=\sum_{n=0}^{\infty}\Pr(T_D>n).$$
Así que el problema se reduce a calcular la probabilidad de que, después de \(n\) selecciones, todavía exista al menos un día sin cubrir.
Un día \(j\) sigue sin cubrir exactamente cuando ninguna selección cayó en \(j\) ni en \(j-1\), porque esos son los únicos dos puntos de inicio que podrían cubrirlo.
Para un subconjunto \(S\) de días definimos
$$N(S)=S\cup(S-1), \qquad S-1=\{i-1 \pmod D : i\in S\}.$$
Si todos los días de \(S\) siguen sin cubrir tras \(n\) selecciones, entonces cada selección debe evitar el conjunto \(N(S)\). Como cada selección es uniforme entre los \(D\) días de inicio,
$$\Pr(A_S(n))=\left(1-\frac{|N(S)|}{D}\right)^n,$$
donde \(A_S(n)\) es el evento de que todos los días de \(S\) permanecen descubiertos tras \(n\) selecciones.
Ahora la inclusión-exclusión sobre los días no cubiertos produce
$$\Pr(T_D>n)=\sum_{\emptyset\ne S\subseteq C_D}(-1)^{|S|+1}\left(1-\frac{|N(S)|}{D}\right)^n,$$
siendo \(C_D\) el ciclo de \(D\) días.
Al insertar esa expresión en la suma de colas e intercambiar el orden de las sumas obtenemos
$$\frac{E(D)}{D}=\sum_{\emptyset\ne S\subseteq C_D}(-1)^{|S|+1}\frac{1}{|N(S)|}.$$
Convertimos el recíproco usando
$$\frac{1}{m}=\int_0^1 x^{m-1}\,dx.$$
Entonces
$$\frac{E(D)}{D}=\int_0^1\frac{1-Z_D(x)}{x}\,dx,$$
con
$$Z_D(x)=\sum_{S\subseteq C_D}(-1)^{|S|}x^{|N(S)|}.$$
El subconjunto vacío aporta el término inicial \(1\), y por eso aparece \(1-Z_D(x)\) en el numerador.
Representamos un subconjunto \(S\) mediante una palabra binaria cíclica \(s_1,\dots,s_D\), donde \(s_i=1\) significa que el día \(i\) se exige como no cubierto. Entonces
$$|N(S)|=\sum_{i=1}^{D}s_i+\sum_{i=1}^{D}(1-s_i)s_{i+1}, \qquad s_{D+1}=s_1.$$
La primera suma cuenta los propios días marcados; la segunda añade un predecesor extra por cada bloque consecutivo de unos en el ciclo.
Esto hace que el peso se factorice localmente:
$$(-1)^{|S|}x^{|N(S)|}=\prod_{i=1}^{D} M_{s_i,s_{i+1}}(x),$$
con la matriz de transferencia
$$M(x)=\begin{pmatrix}1 & -x^2 \\ 1 & -x\end{pmatrix}.$$
Sumar sobre todas las palabras binarias cíclicas equivale a tomar la traza de la \(D\)-ésima potencia:
$$Z_D(x)=\operatorname{tr}(M(x)^D).$$
El polinomio característico de \(M(x)\) es
$$\lambda^2-(1-x)\lambda-x(1-x)=0.$$
Sus autovalores son
$$\lambda_{1,2}(x)=\frac{1-x\pm\sqrt{(1-x)(1+3x)}}{2}.$$
Por lo tanto
$$Z_D(x)=\lambda_1(x)^D+\lambda_2(x)^D,$$
y la fórmula final queda
$$\boxed{\frac{E(D)}{D}=\int_0^1\frac{1-\lambda_1(x)^D-\lambda_2(x)^D}{x}\,dx.}$$
Esa es exactamente la expresión evaluada por la implementación.
La división por \(x\) parece problemática, pero la singularidad en \(x=0\) es removible. Cerca de \(0\),
$$\lambda_1(x)=1-x^2+O(x^3), \qquad \lambda_2(x)=-x+O(x^2),$$
de modo que
$$1-\lambda_1(x)^D-\lambda_2(x)^D = D x^2 + O(x^3).$$
Así, el integrando se comporta como \(Dx+O(x^2)\) y tiende a \(0\). En \(x=1\), ambos autovalores valen \(0\), así que el integrando vale \(1\).
Para cinco días, la traza de la matriz se simplifica a
$$Z_5(x)=1-5x^2+5x^3-x^5.$$
Entonces
$$\frac{E(5)}{5}=\int_0^1\frac{5x^2-5x^3+x^5}{x}\,dx=\int_0^1(5x-5x^2+x^4)\,dx.$$
Integrando término a término:
$$\frac{E(5)}{5}=\frac{5}{2}-\frac{5}{3}+\frac{1}{5}=\frac{31}{30},$$
y por tanto
$$E(5)=\frac{31}{6}.$$
Ese valor coincide con el punto de control exacto usado por las implementaciones. En el caso mínimo no trivial, \(Z_2(x)=1-x^2\), por lo que \(E(2)=1\).
Las implementaciones en C++, Python y Java evalúan el mismo integrando:
$$I_D(x)=\frac{1-\lambda_1(x)^D-\lambda_2(x)^D}{x}.$$
Primero calculan el discriminante \(\sqrt{(1-x)(1+3x)}\), luego los dos autovalores y finalmente el numerador. Los extremos se tratan aparte: el límite en \(x=0\) es \(0\), y en \(x=1\) el valor es \(1\).
La parte delicada es \(1-\lambda_1(x)^D\). Cuando \(\lambda_1(x)\) está muy cerca de \(1\), la resta directa pierde precisión. Para evitarlo, la implementación usa la identidad estable
$$1-\lambda_1(x)^D=-\operatorname{expm1}\!\bigl(D\log \lambda_1(x)\bigr).$$
El segundo autovalor tiene menor magnitud y se eleva directamente a la potencia \(D\).
Después la integral sobre \([0,1]\) se aproxima con la regla de Simpson usando un número par de subintervalos. A partir de una malla inicial, el código duplica repetidamente el número de intervalos y compara dos estimaciones consecutivas de Simpson:
$$|S_{2N}-S_N|\le \text{rtol}\cdot \max(1,|S_{2N}|).$$
Cuando se satisface ese criterio relativo, se devuelve la estimación refinada. La implementación en C++ también puede paralelizar la suma de los nodos interiores de Simpson cuando la malla es grande y comprueba la derivación con los valores \(E(2)=1\), \(E(5)=31/6\) y \(E(365)\approx 1174.3501\).
Si la malla final de Simpson tiene \(N\) subintervalos, una pasada cuesta \(O(N)\) evaluaciones del integrando. Como la malla se duplica geométricamente, el trabajo total hasta llegar a la malla aceptada sigue siendo \(O(N)\). Las versiones monohilo usan \(O(1)\) memoria auxiliar; la versión paralela en C++ usa \(O(T)\) para \(T\) sumas parciales.
把 \(D\) 个日期排成一个环。每次随机选择都等概率地选中某一天，并覆盖两个相邻日期：被选中的那一天以及它在环上的后继。记 \(T_D\) 为“所有日期都至少被覆盖一次”所需的选择次数，则要求的是
$$E(D)=\mathbb{E}[T_D]$$
在大规模目标 \(D=10000\) 时的数值。
C++、Python 和 Java 实现都没有做蒙特卡洛模拟。它们把这个期望值化成一个一维积分，而这个积分的被积函数来自两个步骤：先对“仍未覆盖的日期集合”做容斥，再把环上的组合计数写成一个 \(2\times2\) 转移矩阵问题。
下面把日期按模 \(D\) 编号。一次落在第 \(i\) 天的选择会覆盖 \((i,i+1)\) 这对相邻日期。停止时间 \(T_D\) 就是第一次所有日期都被覆盖的时刻。
对任何取非负整数值的停止时间，都有经典恒等式
$$E(D)=\sum_{n=0}^{\infty}\Pr(T_D>n).$$
因此问题被改写成：做完 \(n\) 次选择后，至少还有一天没有被覆盖的概率是多少？
某一天 \(j\) 仍未被覆盖，当且仅当没有选择落在 \(j\) 上，也没有选择落在 \(j-1\) 上，因为只有这两个起点会覆盖到 \(j\)。
对任意日期子集 \(S\)，定义
$$N(S)=S\cup(S-1), \qquad S-1=\{i-1 \pmod D : i\in S\}.$$
如果 \(S\) 中所有日期在 \(n\) 次选择后都还未覆盖，那么每一次选择都必须避开集合 \(N(S)\)。由于每次选择都在 \(D\) 个起点中等概率发生，所以
$$\Pr(A_S(n))=\left(1-\frac{|N(S)|}{D}\right)^n,$$
其中 \(A_S(n)\) 表示“\(S\) 中所有日期在 \(n\) 次选择后都未覆盖”的事件。
现在对所有未覆盖日期做容斥，可得
$$\Pr(T_D>n)=\sum_{\emptyset\ne S\subseteq C_D}(-1)^{|S|+1}\left(1-\frac{|N(S)|}{D}\right)^n,$$
这里 \(C_D\) 表示长度为 \(D\) 的环。
把上式代回尾和公式，并交换求和顺序，就得到
$$\frac{E(D)}{D}=\sum_{\emptyset\ne S\subseteq C_D}(-1)^{|S|+1}\frac{1}{|N(S)|}.$$
接着用
$$\frac{1}{m}=\int_0^1 x^{m-1}\,dx$$
把倒数写成积分，因此
$$\frac{E(D)}{D}=\int_0^1\frac{1-Z_D(x)}{x}\,dx,$$
其中
$$Z_D(x)=\sum_{S\subseteq C_D}(-1)^{|S|}x^{|N(S)|}.$$
空集贡献的是首项 \(1\)，这就是分子里出现 \(1-Z_D(x)\) 的原因。
把子集 \(S\) 表示成一个环形二进制串 \(s_1,\dots,s_D\)，其中 \(s_i=1\) 表示第 \(i\) 天被要求保持未覆盖。此时有
$$|N(S)|=\sum_{i=1}^{D}s_i+\sum_{i=1}^{D}(1-s_i)s_{i+1}, \qquad s_{D+1}=s_1.$$
第一项计算被选中的那些日期本身，第二项则给每一段连续的 \(1\) 增加一个额外的前驱位置，所以正好得到需要避开的起点总数。
这样权重就能局部分解：
$$(-1)^{|S|}x^{|N(S)|}=\prod_{i=1}^{D} M_{s_i,s_{i+1}}(x),$$
对应的转移矩阵是
$$M(x)=\begin{pmatrix}1 & -x^2 \\ 1 & -x\end{pmatrix}.$$
对所有环形二进制串求和，恰好等于这个矩阵 \(D\) 次幂的迹：
$$Z_D(x)=\operatorname{tr}(M(x)^D).$$
\(M(x)\) 的特征多项式是
$$\lambda^2-(1-x)\lambda-x(1-x)=0.$$
因此两个特征值为
$$\lambda_{1,2}(x)=\frac{1-x\pm\sqrt{(1-x)(1+3x)}}{2}.$$
于是
$$Z_D(x)=\lambda_1(x)^D+\lambda_2(x)^D,$$
从而得到最终闭式
$$\boxed{\frac{E(D)}{D}=\int_0^1\frac{1-\lambda_1(x)^D-\lambda_2(x)^D}{x}\,dx.}$$
这正是实现中实际计算的公式。
虽然公式里有 \(1/x\)，但 \(x=0\) 只是可去奇点。对小 \(x\) 展开可得
$$\lambda_1(x)=1-x^2+O(x^3), \qquad \lambda_2(x)=-x+O(x^2),$$
所以
$$1-\lambda_1(x)^D-\lambda_2(x)^D = D x^2 + O(x^3).$$
于是被积函数是 \(Dx+O(x^2)\)，极限为 \(0\)。在 \(x=1\) 处，两个特征值都变成 \(0\)，因此被积函数的值正好是 \(1\)。这些端点行为与实现中的特判完全一致。
当 \(D=5\) 时，矩阵迹可以化为
$$Z_5(x)=1-5x^2+5x^3-x^5.$$
于是
$$\frac{E(5)}{5}=\int_0^1\frac{5x^2-5x^3+x^5}{x}\,dx=\int_0^1(5x-5x^2+x^4)\,dx.$$
逐项积分得到
$$\frac{E(5)}{5}=\frac{5}{2}-\frac{5}{3}+\frac{1}{5}=\frac{31}{30},$$
因此
$$E(5)=\frac{31}{6}.$$
这与实现里使用的精确校验值完全一致。最小的非平凡例子 \(D=2\) 对应 \(Z_2(x)=1-x^2\)，所以 \(E(2)=1\)。
C++、Python 和 Java 实现都在数值上计算同一个被积函数
$$I_D(x)=\frac{1-\lambda_1(x)^D-\lambda_2(x)^D}{x}.$$
它们先计算判别式 \(\sqrt{(1-x)(1+3x)}\)，再得到两个特征值，然后构造分子。端点单独处理：\(x=0\) 取极限值 \(0\)，\(x=1\) 取值 \(1\)。
数值上最敏感的部分是 \(1-\lambda_1(x)^D\)。当 \(\lambda_1(x)\) 非常接近 \(1\) 时，直接相减会产生严重的消减误差。因此实现改写为
$$1-\lambda_1(x)^D=-\operatorname{expm1}\!\bigl(D\log \lambda_1(x)\bigr),$$
这样在 \(D\log \lambda_1(x)\) 接近 \(0\) 时仍能保持精度。第二个特征值的绝对值更小，因此直接做 \(D\) 次幂即可。
积分区间 \([0,1]\) 采用 Simpson 求积，子区间数必须为偶数。程序从初始网格开始，不断把区间数翻倍，并比较相邻两次 Simpson 结果：
$$|S_{2N}-S_N|\le \text{rtol}\cdot \max(1,|S_{2N}|).$$
一旦满足这个相对误差判据，就返回较细网格上的值。C++ 实现还会在网格足够大时并行累加 Simpson 内部节点，并检查 \(E(2)=1\)、\(E(5)=31/6\)、\(E(365)\approx 1174.3501\) 这几个参考值。
如果最终接受的 Simpson 网格有 \(N\) 个子区间，那么单次积分扫描需要 \(O(N)\) 次被积函数求值。由于网格按几何级数翻倍，直到达到最终网格为止的总工作量仍然是 \(O(N)\)。单线程实现只需要 \(O(1)\) 额外空间；并行的 C++ 版本需要 \(O(T)\) 额外空间来保存 \(T\) 份部分和。
Рассмотрим \(D\) дней, расположенных по кругу. Один случайный выбор равновероятно выбирает день и покрывает два соседних дня: выбранный день и следующий за ним по циклу. Пусть \(T_D\) обозначает число выборов, необходимое для того, чтобы каждый день оказался покрыт хотя бы один раз. Требуется вычислить
$$E(D)=\mathbb{E}[T_D]$$
для большого значения \(D=10000\).
Реализации на C++, Python и Java не моделируют процесс напрямую. Они сводят математическое ожидание к одномерному интегралу, а сам подынтегральный множитель получается из включения-исключения по непокрытым дням и из \(2\times2\) матрицы перехода на цикле.
Пронумеруем дни по модулю \(D\). Выбор дня \(i\) покрывает пару \((i,i+1)\). Время остановки \(T_D\) — это первый момент, когда покрыты все дни.
Для любого неотрицательного целочисленного времени остановки верно
$$E(D)=\sum_{n=0}^{\infty}\Pr(T_D>n).$$
Значит, задача сводится к вероятности того, что после \(n\) выборов хотя бы один день еще не покрыт.
День \(j\) остается непокрытым тогда и только тогда, когда ни один выбор не попал в \(j\) и ни один выбор не попал в \(j-1\), потому что только эти два старта могут покрыть день \(j\).
Для подмножества дней \(S\) определим
$$N(S)=S\cup(S-1), \qquad S-1=\{i-1 \pmod D : i\in S\}.$$
Если после \(n\) выборов все дни из \(S\) все еще непокрыты, то каждый выбор обязан избегать множества \(N(S)\). Так как каждый выбор равновероятен среди \(D\) стартовых дней, получаем
$$\Pr(A_S(n))=\left(1-\frac{|N(S)|}{D}\right)^n,$$
где \(A_S(n)\) — событие, что все дни из \(S\) остаются непокрытыми после \(n\) выборов.
Теперь включение-исключение по непокрытым дням дает
$$\Pr(T_D>n)=\sum_{\emptyset\ne S\subseteq C_D}(-1)^{|S|+1}\left(1-\frac{|N(S)|}{D}\right)^n,$$
где \(C_D\) — цикл из \(D\) дней.
Подставляя это выражение в сумму хвостов и меняя порядок суммирования, получаем
$$\frac{E(D)}{D}=\sum_{\emptyset\ne S\subseteq C_D}(-1)^{|S|+1}\frac{1}{|N(S)|}.$$
Далее используем представление
$$\frac{1}{m}=\int_0^1 x^{m-1}\,dx.$$
Следовательно,
$$\frac{E(D)}{D}=\int_0^1\frac{1-Z_D(x)}{x}\,dx,$$
где
$$Z_D(x)=\sum_{S\subseteq C_D}(-1)^{|S|}x^{|N(S)|}.$$
Пустое множество дает вклад \(1\), поэтому в числителе появляется выражение \(1-Z_D(x)\).
Представим подмножество \(S\) циклическим двоичным словом \(s_1,\dots,s_D\), где \(s_i=1\) означает, что день \(i\) требуется оставить непокрытым. Тогда
$$|N(S)|=\sum_{i=1}^{D}s_i+\sum_{i=1}^{D}(1-s_i)s_{i+1}, \qquad s_{D+1}=s_1.$$
Первая сумма считает сами отмеченные дни, а вторая добавляет ровно одного предшественника для каждого блока подряд идущих единиц на цикле.
Из-за этого вес факторизуется локально:
$$(-1)^{|S|}x^{|N(S)|}=\prod_{i=1}^{D} M_{s_i,s_{i+1}}(x),$$
где матрица перехода имеет вид
$$M(x)=\begin{pmatrix}1 & -x^2 \\ 1 & -x\end{pmatrix}.$$
Суммирование по всем циклическим двоичным словам эквивалентно следу \(D\)-й степени этой матрицы:
$$Z_D(x)=\operatorname{tr}(M(x)^D).$$
Характеристический многочлен матрицы \(M(x)\) равен
$$\lambda^2-(1-x)\lambda-x(1-x)=0.$$
Ее собственные значения:
$$\lambda_{1,2}(x)=\frac{1-x\pm\sqrt{(1-x)(1+3x)}}{2}.$$
Поэтому
$$Z_D(x)=\lambda_1(x)^D+\lambda_2(x)^D,$$
и итоговая формула становится
$$\boxed{\frac{E(D)}{D}=\int_0^1\frac{1-\lambda_1(x)^D-\lambda_2(x)^D}{x}\,dx.}$$
Именно этот точный интеграл вычисляет реализация.
Деление на \(x\) выглядит подозрительно, но при \(x=0\) здесь только устранимая особенность. Для малых \(x\)
$$\lambda_1(x)=1-x^2+O(x^3), \qquad \lambda_2(x)=-x+O(x^2),$$
поэтому
$$1-\lambda_1(x)^D-\lambda_2(x)^D = D x^2 + O(x^3).$$
Следовательно, подынтегральная функция ведет себя как \(Dx+O(x^2)\) и стремится к \(0\). При \(x=1\) оба собственных значения равны \(0\), так что значение подынтегральной функции равно \(1\).
Для пяти дней след матрицы упрощается до
$$Z_5(x)=1-5x^2+5x^3-x^5.$$
Тогда
$$\frac{E(5)}{5}=\int_0^1\frac{5x^2-5x^3+x^5}{x}\,dx=\int_0^1(5x-5x^2+x^4)\,dx.$$
Интегрируя по частям, получаем
$$\frac{E(5)}{5}=\frac{5}{2}-\frac{5}{3}+\frac{1}{5}=\frac{31}{30},$$
откуда
$$E(5)=\frac{31}{6}.$$
Это точно совпадает с контрольным значением в реализациях. Для \(D=2\) имеем \(Z_2(x)=1-x^2\), значит \(E(2)=1\).
Реализации на C++, Python и Java численно вычисляют один и тот же интегранд
$$I_D(x)=\frac{1-\lambda_1(x)^D-\lambda_2(x)^D}{x}.$$
Сначала находится дискриминант \(\sqrt{(1-x)(1+3x)}\), затем оба собственных значения, а после этого числитель. Концы отрезка обработаны отдельно: при \(x=0\) берется предельное значение \(0\), а при \(x=1\) значение равно \(1\).
Самое чувствительное место — выражение \(1-\lambda_1(x)^D\). Когда \(\lambda_1(x)\) очень близко к \(1\), прямое вычитание теряет значащие цифры. Поэтому реализация использует устойчивое преобразование
$$1-\lambda_1(x)^D=-\operatorname{expm1}\!\bigl(D\log \lambda_1(x)\bigr).$$
Второе собственное значение меньше по модулю, поэтому его возводят в степень напрямую.
Далее интеграл на \([0,1]\) считается формулой Симпсона с четным числом интервалов. Начиная с начальной сетки, код последовательно удваивает число интервалов и сравнивает две соседние оценки Симпсона:
$$|S_{2N}-S_N|\le \text{rtol}\cdot \max(1,|S_{2N}|).$$
Как только это относительное условие выполнено, возвращается уточненное значение. Реализация на C++ дополнительно умеет параллельно суммировать внутренние узлы Симпсона на больших сетках и проверяет вывод на значениях \(E(2)=1\), \(E(5)=31/6\) и \(E(365)\approx 1174.3501\).
Если итоговая сетка Симпсона содержит \(N\) подинтервалов, то один проход требует \(O(N)\) вычислений интегранда. Поскольку число интервалов удваивается геометрически, суммарная работа до принятой сетки тоже имеет порядок \(O(N)\). Однопоточные реализации используют \(O(1)\) дополнительной памяти, а параллельная версия на C++ — \(O(T)\) для \(T\) частичных сумм.
نرتب \(D\) يومًا على دائرة. في كل اختيار عشوائي نختار يومًا بتوزيع منتظم، وهذا الاختيار يغطي يومين متجاورين: اليوم المختار واليوم التالي له على الدورة. لنرمز بـ \(T_D\) إلى عدد الاختيارات اللازمة حتى تصبح جميع الأيام مغطاة مرة واحدة على الأقل. المطلوب هو حساب
$$E(D)=\mathbb{E}[T_D]$$
عندما يكون \(D=10000\).
تنفيذات ++C وPython وJava لا تعتمد على المحاكاة المباشرة. بدلًا من ذلك تحوّل قيمة التوقع إلى تكامل أحادي البعد، ويأتي هذا التكامل من مبدأ الاشتمال والاستبعاد على مجموعة الأيام غير المغطاة، ومن مصفوفة انتقال \(2\times2\) على الدورة.
نرقم الأيام بترديد \(D\). إذا وقع الاختيار على اليوم \(i\)، فإنه يغطي الزوج \((i,i+1)\). وزمن التوقف \(T_D\) هو أول لحظة تصبح فيها جميع الأيام مغطاة.
لكل زمن توقف غير سالب ذي قيم صحيحة لدينا
$$E(D)=\sum_{n=0}^{\infty}\Pr(T_D>n).$$
إذًا جوهر المسألة هو إيجاد احتمال أن يبقى يوم واحد على الأقل غير مغطى بعد \(n\) اختيارات.
ويبقى اليوم \(j\) غير مغطى إذا وفقط إذا لم يقع أي اختيار على \(j\) ولم يقع أي اختيار على \(j-1\)، لأن هاتين نقطتي البداية الوحيدتين القادرتين على تغطيته.
لكل مجموعة أيام \(S\) نعرّف
$$N(S)=S\cup(S-1), \qquad S-1=\{i-1 \pmod D : i\in S\}.$$
إذا كانت جميع الأيام في \(S\) ما تزال غير مغطاة بعد \(n\) اختيارات، فلا بد أن تتجنب كل الاختيارات المجموعة \(N(S)\). وبما أن كل اختيار منتظم بين \(D\) بدايات ممكنة، فإن
$$\Pr(A_S(n))=\left(1-\frac{|N(S)|}{D}\right)^n,$$
حيث يعبّر \(A_S(n)\) عن حدث بقاء كل أيام \(S\) غير مغطاة بعد \(n\) اختيارات.
الآن يعطي الاشتمال والاستبعاد على الأيام غير المغطاة الصيغة
$$\Pr(T_D>n)=\sum_{\emptyset\ne S\subseteq C_D}(-1)^{|S|+1}\left(1-\frac{|N(S)|}{D}\right)^n,$$
حيث \(C_D\) هي الدورة المؤلفة من \(D\) يومًا.
بالتعويض في مجموع الذيول وتبديل ترتيب الجمع نحصل على
$$\frac{E(D)}{D}=\sum_{\emptyset\ne S\subseteq C_D}(-1)^{|S|+1}\frac{1}{|N(S)|}.$$
ثم نستخدم الهوية
$$\frac{1}{m}=\int_0^1 x^{m-1}\,dx$$
لتحويل المقلوب إلى تكامل. ومن ثم
$$\frac{E(D)}{D}=\int_0^1\frac{1-Z_D(x)}{x}\,dx,$$
حيث
$$Z_D(x)=\sum_{S\subseteq C_D}(-1)^{|S|}x^{|N(S)|}.$$
المجموعة الخالية تعطي الحد \(1\)، ولذلك يظهر في البسط التعبير \(1-Z_D(x)\).
نمثل المجموعة \(S\) بكلمة ثنائية دورية \(s_1,\dots,s_D\)، حيث \(s_i=1\) يعني أن اليوم \(i\) مطلوب أن يبقى غير مغطى. عندها
$$|N(S)|=\sum_{i=1}^{D}s_i+\sum_{i=1}^{D}(1-s_i)s_{i+1}, \qquad s_{D+1}=s_1.$$
المجموع الأول يحسب الأيام المحددة نفسها، والمجموع الثاني يضيف سلفًا واحدًا لكل كتلة متتالية من القيم \(1\) على الدورة.
وبذلك يتفكك الوزن إلى حاصل ضرب محلي:
$$(-1)^{|S|}x^{|N(S)|}=\prod_{i=1}^{D} M_{s_i,s_{i+1}}(x),$$
مع مصفوفة انتقال تساوي
$$M(x)=\begin{pmatrix}1 & -x^2 \\ 1 & -x\end{pmatrix}.$$
وجمع الأوزان على جميع الكلمات الثنائية الدورية يساوي بالضبط أثر القوة \(D\) لهذه المصفوفة:
$$Z_D(x)=\operatorname{tr}(M(x)^D).$$
كثير الحدود المميز للمصفوفة \(M(x)\) هو
$$\lambda^2-(1-x)\lambda-x(1-x)=0.$$
ومن ثم تكون القيمتان الذاتيتان
$$\lambda_{1,2}(x)=\frac{1-x\pm\sqrt{(1-x)(1+3x)}}{2}.$$
إذًا
$$Z_D(x)=\lambda_1(x)^D+\lambda_2(x)^D,$$
ونصل إلى الصيغة النهائية
$$\boxed{\frac{E(D)}{D}=\int_0^1\frac{1-\lambda_1(x)^D-\lambda_2(x)^D}{x}\,dx.}$$
وهذه هي الصيغة الدقيقة التي تحسبها التنفيذات عدديًا.
قد يبدو أن القسمة على \(x\) تسبب مشكلة، لكن عند \(x=0\) تكون لدينا مفردة قابلة للإزالة فقط. عندما يكون \(x\) صغيرًا نحصل على
$$\lambda_1(x)=1-x^2+O(x^3), \qquad \lambda_2(x)=-x+O(x^2),$$
ومن ثم
$$1-\lambda_1(x)^D-\lambda_2(x)^D = D x^2 + O(x^3).$$
إذن integrand يتصرف مثل \(Dx+O(x^2)\) ويؤول إلى الصفر. وعند \(x=1\) تختفي كلتا القيمتين الذاتيتين، فيكون مقدار integrand مساويًا لـ \(1\). هذا يفسر المعالجة الخاصة للطرفين في التنفيذ.
عندما \(D=5\)، يتبسط أثر المصفوفة إلى
$$Z_5(x)=1-5x^2+5x^3-x^5.$$
وعليه
$$\frac{E(5)}{5}=\int_0^1\frac{5x^2-5x^3+x^5}{x}\,dx=\int_0^1(5x-5x^2+x^4)\,dx.$$
وبالتكامل حدًا حدًا نجد
$$\frac{E(5)}{5}=\frac{5}{2}-\frac{5}{3}+\frac{1}{5}=\frac{31}{30},$$
أي أن
$$E(5)=\frac{31}{6}.$$
وهذا يطابق تمامًا قيمة الفحص المعروفة في التنفيذات. أما في الحالة الأصغر غير التافهة \(D=2\)، فلدينا \(Z_2(x)=1-x^2\) وبالتالي \(E(2)=1\).
تنفيذات ++C وPython وJava كلها تقيّم نفس integrand العددي
$$I_D(x)=\frac{1-\lambda_1(x)^D-\lambda_2(x)^D}{x}.$$
تبدأ بحساب المميز \(\sqrt{(1-x)(1+3x)}\)، ثم القيمتين الذاتيتين، ثم البسط. ويُتعامل مع الطرفين مباشرة: القيمة الحدية عند \(x=0\) هي \(0\)، وعند \(x=1\) هي \(1\).
الجزء الحساس عدديًا هو \(1-\lambda_1(x)^D\). عندما تكون \(\lambda_1(x)\) قريبة جدًا من \(1\)، فإن الطرح المباشر يفقد عددًا من الخانات المهمة. لذلك تستخدم التنفيذات الصيغة المستقرة
$$1-\lambda_1(x)^D=-\operatorname{expm1}\!\bigl(D\log \lambda_1(x)\bigr).$$
أما القيمة الذاتية الثانية فهي أصغر من حيث المقدار، ولذلك تُرفع إلى القوة \(D\) مباشرة.
بعد ذلك يُحسب التكامل على \([0,1]\) بقاعدة Simpson مع عدد زوجي من الفواصل. يبدأ البرنامج بشبكة أولية، ثم يضاعف عدد الفواصل مرارًا ويقارن تقريبين متتاليين:
$$|S_{2N}-S_N|\le \text{rtol}\cdot \max(1,|S_{2N}|).$$
متى تحقق هذا الشرط النسبي أعاد القيمة المحسّنة. كما أن تنفيذ ++C يستطيع جمع نقاط Simpson الداخلية على نحو متوازٍ عندما تصبح الشبكة كبيرة، ويتحقق من النتائج المعروفة \(E(2)=1\)، و\(E(5)=31/6\)، و\(E(365)\approx 1174.3501\).
إذا احتوت شبكة Simpson النهائية على \(N\) فاصلًا فرعيًا، فإن مرورًا واحدًا يحتاج إلى \(O(N)\) من تقييمات integrand. وبما أن الشبكة تتضاعف هندسيًا، فإن مجموع العمل حتى الشبكة المقبولة يبقى من الرتبة \(O(N)\). التنفيذات أحادية الخيط تستخدم ذاكرة إضافية \(O(1)\)، أما النسخة المتوازية في ++C فتحتاج \(O(T)\) لتجميع \(T\) مجاميع جزئية.