Problem Summary

Birthdays are modeled as independent uniform choices on a cyclic year of \(N\) days. Let \(T\) be the first time that, among the people already seen, there exist \(k\) birthdays contained in some block of \(D+1\) consecutive calendar days. The task is to compute the exact expectation \(\mathbb E[T]\) for the Earth case \((N,k,D)=(365,4,7)\), while also matching smaller checkpoint cases.

Mathematical Approach

The implementations do not simulate random trials. Instead they count, for each possible population size \(n\), all birthday placements that still avoid a forbidden \(k\)-cluster, and then convert those counts into survival probabilities.

Step 1: Rewrite the expectation as a survival sum

Because \(T\) is a positive integer-valued stopping time, we may use the standard identity

$$\mathbb E[T]=\sum_{n\ge 0}\Pr(T>n).$$

So the problem becomes: for each \(n\), what is the probability that after placing \(n\) people, no set of \(k\) birthdays lies inside any span of \(D\) days?

Step 2: Convert the event into a cyclic window constraint

Let \(c_i\) be the number of birthdays on day \(i\), with indices interpreted cyclically modulo \(N\). Write

$$w=D+1.$$

Then the event \(T>n\) is equivalent to saying that every cyclic window of length \(w\) contains at most \(k-1\) people:

$$c_i+c_{i+1}+\cdots+c_{i+D}\le k-1 \qquad \text{for every } i.$$

This formulation is exact, including windows that wrap across the end of the year, because the year is treated as circular.

Step 3: Count valid occupancy vectors with exponential weights

Fix \(n\). A vector \((c_1,\dots,c_N)\) with total \(n\) represents how many of the \(n\) labeled people land on each day. The number of assignments producing that occupancy vector is the multinomial count

$$\frac{n!}{c_1!\cdots c_N!}.$$

If \(\mathcal V_n\) denotes the set of occupancy vectors satisfying the window constraint and \(c_1+\cdots+c_N=n\), then

$$\Pr(T>n)=\frac{n!}{N^n}\sum_{(c_1,\dots,c_N)\in \mathcal V_n}\frac{1}{c_1!\cdots c_N!}.$$

This is why exponential generating functions appear naturally: the weight of placing \(c\) people on one day is \(x^c/c!\).

Step 4: Use a finite state space based on the last \(D\) days

When we move from one day to the next, only the previous \(D\) day-counts matter. A state can therefore be written as

$$s=(a_1,\dots,a_D),\qquad a_j\ge 0,\qquad a_1+\cdots+a_D\le k-1.$$

It records the counts on the last \(D\) days before the new day is appended. If the current state sum is \(\sigma\), then the next day may receive any count \(c\) with

$$0\le c\le k-1-\sigma.$$

The next state is obtained by dropping the oldest entry and appending \(c\):

$$s'=(a_2,\dots,a_D,c).$$

The generating-function weight of this transition is

$$\frac{x^c}{c!}.$$

Step 5: Enforce the cyclic year by taking a trace

If we assemble the transition weights into a transfer matrix \(M(x)\), then one pass over \(N\) days gives \(M(x)^N\). However, windows crossing New Year must obey the same rule as internal windows, so the profile at the end of the year must match the profile at the beginning. That is exactly the diagonal condition encoded by the trace:

$$G(x)=\operatorname{tr}(M(x)^N).$$

If

$$A_n=[x^n]\,G(x),$$

then \(A_n\) is the total exponential weight of all valid cyclic birthday-count configurations with \(n\) people, and therefore

$$\Pr(T>n)=A_n\frac{n!}{N^n}.$$

Step 6: The sum is finite

There is a clean upper bound on the possible population size before a violation becomes unavoidable. Summing the \(N\) cyclic window inequalities gives

$$w(c_1+\cdots+c_N)\le N(k-1),$$

because each person is counted in exactly \(w\) windows. Hence any valid configuration must satisfy

$$n\le \left\lfloor\frac{N(k-1)}{D+1}\right\rfloor.$$

So the expectation is actually a finite sum:

$$\boxed{\mathbb E[T]=\sum_{n=0}^{\left\lfloor N(k-1)/(D+1)\right\rfloor} A_n\frac{n!}{N^n}.}$$

Worked Example: \((N,k,D)=(10,3,1)\)

Here \(w=2\), so survival means that every pair of consecutive days contains at most \(2\) people. A state consists of the count on the immediately previous day, so the possible states are \(0,1,2\).

From state \(0\), the next day may contain \(0,1,\) or \(2\) people. From state \(1\), the next day may contain \(0\) or \(1\) person. From state \(2\), the next day must contain \(0\) people. With the state order \((0,1,2)\), the transfer matrix is

$$M(x)=\begin{pmatrix} 1 & x & x^2/2 \\ 1 & x & 0 \\ 1 & 0 & 0 \end{pmatrix}.$$

The cyclic generating function is \(G(x)=\operatorname{tr}(M(x)^{10})\). Extracting coefficients \(A_n\), converting them by \(A_n n!/10^n\), and summing over \(n\) yields

$$\mathbb E[T]=5.78688636,$$

which is the checkpoint reproduced by the implementations.

How the Code Works

The C++, Python, and Java implementations all follow the same recurrence. They first precompute the reciprocals \(1/c!\) for \(0\le c\le k-1\), because those are the transition weights in the exponential generating function. Next they enumerate every valid profile of length \(D\) whose total is at most \(k-1\), and for each profile they precompute all allowed next profiles together with the number of newly added people.

After that, the implementation runs a layered dynamic program for exactly \(N\) days. For a fixed starting profile it stores, for every current profile and every total population \(n\), the coefficient accumulated so far. After the \(N\)-th day it reads only the coefficient vector that returned to the same starting profile; summing these diagonal contributions over all starting profiles is exactly the trace \(\operatorname{tr}(M(x)^N)\). Finally, it multiplies the resulting coefficient of degree \(n\) by \(n!/N^n\) and adds the terms to obtain \(\mathbb E[T]\). The C++ implementation also distributes different starting profiles across worker threads, but the mathematical recurrence is identical in all three languages.

Complexity Analysis

The number of states is

$$S=\#\{(a_1,\dots,a_D)\in \mathbb Z_{\ge 0}^D: a_1+\cdots+a_D\le k-1\}=\binom{D+k-1}{D}.$$

Each state has at most \(k\) outgoing transitions, since the next day count can range only from \(0\) to \(k-1-\sigma\). The straightforward trace computation used here runs the \(N\)-day dynamic program once for each starting state, so its time cost is

$$O\!\left(N\cdot S^2\cdot k\cdot \left\lfloor\frac{N(k-1)}{D+1}\right\rfloor\right).$$

The memory usage per worker is

$$O\!\left(S\cdot \left\lfloor\frac{N(k-1)}{D+1}\right\rfloor\right),$$

because only the current and next DP layers are stored.

Footnotes and References

  1. Problem page: https://projecteuler.net/problem=584
  2. Birthday problem background: Wikipedia - Birthday problem
  3. Expectation as a survival sum: Wikipedia - Expected value
  4. Exponential generating functions: Wikipedia - Exponential generating function
  5. Transfer-matrix methods: Wikipedia - Transfer-matrix method

Problemzusammenfassung

Geburtstage werden als unabhängige gleichverteilte Zufallswerte auf einem zyklischen Jahr mit \(N\) Tagen modelliert. \(T\) sei der erste Zeitpunkt, an dem es unter den bisher eingetroffenen Personen \(k\) Geburtstage gibt, die alle in einem Block von \(D+1\) aufeinanderfolgenden Kalendertagen liegen. Gesucht ist der exakte Erwartungswert \(\mathbb E[T]\) für den Erd-Fall \((N,k,D)=(365,4,7)\), wobei auch kleinere Kontrollfälle exakt getroffen werden.

Mathematischer Ansatz

Die Implementierungen simulieren nicht, sondern zählen für jede mögliche Personenzahl \(n\) alle Geburtstagsbelegungen, die noch keinen verbotenen \(k\)-Cluster enthalten. Daraus werden dann die Überlebenswahrscheinlichkeiten berechnet.

Schritt 1: Erwartungswert als Summe von Überlebenswahrscheinlichkeiten

Da \(T\) eine positive ganzzahlige Stoppzeit ist, gilt

$$\mathbb E[T]=\sum_{n\ge 0}\Pr(T>n).$$

Damit reduziert sich die Aufgabe auf die Frage: Wie groß ist für jedes \(n\) die Wahrscheinlichkeit, dass nach \(n\) Personen noch kein Satz von \(k\) Geburtstagen innerhalb eines Zeitraums von \(D\) Tagen aufgetreten ist?

Schritt 2: Formuliere das Ereignis als zyklische Fensterbedingung

Sei \(c_i\) die Anzahl der Geburtstage am Tag \(i\), wobei die Indizes zyklisch modulo \(N\) gelesen werden. Setze

$$w=D+1.$$

Dann ist \(T>n\) genau dann erfüllt, wenn jedes zyklische Fenster der Länge \(w\) höchstens \(k-1\) Personen enthält:

$$c_i+c_{i+1}+\cdots+c_{i+D}\le k-1 \qquad \text{für alle } i.$$

Diese Formulierung behandelt auch die Fenster über den Jahreswechsel hinweg korrekt, weil das Jahr als Kreis betrachtet wird.

Schritt 3: Gültige Belegungsvektoren mit exponentiellen Gewichten zählen

Fixiere \(n\). Ein Vektor \((c_1,\dots,c_N)\) mit Summe \(n\) beschreibt, wie viele der \(n\) markierten Personen auf die einzelnen Tage fallen. Die Anzahl der Zuordnungen mit genau diesem Belegungsmuster ist

$$\frac{n!}{c_1!\cdots c_N!}.$$

Ist \(\mathcal V_n\) die Menge aller solchen Vektoren, die zusätzlich die Fensterbedingung erfüllen, so gilt

$$\Pr(T>n)=\frac{n!}{N^n}\sum_{(c_1,\dots,c_N)\in \mathcal V_n}\frac{1}{c_1!\cdots c_N!}.$$

Darum sind exponentielle Erzeugendenfunktionen das passende Werkzeug: Ein Tag mit \(c\) Personen trägt das Gewicht \(x^c/c!\).

Schritt 4: Endlicher Zustandsraum über die letzten \(D\) Tage

Beim Übergang zum nächsten Tag sind nur die Belegungen der letzten \(D\) Tage relevant. Ein Zustand hat also die Form

$$s=(a_1,\dots,a_D),\qquad a_j\ge 0,\qquad a_1+\cdots+a_D\le k-1.$$

Er speichert die Besetzungen der letzten \(D\) Tage vor dem neuen Tag. Hat der aktuelle Zustand die Summe \(\sigma\), dann darf der nächste Tag jeden Wert \(c\) annehmen mit

$$0\le c\le k-1-\sigma.$$

Der Folgezustand entsteht, indem der älteste Eintrag entfernt und \(c\) angehängt wird:

$$s'=(a_2,\dots,a_D,c).$$

Das Gewicht dieses Übergangs in der Erzeugendenfunktion ist

$$\frac{x^c}{c!}.$$

Schritt 5: Zyklisches Jahr durch die Spur erzwingen

Fasst man die Übergänge in einer Transfermatrix \(M(x)\) zusammen, dann beschreibt ein Durchlauf über \(N\) Tage die Matrix \(M(x)^N\). Fenster, die über den Jahreswechsel laufen, müssen dieselbe Bedingung erfüllen wie innere Fenster; deshalb muss das Endprofil mit dem Anfangsprofil übereinstimmen. Genau das erfasst die Spur:

$$G(x)=\operatorname{tr}(M(x)^N).$$

Mit

$$A_n=[x^n]\,G(x)$$

ist \(A_n\) das gesamte exponentielle Gewicht aller gültigen zyklischen Belegungen mit \(n\) Personen, also

$$\Pr(T>n)=A_n\frac{n!}{N^n}.$$

Schritt 6: Die Summe ist endlich

Es gibt eine direkte obere Schranke für die mögliche Personenzahl vor der ersten unvermeidlichen Verletzung. Summiert man die \(N\) Fensterungleichungen, erhält man

$$w(c_1+\cdots+c_N)\le N(k-1),$$

denn jede Person wird in genau \(w\) Fenstern mitgezählt. Daher muss jede gültige Konfiguration

$$n\le \left\lfloor\frac{N(k-1)}{D+1}\right\rfloor$$

erfüllen. Somit wird der Erwartungswert zu einer endlichen Summe:

$$\boxed{\mathbb E[T]=\sum_{n=0}^{\left\lfloor N(k-1)/(D+1)\right\rfloor} A_n\frac{n!}{N^n}.}$$

Durchgerechnetes Beispiel: \((N,k,D)=(10,3,1)\)

Hier ist \(w=2\), also bedeutet Überleben: Jedes Paar benachbarter Tage enthält höchstens \(2\) Personen. Ein Zustand besteht nur aus der Belegung des unmittelbar vorherigen Tages, also sind die möglichen Zustände \(0,1,2\).

Aus Zustand \(0\) darf der nächste Tag \(0,1\) oder \(2\) Personen enthalten. Aus Zustand \(1\) sind \(0\) oder \(1\) erlaubt. Aus Zustand \(2\) ist nur \(0\) erlaubt. In der Zustandsreihenfolge \((0,1,2)\) ergibt sich die Transfermatrix

$$M(x)=\begin{pmatrix} 1 & x & x^2/2 \\ 1 & x & 0 \\ 1 & 0 & 0 \end{pmatrix}.$$

Die zyklische Erzeugendenfunktion ist \(G(x)=\operatorname{tr}(M(x)^{10})\). Wenn man die Koeffizienten \(A_n\) bestimmt, mit \(A_n n!/10^n\) in Wahrscheinlichkeiten umrechnet und über \(n\) summiert, erhält man

$$\mathbb E[T]=5.78688636,$$

genau den Kontrollwert der Implementierungen.

Wie der Code arbeitet

Die C++-, Python- und Java-Implementierungen verwenden dieselbe Rekursion. Zuerst werden die Kehrwerte \(1/c!\) für \(0\le c\le k-1\) vorbereitet, weil dies die Übergangsgewichte der exponentiellen Erzeugendenfunktion sind. Danach werden alle zulässigen Profile der Länge \(D\) mit Gesamtsumme höchstens \(k-1\) erzeugt, und für jedes Profil werden alle erlaubten Nachfolger samt Anzahl der neu hinzugefügten Personen vorab gespeichert.

Anschließend läuft ein geschichtetes dynamisches Programm über genau \(N\) Tage. Für ein festes Startprofil speichert es für jedes aktuelle Profil und jede Gesamtpersonenzahl \(n\) den bisher angesammelten Koeffizienten. Nach dem \(N\)-ten Tag wird nur der Koeffizientenvektor ausgelesen, der wieder im gleichen Startprofil endet; die Summe über alle Startprofile ist exakt die Spur \(\operatorname{tr}(M(x)^N)\). Zum Schluss wird der Koeffizient vom Grad \(n\) mit \(n!/N^n\) multipliziert und alles zu \(\mathbb E[T]\) aufsummiert. Die C++-Version verteilt verschiedene Startprofile zusätzlich auf mehrere Threads, mathematisch bleibt das Verfahren aber in allen drei Sprachen identisch.

Komplexitätsanalyse

Die Zahl der Zustände ist

$$S=\#\{(a_1,\dots,a_D)\in \mathbb Z_{\ge 0}^D: a_1+\cdots+a_D\le k-1\}=\binom{D+k-1}{D}.$$

Jeder Zustand hat höchstens \(k\) ausgehende Übergänge, weil die Belegung des nächsten Tages nur zwischen \(0\) und \(k-1-\sigma\) liegen kann. Die hier benutzte direkte Spur-Berechnung führt das \(N\)-Tage-DP für jedes Startprofil einmal aus. Daraus ergibt sich der Zeitaufwand

$$O\!\left(N\cdot S^2\cdot k\cdot \left\lfloor\frac{N(k-1)}{D+1}\right\rfloor\right).$$

Der Speicherbedarf pro Worker ist

$$O\!\left(S\cdot \left\lfloor\frac{N(k-1)}{D+1}\right\rfloor\right),$$

weil nur die aktuelle und die nächste DP-Schicht gehalten werden.

Fußnoten und Quellen

  1. Problemseite: https://projecteuler.net/problem=584
  2. Hintergrund zum Birthday Problem: Wikipedia - Birthday problem
  3. Erwartungswert als Überlebenssumme: Wikipedia - Expected value
  4. Exponentiale Erzeugendenfunktionen: Wikipedia - Exponential generating function
  5. Transfermatrix-Methode: Wikipedia - Transfer-matrix method

Problem Özeti

Doğum günleri, \(N\) günlük döngüsel bir yıl üzerinde bağımsız ve eşit olasılıklı seçimler olarak modellenir. \(T\), o ana kadar gelen insanlar arasında \(D+1\) ardışık takvim günü içinde kalan \(k\) doğum gününün ilk kez görüldüğü giriş sırası olsun. Amaç, Dünya örneği \((N,k,D)=(365,4,7)\) için \(\mathbb E[T]\) değerini tam olarak hesaplamak ve daha küçük kontrol örneklerini de doğru üretmektir.

Matematiksel Yaklaşım

İmplementasyonlar rastgele benzetim yapmaz. Bunun yerine her olası \(n\) kişi sayısı için henüz yasaklı \(k\)-kümesini oluşturmayan tüm doğum günü yerleşimlerini sayar, sonra bu sayımları hayatta kalma olasılıklarına dönüştürür.

Adım 1: Beklentiyi hayatta kalma toplamı olarak yaz

\(T\) pozitif tamsayı değerli bir durma zamanı olduğundan şu özdeşlik geçerlidir:

$$\mathbb E[T]=\sum_{n\ge 0}\Pr(T>n).$$

Dolayısıyla asıl soru şudur: Her \(n\) için, ilk \(n\) kişiden sonra hâlâ \(D\) gün içinde toplanmış \(k\) doğum günü oluşmamış olma olasılığı nedir?

Adım 2: Olayı döngüsel pencere kısıtına çevir

\(c_i\), \(i\). gündeki doğum günü sayısı olsun; indisler \(N\) modunda döngüsel okunur. Şunu tanımlayalım:

$$w=D+1.$$

O zaman \(T>n\) olayı, uzunluğu \(w\) olan her döngüsel pencerenin en fazla \(k-1\) kişi içermesine denktir:

$$c_i+c_{i+1}+\cdots+c_{i+D}\le k-1 \qquad \text{her } i \text{ için}.$$

Yıl sonunu aşan pencereler de bu tanımın içindedir; çünkü takvim bir çember olarak ele alınır.

Adım 3: Geçerli doluluk vektörlerini üstel ağırlıklarla say

\(n\) sabit olsun. Toplamı \(n\) olan \((c_1,\dots,c_N)\) vektörü, etiketli \(n\) kişinin hangi günlere düştüğünü gösterir. Bu doluluk vektörünü veren atama sayısı çokterimli katsayıdır:

$$\frac{n!}{c_1!\cdots c_N!}.$$

\(\mathcal V_n\), hem toplamı \(n\) olan hem de pencere kısıtını sağlayan tüm vektörler kümesi olsun. O zaman

$$\Pr(T>n)=\frac{n!}{N^n}\sum_{(c_1,\dots,c_N)\in \mathcal V_n}\frac{1}{c_1!\cdots c_N!}.$$

Üstel üreteç fonksiyonunun nedeni tam olarak budur: bir güne \(c\) kişi yerleştirmenin ağırlığı \(x^c/c!\) olur.

Adım 4: Son \(D\) güne dayalı sonlu durum uzayı kur

Bir günden sonraki güne geçerken yalnızca son \(D\) günün sayıları önemlidir. Bu yüzden durum şu şekilde yazılabilir:

$$s=(a_1,\dots,a_D),\qquad a_j\ge 0,\qquad a_1+\cdots+a_D\le k-1.$$

Bu vektör, yeni gün eklenmeden önceki son \(D\) günün doluluklarını tutar. Durum toplamı \(\sigma\) ise yeni gün için seçilebilecek \(c\) değerleri

$$0\le c\le k-1-\sigma$$

aralığındadır. Yeni durum, en eski değeri atıp \(c\)'yi sona ekleyerek elde edilir:

$$s'=(a_2,\dots,a_D,c).$$

Bu geçişin üreteç fonksiyonundaki ağırlığı

$$\frac{x^c}{c!}$$

şeklindedir.

Adım 5: Döngüsel yılı iz işlemiyle dayat

Geçişleri \(M(x)\) adlı bir transfer matrisi içinde toplarsak, \(N\) gün boyunca ilerlemek \(M(x)^N\) verir. Ancak yıl sonunu geçen pencerelerin de aynı kuralı sağlaması gerekir; bu yüzden yıl sonundaki profil yıl başındaki profille aynı olmalıdır. Bu koşul tam olarak izin diyagonalinden, yani izden gelir:

$$G(x)=\operatorname{tr}(M(x)^N).$$

Eğer

$$A_n=[x^n]\,G(x)$$

ise, \(A_n\) \(n\) kişilik tüm geçerli döngüsel doğum günü yapılandırmalarının toplam üstel ağırlığıdır. Dolayısıyla

$$\Pr(T>n)=A_n\frac{n!}{N^n}.$$

Adım 6: Toplam aslında sonludur

İhlal kaçınılmaz olmadan önce mümkün olan kişi sayısına doğrudan bir üst sınır vardır. \(N\) pencere eşitsizliğinin hepsini toplarsak

$$w(c_1+\cdots+c_N)\le N(k-1)$$

elde ederiz; çünkü her kişi tam olarak \(w\) pencereye katkı yapar. Bu nedenle her geçerli konfigürasyon

$$n\le \left\lfloor\frac{N(k-1)}{D+1}\right\rfloor$$

koşulunu sağlamalıdır. Böylece beklenti sonlu bir toplam olur:

$$\boxed{\mathbb E[T]=\sum_{n=0}^{\left\lfloor N(k-1)/(D+1)\right\rfloor} A_n\frac{n!}{N^n}.}$$

Çözümlü Örnek: \((N,k,D)=(10,3,1)\)

Burada \(w=2\) olduğu için hayatta kalma koşulu şudur: art arda gelen her iki gün toplamda en fazla \(2\) kişi içersin. Durum yalnızca bir önceki günün sayısından oluşur; dolayısıyla olası durumlar \(0,1,2\)'dir.

Durum \(0\)'dan sonra yeni gün \(0,1\) veya \(2\) kişi alabilir. Durum \(1\)'den sonra yalnızca \(0\) veya \(1\) kişi eklenebilir. Durum \(2\)'den sonra yeni gün zorunlu olarak \(0\) olur. Durum sırası \((0,1,2)\) seçilirse transfer matrisi

$$M(x)=\begin{pmatrix} 1 & x & x^2/2 \\ 1 & x & 0 \\ 1 & 0 & 0 \end{pmatrix}$$

şeklindedir. Döngüsel üreteç fonksiyonu \(G(x)=\operatorname{tr}(M(x)^{10})\) olur. Buradan çıkan \(A_n\) katsayılarını \(A_n n!/10^n\) ile olasılığa çevirip topladığımızda

$$\mathbb E[T]=5.78688636$$

elde edilir; bu da implementasyonların ürettiği kontrol değeridir.

Kod Nasıl Çalışır

C++, Python ve Java implementasyonları aynı yinelemeyi uygular. Önce \(0\le c\le k-1\) için \(1/c!\) değerleri önceden hesaplanır; çünkü bunlar üstel üreteç fonksiyonundaki geçiş ağırlıklarıdır. Sonra toplamı en fazla \(k-1\) olan tüm uzunluk-\(D\) profilleri üretilir ve her profil için hangi sonraki profillere gidilebildiğiyle birlikte yeni eklenen kişi sayısı saklanır.

Bundan sonra tam \(N\) gün süren katmanlı bir dinamik program çalıştırılır. Sabit bir başlangıç profili için, her anki profil ve her toplam kişi sayısı \(n\) adına o ana kadar biriken katsayı tutulur. \(N\). günün sonunda yalnızca yeniden aynı başlangıç profiline dönen katsayı vektörü okunur; tüm başlangıç profilleri üzerinden toplamak tam olarak \(\operatorname{tr}(M(x)^N)\) demektir. Son aşamada derece \(n\) katsayısı \(n!/N^n\) ile çarpılır ve bütün terimler toplanarak \(\mathbb E[T]\) elde edilir. C++ implementasyonu başlangıç profillerini iş parçacıklarına dağıtır; fakat matematiksel yineleme üç dilde de aynıdır.

Karmaşıklık Analizi

Durum sayısı

$$S=\#\{(a_1,\dots,a_D)\in \mathbb Z_{\ge 0}^D: a_1+\cdots+a_D\le k-1\}=\binom{D+k-1}{D}$$

olur. Her durumun en fazla \(k\) çıkışı vardır; çünkü yeni gün sayısı yalnızca \(0\) ile \(k-1-\sigma\) arasında değişebilir. Burada kullanılan doğrudan iz hesabı, \(N\) günlük DP'yi her başlangıç profili için bir kez çalıştırır. Bu yüzden zaman maliyeti

$$O\!\left(N\cdot S^2\cdot k\cdot \left\lfloor\frac{N(k-1)}{D+1}\right\rfloor\right)$$

mertebesindedir. İşçi başına bellek kullanımı ise

$$O\!\left(S\cdot \left\lfloor\frac{N(k-1)}{D+1}\right\rfloor\right)$$

kadardır; çünkü yalnızca mevcut ve bir sonraki DP katmanı tutulur.

Dipnotlar ve Kaynaklar

  1. Problem sayfası: https://projecteuler.net/problem=584
  2. Birthday problem arka planı: Wikipedia - Birthday problem
  3. Beklentinin hayatta kalma toplamı olarak yazılması: Wikipedia - Expected value
  4. Üstel üreteç fonksiyonları: Wikipedia - Exponential generating function
  5. Transfer matrisi yöntemi: Wikipedia - Transfer-matrix method

Resumen del Problema

Los cumpleaños se modelan como elecciones independientes y uniformes sobre un año cíclico de \(N\) días. Sea \(T\) el primer instante en que, entre las personas ya observadas, existen \(k\) cumpleaños contenidos en algún bloque de \(D+1\) días consecutivos del calendario. La meta es calcular exactamente \(\mathbb E[T]\) para el caso terrestre \((N,k,D)=(365,4,7)\), reproduciendo también los casos de control más pequeños.

Enfoque Matemático

Las implementaciones no recurren a simulación aleatoria. Para cada población posible \(n\), cuentan todas las asignaciones de cumpleaños que todavía evitan un conglomerado prohibido de tamaño \(k\), y luego convierten ese conteo en probabilidades de supervivencia.

Paso 1: Reescribir la esperanza como suma de supervivencias

Como \(T\) es un tiempo de parada entero positivo, vale la identidad

$$\mathbb E[T]=\sum_{n\ge 0}\Pr(T>n).$$

Así, el problema pasa a ser: para cada \(n\), ¿cuál es la probabilidad de que, tras ubicar \(n\) personas, todavía no exista un conjunto de \(k\) cumpleaños dentro de una distancia de \(D\) días?

Paso 2: Traducir el evento a una restricción de ventanas cíclicas

Sea \(c_i\) el número de cumpleaños en el día \(i\), interpretando los índices cíclicamente módulo \(N\). Definimos

$$w=D+1.$$

Entonces \(T>n\) equivale a exigir que toda ventana cíclica de longitud \(w\) contenga como máximo \(k-1\) personas:

$$c_i+c_{i+1}+\cdots+c_{i+D}\le k-1 \qquad \text{para todo } i.$$

La formulación también cubre las ventanas que cruzan el final del año, porque el calendario se trata como un ciclo.

Paso 3: Contar vectores de ocupación válidos con pesos exponenciales

Fijemos \(n\). Un vector \((c_1,\dots,c_N)\) con suma total \(n\) describe cuántas de las \(n\) personas etiquetadas caen en cada día. El número de asignaciones que producen exactamente ese vector es

$$\frac{n!}{c_1!\cdots c_N!}.$$

Si \(\mathcal V_n\) es el conjunto de todos esos vectores que además cumplen la restricción de ventanas, entonces

$$\Pr(T>n)=\frac{n!}{N^n}\sum_{(c_1,\dots,c_N)\in \mathcal V_n}\frac{1}{c_1!\cdots c_N!}.$$

Por eso aparece de manera natural la función generadora exponencial: colocar \(c\) personas en un día aporta el peso \(x^c/c!\).

Paso 4: Construir un espacio de estados finito con los últimos \(D\) días

Al avanzar un día, sólo importa lo ocurrido en los \(D\) días anteriores. Por tanto, un estado puede escribirse como

$$s=(a_1,\dots,a_D),\qquad a_j\ge 0,\qquad a_1+\cdots+a_D\le k-1.$$

Ese perfil guarda las ocupaciones de los últimos \(D\) días antes de añadir el nuevo día. Si la suma actual es \(\sigma\), el nuevo día puede recibir cualquier valor \(c\) tal que

$$0\le c\le k-1-\sigma.$$

El estado siguiente se obtiene eliminando la entrada más antigua y añadiendo \(c\):

$$s'=(a_2,\dots,a_D,c).$$

El peso de esta transición en la función generadora es

$$\frac{x^c}{c!}.$$

Paso 5: Imponer el carácter cíclico del año mediante la traza

Si reunimos estas transiciones en una matriz de transferencia \(M(x)\), entonces recorrer \(N\) días corresponde a \(M(x)^N\). Sin embargo, las ventanas que atraviesan el Año Nuevo deben obedecer exactamente la misma regla, así que el perfil final debe coincidir con el inicial. Esa condición diagonal es precisamente la traza:

$$G(x)=\operatorname{tr}(M(x)^N).$$

Si definimos

$$A_n=[x^n]\,G(x),$$

entonces \(A_n\) es el peso exponencial total de todas las configuraciones cíclicas válidas con \(n\) personas, y por tanto

$$\Pr(T>n)=A_n\frac{n!}{N^n}.$$

Paso 6: La suma final es finita

Existe una cota superior limpia para la población antes de que la violación sea inevitable. Si sumamos las \(N\) desigualdades de ventana, obtenemos

$$w(c_1+\cdots+c_N)\le N(k-1),$$

porque cada persona es contada exactamente en \(w\) ventanas. Así, toda configuración válida debe satisfacer

$$n\le \left\lfloor\frac{N(k-1)}{D+1}\right\rfloor.$$

En consecuencia, la esperanza se convierte en una suma finita:

$$\boxed{\mathbb E[T]=\sum_{n=0}^{\left\lfloor N(k-1)/(D+1)\right\rfloor} A_n\frac{n!}{N^n}.}$$

Ejemplo desarrollado: \((N,k,D)=(10,3,1)\)

Aquí \(w=2\), de modo que sobrevivir significa que cada par de días consecutivos contiene como máximo \(2\) personas. Un estado sólo necesita guardar la ocupación del día inmediatamente anterior, así que los estados posibles son \(0,1,2\).

Desde el estado \(0\), el día nuevo puede contener \(0,1\) o \(2\) personas. Desde el estado \(1\), puede contener \(0\) o \(1\). Desde el estado \(2\), necesariamente contiene \(0\). En el orden de estados \((0,1,2)\), la matriz de transferencia es

$$M(x)=\begin{pmatrix} 1 & x & x^2/2 \\ 1 & x & 0 \\ 1 & 0 & 0 \end{pmatrix}.$$

La función generadora cíclica es \(G(x)=\operatorname{tr}(M(x)^{10})\). Al extraer los coeficientes \(A_n\), convertirlos con \(A_n n!/10^n\) y sumar sobre \(n\), se obtiene

$$\mathbb E[T]=5.78688636,$$

que coincide con el valor de control reproducido por las implementaciones.

Cómo Funciona el Código

Las implementaciones en C++, Python y Java siguen exactamente la misma recurrencia. Primero precalculan los recíprocos \(1/c!\) para \(0\le c\le k-1\), ya que ésos son los pesos de transición en la función generadora exponencial. Luego enumeran todos los perfiles válidos de longitud \(D\) cuya suma total no supera \(k-1\), y para cada perfil precalculan todos los perfiles siguientes permitidos junto con la cantidad de personas añadidas.

Después ejecutan una programación dinámica por capas durante exactamente \(N\) días. Para un perfil inicial fijo se guarda, para cada perfil actual y cada población total \(n\), el coeficiente acumulado hasta ese momento. Tras el día \(N\), sólo se lee el vector de coeficientes que regresa al mismo perfil inicial; sumar esas contribuciones diagonales sobre todos los perfiles iniciales equivale exactamente a \(\operatorname{tr}(M(x)^N)\). Finalmente, el coeficiente de grado \(n\) se multiplica por \(n!/N^n\) y se suman todos los términos para obtener \(\mathbb E[T]\). La versión en C++ reparte distintos perfiles iniciales entre varios hilos, pero la matemática es la misma en los tres lenguajes.

Análisis de Complejidad

El número de estados es

$$S=\#\{(a_1,\dots,a_D)\in \mathbb Z_{\ge 0}^D: a_1+\cdots+a_D\le k-1\}=\binom{D+k-1}{D}.$$

Cada estado tiene como máximo \(k\) transiciones salientes, ya que la ocupación del día nuevo sólo puede variar entre \(0\) y \(k-1-\sigma\). El cálculo directo de la traza usado aquí ejecuta la programación dinámica de \(N\) días una vez por cada estado inicial, así que el tiempo total es

$$O\!\left(N\cdot S^2\cdot k\cdot \left\lfloor\frac{N(k-1)}{D+1}\right\rfloor\right).$$

La memoria por trabajador es

$$O\!\left(S\cdot \left\lfloor\frac{N(k-1)}{D+1}\right\rfloor\right),$$

porque sólo se conservan la capa actual y la siguiente de la programación dinámica.

Notas y Referencias

  1. Página del problema: https://projecteuler.net/problem=584
  2. Contexto del problema de cumpleaños: Wikipedia - Birthday problem
  3. Esperanza como suma de supervivencias: Wikipedia - Expected value
  4. Funciones generadoras exponenciales: Wikipedia - Exponential generating function
  5. Método de la matriz de transferencia: Wikipedia - Transfer-matrix method

问题概述

把生日视为在一个长度为 \(N\) 的环形年份上独立且均匀地落点。设 \(T\) 为第一个满足如下条件的到达人数:在已经进入的人中,存在 \(k\) 个生日同时落在某个长度为 \(D+1\) 的连续日历区间内。目标是精确计算地球参数 \((N,k,D)=(365,4,7)\) 时的 \(\mathbb E[T]\),并且与较小的校验样例保持一致。

数学方法

这些实现并不是做随机模拟,而是对每一个可能的人数 \(n\) 统计“仍然没有出现违例簇”的所有生日分布方式,再把这些计数转成生存概率,最后求出期望。

步骤 1:把期望写成生存概率之和

由于 \(T\) 是一个取正整数值的停时,可以使用标准恒等式

$$\mathbb E[T]=\sum_{n\ge 0}\Pr(T>n).$$

因此问题转化为:对于每个 \(n\),在前 \(n\) 个人的生日都放好之后,仍然没有出现“某 \(D+1\) 天窗口内至少有 \(k\) 个生日”这一事件的概率是多少。

步骤 2:把事件改写成环形滑动窗口约束

记 \(c_i\) 为第 \(i\) 天上的生日人数,指标按模 \(N\) 的方式循环读取。设

$$w=D+1.$$

那么事件 \(T>n\) 等价于每个长度为 \(w\) 的环形窗口都至多包含 \(k-1\) 个人:

$$c_i+c_{i+1}+\cdots+c_{i+D}\le k-1 \qquad \text{对所有 } i.$$

这个写法天然包含了跨年窗口,因为年份本身被看作一个圆环,而不是线段。

步骤 3:用指数型权重统计合法占用向量

固定总人数 \(n\)。若 \((c_1,\dots,c_N)\) 的和为 \(n\),它就描述了这 \(n\) 个带标签的人分别落在各个日期上的人数。给定这样一个占用向量,产生它的标号分配总数是

$$\frac{n!}{c_1!\cdots c_N!}.$$

若 \(\mathcal V_n\) 表示所有总和为 \(n\) 且满足窗口约束的占用向量集合,则有

$$\Pr(T>n)=\frac{n!}{N^n}\sum_{(c_1,\dots,c_N)\in \mathcal V_n}\frac{1}{c_1!\cdots c_N!}.$$

这正是指数生成函数出现的原因:如果某一天放入 \(c\) 个人,那么它对应的权重正好是 \(x^c/c!\)。

步骤 4:只保留最近 \(D\) 天即可形成有限状态

当日历从一天推进到下一天时,是否会违规只取决于前面最近的 \(D\) 天,因此状态可以写成

$$s=(a_1,\dots,a_D),\qquad a_j\ge 0,\qquad a_1+\cdots+a_D\le k-1.$$

这个向量记录的是新的一天加入之前,最近 \(D\) 天各自已有多少人。如果当前状态的总和是 \(\sigma\),那么新的一天可以放入的生日人数 \(c\) 必须满足

$$0\le c\le k-1-\sigma.$$

加入这一天以后,新的状态就是把最旧的计数删除,再把 \(c\) 追加到末尾:

$$s'=(a_2,\dots,a_D,c).$$

这一步在生成函数中的权重为

$$\frac{x^c}{c!}.$$

步骤 5:通过迹来强制满足环形年份条件

把这些状态转移放进一个传递矩阵 \(M(x)\) 之后,沿着年份推进 \(N\) 天就对应于 \(M(x)^N\)。但这里不是普通线性时间轴,而是环形年份,因此跨越年末与年初的窗口也必须满足同样的约束。要做到这一点,年终状态必须与年初状态相同,这恰好对应矩阵的对角线贡献,也就是迹:

$$G(x)=\operatorname{tr}(M(x)^N).$$

如果记

$$A_n=[x^n]\,G(x),$$

那么 \(A_n\) 就是所有合法环形生日配置在总人数为 \(n\) 时的指数型总权重,于是

$$\Pr(T>n)=A_n\frac{n!}{N^n}.$$

步骤 6:最终求和实际上是有限的

在没有违反约束的前提下,总人数存在一个直接的上界。把全部 \(N\) 个窗口不等式相加,可得

$$w(c_1+\cdots+c_N)\le N(k-1),$$

因为每一个人恰好会被计算进 \(w\) 个长度为 \(w\) 的窗口。于是任何合法配置都必须满足

$$n\le \left\lfloor\frac{N(k-1)}{D+1}\right\rfloor.$$

因此期望值不是无穷级数,而是一个有限和:

$$\boxed{\mathbb E[T]=\sum_{n=0}^{\left\lfloor N(k-1)/(D+1)\right\rfloor} A_n\frac{n!}{N^n}.}$$

例子:\((N,k,D)=(10,3,1)\)

此时 \(w=2\),所以“仍未触发”意味着任何相邻两天合起来最多只有 \(2\) 个人。状态只需要记录前一天的人数,因此可能状态只有 \(0,1,2\) 三种。

若当前状态为 \(0\),则新的一天可以放 \(0,1,2\) 个人;若当前状态为 \(1\),则只能放 \(0\) 或 \(1\) 人;若当前状态为 \(2\),则新的一天只能放 \(0\) 人。按照状态顺序 \((0,1,2)\),传递矩阵为

$$M(x)=\begin{pmatrix} 1 & x & x^2/2 \\ 1 & x & 0 \\ 1 & 0 & 0 \end{pmatrix}.$$

于是环形生成函数是 \(G(x)=\operatorname{tr}(M(x)^{10})\)。从中提取各个 \(A_n\),再按 \(A_n n!/10^n\) 转成生存概率并求和,就得到

$$\mathbb E[T]=5.78688636,$$

这正是实现所复现的校验值。

代码如何工作

C++、Python 和 Java 实现采用完全相同的递推结构。首先预计算 \(0\le c\le k-1\) 的 \(1/c!\),因为这正是指数生成函数中的转移权重。然后枚举所有长度为 \(D\)、总和不超过 \(k-1\) 的合法状态轮廓,并为每个状态预先列出所有允许的下一状态以及新增人数。

接下来,程序对固定的起始状态执行恰好 \(N\) 层的动态规划。对于每一个当前状态以及每一个总人数 \(n\),它维护到达该位置时累积出来的系数。经过第 \(N\) 天之后,只读取“回到同一起始状态”的那一条系数向量;把所有起始状态的这类对角贡献加起来,正好就是 \(\operatorname{tr}(M(x)^N)\)。最后再把次数为 \(n\) 的系数乘上 \(n!/N^n\),并把所有项相加得到 \(\mathbb E[T]\)。其中 C++ 版本会把不同的起始状态分配给多个线程,但三种语言的数学核心完全一致。

复杂度分析

状态总数为

$$S=\#\{(a_1,\dots,a_D)\in \mathbb Z_{\ge 0}^D: a_1+\cdots+a_D\le k-1\}=\binom{D+k-1}{D}.$$

每个状态的出边数最多为 \(k\),因为新的一天人数只能在 \(0\) 到 \(k-1-\sigma\) 之间变化。当前这种直接求迹的做法,会对每个起始状态各运行一次 \(N\) 天动态规划,因此时间复杂度为

$$O\!\left(N\cdot S^2\cdot k\cdot \left\lfloor\frac{N(k-1)}{D+1}\right\rfloor\right).$$

每个工作单元的空间复杂度为

$$O\!\left(S\cdot \left\lfloor\frac{N(k-1)}{D+1}\right\rfloor\right),$$

因为程序只保留当前层和下一层两份动态规划数组。

注释与参考资料

  1. 题目页面:https://projecteuler.net/problem=584
  2. 生日问题背景:Wikipedia - Birthday problem
  3. 把期望写成生存概率之和:Wikipedia - Expected value
  4. 指数生成函数:Wikipedia - Exponential generating function
  5. 传递矩阵方法:Wikipedia - Transfer-matrix method

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

Дни рождения моделируются как независимые равномерные попадания в циклический год длины \(N\). Обозначим через \(T\) номер человека, после появления которого впервые находятся \(k\) дней рождения, лежащие внутри некоторого блока из \(D+1\) последовательных календарных дней. Требуется точно вычислить \(\mathbb E[T]\) для земного случая \((N,k,D)=(365,4,7)\) и при этом воспроизвести меньшие контрольные примеры.

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

Реализация не опирается на моделирование. Вместо этого для каждого возможного числа людей \(n\) она подсчитывает все размещения дней рождения, в которых ещё нет запрещённого кластера из \(k\) попаданий, а затем переводит эти подсчёты в вероятности выживания.

Шаг 1: Представим ожидание как сумму вероятностей выживания

Так как \(T\) принимает положительные целые значения, справедлива стандартная формула

$$\mathbb E[T]=\sum_{n\ge 0}\Pr(T>n).$$

Значит, задача сводится к следующему: какова для каждого \(n\) вероятность того, что после размещения \(n\) человек ещё не возникло \(k\) дней рождения внутри промежутка длины \(D\) дней?

Шаг 2: Заменим событие циклическим оконным ограничением

Пусть \(c_i\) обозначает число дней рождения в день \(i\), причём индексы понимаются по модулю \(N\). Введём

$$w=D+1.$$

Тогда условие \(T>n\) эквивалентно тому, что каждое циклическое окно длины \(w\) содержит не более \(k-1\) человек:

$$c_i+c_{i+1}+\cdots+c_{i+D}\le k-1 \qquad \text{для всех } i.$$

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

Шаг 3: Считаем допустимые векторы заполнения с экспоненциальными весами

Зафиксируем \(n\). Вектор \((c_1,\dots,c_N)\) с суммой \(n\) описывает, сколько из \(n\) размеченных людей попало на каждый день. Число назначений, порождающих именно этот вектор, равно

$$\frac{n!}{c_1!\cdots c_N!}.$$

Если через \(\mathcal V_n\) обозначить множество всех таких векторов, которые одновременно удовлетворяют оконному ограничению, то

$$\Pr(T>n)=\frac{n!}{N^n}\sum_{(c_1,\dots,c_N)\in \mathcal V_n}\frac{1}{c_1!\cdots c_N!}.$$

Именно поэтому естественно возникает экспоненциальная производящая функция: размещение \(c\) человек в один день даёт вес \(x^c/c!\).

Шаг 4: Вводим конечное пространство состояний по последним \(D\) дням

При переходе к следующему дню важны только заполнения последних \(D\) дней. Поэтому состояние можно записать как

$$s=(a_1,\dots,a_D),\qquad a_j\ge 0,\qquad a_1+\cdots+a_D\le k-1.$$

Этот профиль хранит количества людей в последних \(D\) днях перед добавлением нового дня. Если сумма в текущем состоянии равна \(\sigma\), то в новый день можно поставить любое число \(c\), удовлетворяющее

$$0\le c\le k-1-\sigma.$$

Следующее состояние получается сдвигом и дописыванием \(c\):

$$s'=(a_2,\dots,a_D,c).$$

Вес такого перехода в производящей функции равен

$$\frac{x^c}{c!}.$$

Шаг 5: Цикличность года задаётся следом матрицы

Если собрать переходы в матрицу переноса \(M(x)\), то проход по \(N\) дням соответствует \(M(x)^N\). Но окна, проходящие через конец и начало года, должны удовлетворять тем же ограничениям, что и внутренние окна. Поэтому конечный профиль обязан совпадать с начальным; это как раз диагональное условие, то есть след:

$$G(x)=\operatorname{tr}(M(x)^N).$$

Если

$$A_n=[x^n]\,G(x),$$

то \(A_n\) есть полный экспоненциальный вес всех допустимых циклических конфигураций с \(n\) людьми, а значит

$$\Pr(T>n)=A_n\frac{n!}{N^n}.$$

Шаг 6: Итоговая сумма конечна

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

$$w(c_1+\cdots+c_N)\le N(k-1),$$

потому что каждый человек учитывается ровно в \(w\) окнах. Следовательно, любая допустимая конфигурация обязана удовлетворять

$$n\le \left\lfloor\frac{N(k-1)}{D+1}\right\rfloor.$$

Значит, ожидание представляется конечной суммой:

$$\boxed{\mathbb E[T]=\sum_{n=0}^{\left\lfloor N(k-1)/(D+1)\right\rfloor} A_n\frac{n!}{N^n}.}$$

Разобранный пример: \((N,k,D)=(10,3,1)\)

Здесь \(w=2\), поэтому условие выживания означает, что в любой паре соседних дней суммарно не более \(2\) человек. Состояние хранит только число людей в предыдущий день, так что возможны состояния \(0,1,2\).

Из состояния \(0\) в новый день можно добавить \(0,1\) или \(2\) человека. Из состояния \(1\) разрешены \(0\) или \(1\). Из состояния \(2\) разрешён только \(0\). В порядке состояний \((0,1,2)\) матрица переноса имеет вид

$$M(x)=\begin{pmatrix} 1 & x & x^2/2 \\ 1 & x & 0 \\ 1 & 0 & 0 \end{pmatrix}.$$

Циклическая производящая функция равна \(G(x)=\operatorname{tr}(M(x)^{10})\). После извлечения коэффициентов \(A_n\), преобразования по формуле \(A_n n!/10^n\) и суммирования по \(n\) получаем

$$\mathbb E[T]=5.78688636,$$

что совпадает с контрольным значением из реализаций.

Как работает код

Реализации на C++, Python и Java используют одну и ту же рекуррентную схему. Сначала они заранее вычисляют величины \(1/c!\) для \(0\le c\le k-1\), потому что именно эти числа являются весами переходов в экспоненциальной производящей функции. Затем перечисляются все допустимые профили длины \(D\) с суммой не более \(k-1\), и для каждого профиля заранее строится список разрешённых следующих профилей вместе с числом добавленных людей.

После этого запускается послойная динамика ровно на \(N\) дней. Для фиксированного начального профиля хранится коэффициент для каждой текущей вершины состояний и каждого общего числа людей \(n\). После \(N\)-го дня считывается только тот вектор коэффициентов, который вернулся в тот же начальный профиль; сумма всех таких диагональных вкладов по всем начальными профилям и есть \(\operatorname{tr}(M(x)^N)\). На последнем шаге коэффициент степени \(n\) умножается на \(n!/N^n\), и все слагаемые складываются в \(\mathbb E[T]\). В версии на C++ разные начальные профили дополнительно распределяются по потокам, но математически три реализации эквивалентны.

Анализ сложности

Число состояний равно

$$S=\#\{(a_1,\dots,a_D)\in \mathbb Z_{\ge 0}^D: a_1+\cdots+a_D\le k-1\}=\binom{D+k-1}{D}.$$

У каждого состояния не более \(k\) исходящих переходов, поскольку число людей в новом дне может принимать значения только от \(0\) до \(k-1-\sigma\). Прямое вычисление следа, используемое здесь, запускает \(N\)-дневную динамику отдельно для каждого начального состояния, поэтому временная сложность составляет

$$O\!\left(N\cdot S^2\cdot k\cdot \left\lfloor\frac{N(k-1)}{D+1}\right\rfloor\right).$$

Память на один рабочий процесс равна

$$O\!\left(S\cdot \left\lfloor\frac{N(k-1)}{D+1}\right\rfloor\right),$$

так как хранятся только текущий и следующий слои динамики.

Сноски и ссылки

  1. Страница задачи: https://projecteuler.net/problem=584
  2. Фон: задача о днях рождения: Wikipedia - Birthday problem
  3. Ожидание как сумма вероятностей выживания: Wikipedia - Expected value
  4. Экспоненциальные производящие функции: Wikipedia - Exponential generating function
  5. Метод матрицы переноса: Wikipedia - Transfer-matrix method

ملخص المسألة

تُنمذج أعياد الميلاد على أنها اختيارات مستقلة ومتساوية الاحتمال على سنة دائرية طولها \(N\) يومًا. نرمز بـ \(T\) إلى أول عدد من الداخلين تصبح عنده هناك \(k\) أعياد ميلاد تقع كلها داخل كتلة من \(D+1\) أيام تقويمية متتالية. المطلوب هو حساب \(\mathbb E[T]\) بدقة للحالة الأرضية \((N,k,D)=(365,4,7)\)، مع مطابقة حالات التحقق الأصغر أيضًا.

المنهج الرياضي

لا تعتمد الشيفرات على المحاكاة العشوائية. بدلًا من ذلك فهي تعد، لكل عدد محتمل من الأشخاص \(n\)، جميع توزيعات أعياد الميلاد التي ما زالت تتجنب ظهور عنقود محظور حجمه \(k\)، ثم تحوّل هذا العد إلى احتمالات بقاء، ومنها إلى القيمة المتوقعة.

الخطوة 1: كتابة التوقع على صورة مجموع احتمالات البقاء

بما أن \(T\) زمن توقف يأخذ قيمًا صحيحة موجبة، فلدينا الهوية المعروفة

$$\mathbb E[T]=\sum_{n\ge 0}\Pr(T>n).$$

إذًا يكفي أن نحسب، لكل \(n\)، احتمال أن أول \(n\) أشخاص لم يكوّنوا بعد مجموعة من \(k\) أعياد ميلاد داخل مسافة زمنية قدرها \(D\) أيام.

الخطوة 2: تحويل الحدث إلى قيد على نوافذ دائرية

ليكن \(c_i\) عدد أعياد الميلاد في اليوم \(i\)، مع قراءة الفهارس دوريًا بترديد \(N\). نضع

$$w=D+1.$$

عندئذ يكون الشرط \(T>n\) مكافئًا لكون كل نافذة دائرية طولها \(w\) تحتوي على الأكثر \(k-1\) شخصًا:

$$c_i+c_{i+1}+\cdots+c_{i+D}\le k-1 \qquad \text{for all } i.$$

وهذا يضم تلقائيًا النوافذ التي تعبر نهاية السنة، لأن السنة تُعامل كدائرة لا كقطعة مستقيمة.

الخطوة 3: عد متجهات الإشغال الصحيحة بأوزان أسية

ثبّت \(n\). المتجه \((c_1,\dots,c_N)\) الذي مجموعه \(n\) يصف كيف توزّع \(n\) شخصًا مميّزين على الأيام المختلفة. وعدد التعيينات التي تعطي هذا المتجه هو

$$\frac{n!}{c_1!\cdots c_N!}.$$

إذا كانت \(\mathcal V_n\) هي مجموعة جميع هذه المتجهات التي تحقق أيضًا قيد النوافذ، فإن

$$\Pr(T>n)=\frac{n!}{N^n}\sum_{(c_1,\dots,c_N)\in \mathcal V_n}\frac{1}{c_1!\cdots c_N!}.$$

ولهذا تظهر الدوال المولدة الأسية بصورة طبيعية: وضع \(c\) أشخاص في يوم واحد يعطي الوزن \(x^c/c!\).

الخطوة 4: بناء فضاء حالات منتهٍ اعتمادًا على آخر \(D\) أيام

عند الانتقال من يوم إلى اليوم التالي، فإن المعلومات المهمة الوحيدة هي أعداد الأشخاص في آخر \(D\) أيام. لذلك يمكن كتابة الحالة على الصورة

$$s=(a_1,\dots,a_D),\qquad a_j\ge 0,\qquad a_1+\cdots+a_D\le k-1.$$

هذا الملف العددي يحفظ الإشغال في آخر \(D\) أيام قبل إضافة اليوم الجديد. إذا كان مجموع الحالة الحالية هو \(\sigma\)، فإن اليوم الجديد يمكن أن يستقبل أي عدد \(c\) يحقق

$$0\le c\le k-1-\sigma.$$

وتُبنى الحالة التالية بحذف أقدم قيمة ثم إلحاق \(c\) في النهاية:

$$s'=(a_2,\dots,a_D,c).$$

أما وزن هذا الانتقال في الدالة المولدة فهو

$$\frac{x^c}{c!}.$$

الخطوة 5: فرض دورية السنة باستعمال الأثر

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

$$G(x)=\operatorname{tr}(M(x)^N).$$

وإذا كتبنا

$$A_n=[x^n]\,G(x),$$

فإن \(A_n\) يمثل مجموع الأوزان الأسية لكل التوزيعات الدائرية الصحيحة التي تضم \(n\) شخصًا، ومن ثم

$$\Pr(T>n)=A_n\frac{n!}{N^n}.$$

الخطوة 6: المجموع النهائي منتهٍ

هناك حد أعلى مباشر لعدد الأشخاص الممكن قبل أن يصبح الخرق حتميًا. إذا جمعنا متراجحات النوافذ كلها وعددها \(N\)، نحصل على

$$w(c_1+\cdots+c_N)\le N(k-1),$$

لأن كل شخص يُحسب بالضبط في \(w\) نوافذ. لذلك يجب أن يحقق كل توزيع صحيح

$$n\le \left\lfloor\frac{N(k-1)}{D+1}\right\rfloor.$$

وعليه يصبح التوقع مجموعًا منتهيًا:

$$\boxed{\mathbb E[T]=\sum_{n=0}^{\left\lfloor N(k-1)/(D+1)\right\rfloor} A_n\frac{n!}{N^n}.}$$

مثال محلول: \((N,k,D)=(10,3,1)\)

هنا \(w=2\)، لذلك فإن شرط البقاء يعني أن كل زوج من الأيام المتجاورة يحتوي على الأكثر \(2\) شخصًا. الحالة تحتاج فقط إلى حفظ عدد الأشخاص في اليوم السابق مباشرة، ومن ثم فالحالات الممكنة هي \(0,1,2\).

من الحالة \(0\) يمكن أن يحتوي اليوم الجديد على \(0\) أو \(1\) أو \(2\) شخص. ومن الحالة \(1\) يمكنه أن يحتوي على \(0\) أو \(1\) فقط. ومن الحالة \(2\) لا بد أن يحتوي اليوم الجديد على \(0\). إذا أخذنا ترتيب الحالات \((0,1,2)\)، فإن مصفوفة الانتقال هي

$$M(x)=\begin{pmatrix} 1 & x & x^2/2 \\ 1 & x & 0 \\ 1 & 0 & 0 \end{pmatrix}.$$

ومن ثم تكون الدالة المولدة الدائرية هي \(G(x)=\operatorname{tr}(M(x)^{10})\). وبعد استخراج معاملات \(A_n\)، وتحويلها بواسطة \(A_n n!/10^n\)، ثم جمعها على \(n\)، نحصل على

$$\mathbb E[T]=5.78688636,$$

وهو تمامًا مقدار التحقق الذي تعيده الشيفرات.

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

تتبع تطبيقات C++ وPython وJava العلاقة نفسها. في البداية تُحسب القيم \(1/c!\) مسبقًا لكل \(0\le c\le k-1\)، لأن هذه هي أوزان الانتقال في الدالة المولدة الأسية. بعد ذلك تُحصى كل الملفات الحالة ذات الطول \(D\) والتي لا يتجاوز مجموعها \(k-1\)، ويُبنى لكل حالة جدول مسبق بالحالات التالية المسموح بها مع عدد الأشخاص المضافين في كل انتقال.

ثم تُنفذ برمجة ديناميكية طبقية لمدة \(N\) يومًا بالضبط. بالنسبة إلى حالة ابتدائية ثابتة، تحتفظ الشيفرة بمعامل لكل حالة حالية ولكل عدد كلي من الأشخاص \(n\). وبعد اليوم رقم \(N\)، لا يُقرأ إلا متجه المعاملات الذي عاد إلى الحالة الابتدائية نفسها؛ وجمع هذه المساهمات القطرية على جميع الحالات الابتدائية يعطينا بالضبط \(\operatorname{tr}(M(x)^N)\). وفي النهاية يُضرب معامل الدرجة \(n\) في \(n!/N^n\) ثم تُجمع الحدود كلها للحصول على \(\mathbb E[T]\). إصدار C++ يوزع الحالات الابتدائية على عدة خيوط، لكن البنية الرياضية نفسها مشتركة بين اللغات الثلاث.

تحليل التعقيد

عدد الحالات يساوي

$$S=\#\{(a_1,\dots,a_D)\in \mathbb Z_{\ge 0}^D: a_1+\cdots+a_D\le k-1\}=\binom{D+k-1}{D}.$$

ولكل حالة في أقصى تقدير \(k\) انتقالات خارجة، لأن عدد الأشخاص في اليوم الجديد لا يمكن أن يتجاوز \(k-1-\sigma\). والحساب المباشر للأثر كما هو مستخدم هنا يشغّل البرمجة الديناميكية ذات \(N\) يومًا مرة واحدة لكل حالة ابتدائية، لذا يكون الزمن

$$O\!\left(N\cdot S^2\cdot k\cdot \left\lfloor\frac{N(k-1)}{D+1}\right\rfloor\right).$$

أما الذاكرة لكل عامل فهي

$$O\!\left(S\cdot \left\lfloor\frac{N(k-1)}{D+1}\right\rfloor\right),$$

لأن المطلوب حفظ الطبقة الحالية والطبقة التالية فقط.

هوامش ومراجع

  1. صفحة المسألة: https://projecteuler.net/problem=584
  2. خلفية مسألة أعياد الميلاد: Wikipedia - Birthday problem
  3. التوقع بوصفه مجموع احتمالات البقاء: Wikipedia - Expected value
  4. الدوال المولدة الأسية: Wikipedia - Exponential generating function
  5. طريقة مصفوفة الانتقال: Wikipedia - Transfer-matrix method