Detailed mathematical approach

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