Detailed mathematical approach

Problem Summary

There are \(n\) distinct values, and each value appears exactly four times in a shuffled deck of \(4n\) cards. As cards are revealed one by one, a value-group is called open if it has appeared at least once but fewer than four times. If \(M\) denotes the largest number of open groups seen at any moment, the goal is to compute \(\mathbb{E}[M]\) exactly for \(n=60\), then round the result to eight decimal places.

A direct simulation would only produce an approximation. The implementation instead performs an exact probability dynamic program over aggregated deck states, using symmetry to ignore the actual labels and track only how many values are in each partial-completion stage.

Mathematical Approach

The essential simplification is exchangeability: once we know how many values have been seen \(0\), \(1\), \(2\), or \(3\) times, the precise identities of those values no longer matter for future probabilities.

Step 1: Compress the state by revealed multiplicities

Let \(c_0,c_1,c_2,c_3\) be the numbers of values for which exactly \(0,1,2,3\) copies have already been revealed. Values with four revealed copies are closed, so they do not need their own tracked bucket.

The number of currently open groups is therefore

$p=c_1+c_2+c_3.$

The number of already closed values is

$c_4=n-(c_0+c_1+c_2+c_3).$

This compressed description is exact. Two partial decks with the same quadruple \((c_0,c_1,c_2,c_3)\) have the same number of hidden cards of every type, hence the same transition probabilities.

Step 2: Derive the four transition probabilities

From state \((c_0,c_1,c_2,c_3)\), the number of cards still hidden is

$r=4c_0+3c_1+2c_2+c_3.$

Exactly four kinds of next draw can occur:

$\begin{aligned} (c_0,c_1,c_2,c_3)&\to(c_0-1,c_1+1,c_2,c_3) &&\text{with probability } \frac{4c_0}{r},\\ (c_0,c_1,c_2,c_3)&\to(c_0,c_1-1,c_2+1,c_3) &&\text{with probability } \frac{3c_1}{r},\\ (c_0,c_1,c_2,c_3)&\to(c_0,c_1,c_2-1,c_3+1) &&\text{with probability } \frac{2c_2}{r},\\ (c_0,c_1,c_2,c_3)&\to(c_0,c_1,c_2,c_3-1) &&\text{with probability } \frac{c_3}{r}. \end{aligned}$

The first copy of a value opens a new group, the second and third copies keep that group open, and the fourth copy closes it.

Step 3: Replace the maximum by a threshold event

For a fixed integer threshold \(t\), define

$Q_t=\mathbb{P}(M<t).$

Instead of computing the full distribution of \(M\) directly, the dynamic program keeps only those paths for which the open-group count always remains strictly below \(t\). If a transition would produce a new open count \(p'\ge t\), that transition is discarded immediately.

After all \(4n\) cards have been processed, the surviving probability mass is exactly \(Q_t\).

Step 4: Recover the expectation from tail probabilities

Since \(M\) is an integer between \(1\) and \(n\), the tail-sum identity gives

$\mathbb{E}[M]=\sum_{t=1}^{n}\mathbb{P}(M\ge t)=\sum_{t=1}^{n}\left(1-Q_t\right).$

So the exact answer is obtained by running the thresholded DP for every \(t=1,2,\dots,n\), then summing the complements \(1-Q_t\).

Step 5: Why the compression is sufficient

No future probability depends on whether a particular value is, say, the number \(7\) or the number \(43\). What matters is only how many values contribute four unseen cards, three unseen cards, two unseen cards, or one unseen card. That symmetry turns the original deck process into a finite Markov process on aggregated count states.

Worked Example: \(n=2\)

With only two values, the maximum \(M\) can be either \(1\) or \(2\). Therefore

$\mathbb{E}[M]=1+\mathbb{P}(M=2)=2-\mathbb{P}(M<2).$

To keep \(M<2\), the second value must not appear before the first value has already been completed. Starting from \((2,0,0,0)\), the only safe chain of states is

$ (2,0,0,0)\to(1,1,0,0)\to(1,0,1,0)\to(1,0,0,1)\to(1,0,0,0), $

because opening the second value any earlier would create two open groups. The corresponding probability is

$ 1\cdot\frac{3}{7}\cdot\frac{2}{6}\cdot\frac{1}{5}=\frac{1}{35}. $

Hence

$ \mathbb{P}(M<2)=\frac{1}{35},\qquad \mathbb{E}[M]=2-\frac{1}{35}=\frac{69}{35}=1.97142857\ldots $

This is the same checkpoint value used by the implementations before the full \(n=60\) computation.

How the Code Works

The C++, Python, and Java implementations all follow the same structure. For one threshold \(t\), they keep a probability map for the current draw depth and another for the next depth. Each compressed state is encoded into a compact integer key so that hash-based containers can merge probability contributions efficiently.

At each step, the implementation decodes one state, computes the remaining-card count \(r\), applies the four transitions above, and drops every successor whose open-group count would reach or exceed the threshold. After \(4n\) layers, the probability attached to the terminal empty state is exactly \(Q_t=\mathbb{P}(M<t)\). An outer loop repeats this for all thresholds from \(1\) to \(n\), accumulates \(1-Q_t\), and finally rounds the result to eight decimal places.

Complexity Analysis

For a fixed threshold \(t\), the DP has \(4n\) layers and stores only states with fewer than \(t\) open groups. At a fixed layer, the draw depth determines one linear relation among the state counts, so the remaining feasible configurations are described by \(O(t^3)\) possibilities. This yields \(O(nt^3)\) time for one threshold and \(O(t^3)\) memory for the two active layers. Summing over \(t=1,\dots,n\) gives an overall worst-case bound of \(O(n^5)\) time and \(O(n^3)\) memory. In practice the active frontier is much smaller, because unreachable states never enter the hash maps and the threshold pruning is strong.

Footnotes and References

1. Project Euler problem page: https://projecteuler.net/problem=815
2. Markov chain overview: Wikipedia - Markov chain
3. Expected value and tail-sum identities: Wikipedia - Expected value
4. Dynamic programming overview: Wikipedia - Dynamic programming
5. Stars and bars for counting constrained tuples: Wikipedia - Stars and bars