Detailed mathematical approach

Problem Summary

A uniformly shuffled standard deck has \(52\) cards, arranged as \(13\) ranks with \(4\) cards of each rank. A rank is called perfect if no two cards of that rank appear in adjacent positions in the shuffled deck. If \(Y\) denotes the number of perfect ranks, the problem asks for

$\Pr\!\left(Y\in\{2,3,5,7,11,13\}\right).$

Brute force over all \(52!\) permutations is hopeless, so the implementations reorganize the question into two layers: first compute the probability that every rank in a fixed subset is perfect, then recover the full distribution of \(Y\) from those subset probabilities.

Mathematical Approach

The key observation is that the exact identities of the tracked ranks do not matter; only how many ranks are being tracked matters. That symmetry makes it possible to solve thirteen smaller probability problems and then combine them.

Step 1: Fix a subset of ranks

Choose a subset \(S\) of \(j\) ranks, where \(0\le j\le 13\), and define

$q_j=\Pr(\text{every rank in }S\text{ is perfect}).$

Because every rank has the same multiplicity and every shuffle is uniform, \(q_j\) depends only on \(j\), not on the actual labels of the ranks in \(S\).

If we can evaluate \(q_j\) for all \(j\), then we can rebuild the distribution of the random variable \(Y\).

Step 2: Turn subset probabilities into binomial moments

For each fixed \(j\)-subset \(S\), let \(I_S\) be the indicator that all ranks in \(S\) are perfect. In any particular shuffle with exactly \(Y\) perfect ranks, the number of \(j\)-subsets that satisfy this is \(\binom{Y}{j}\). Therefore

$\sum_{\lvert S\rvert=j} I_S=\binom{Y}{j}.$

Taking expectations and using symmetry gives

$B_j=\binom{13}{j}q_j=\mathbb{E}\!\left[\binom{Y}{j}\right].$

Now write

$D_k=\Pr(Y=k),\qquad 0\le k\le 13.$

Then the \(j\)-th binomial moment satisfies

$B_j=\sum_{k=j}^{13}\binom{k}{j}D_k.$

This is an upper-triangular linear system, so once the \(B_j\) values are known, the exact probabilities \(D_k\) can be recovered by backward substitution.

Step 3: Reveal the deck sequentially and compress the state

To compute \(q_j\), imagine revealing the deck from left to right. Cards whose ranks are outside \(S\) are irrelevant except that they can separate two tracked cards and thereby prevent an adjacency violation.

At any point, the state can be summarized by:

$u=\text{number of remaining cards from untracked ranks},$

$a_r=\text{number of tracked ranks with exactly }r\text{ unseen cards left},\qquad r=1,2,3,4,$

and

$s=\text{remaining multiplicity of the most recently drawn tracked rank},$

with the convention \(s=0\) if the previous card was untracked or if the previous tracked rank has already been exhausted.

This compressed description is sufficient because ranks in the same class are exchangeable. The only rank that needs to be distinguished is the most recently drawn tracked rank, since drawing it again immediately would create an adjacent equal-rank pair and destroy perfection for that rank.

The initial state for a subset of size \(j\) is

$\Psi_j(52-4j,0,0,0,j,0),$

where \(\Psi_j\) denotes the probability of eventually finishing the reveal with every tracked rank still perfect.

Step 4: Write the transition probabilities

From a state \((u,a_1,a_2,a_3,a_4,s)\), the number of cards still unseen is

$R=u+a_1+2a_2+3a_3+4a_4.$

If \(R=0\), the tracked subset has survived the whole shuffle without any forbidden adjacency, so

$\Psi_j(0,0,0,0,0,0)=1.$

If an untracked card is drawn, this happens with probability \(u/R\), the count \(u\) drops by \(1\), and the active restriction disappears because an untracked card separates future tracked cards from the previously drawn tracked rank.

If a tracked rank currently belongs to class \(r\), then each such rank contributes \(r\) possible next cards. However, if \(s=r\), one of those ranks is the forbidden rank that cannot be repeated immediately, so the number of eligible ranks in class \(r\) is

$a_r-\mathbf{1}_{\{s=r\}}.$

Hence the probability of drawing a tracked card from class \(r\) is

$\frac{r\left(a_r-\mathbf{1}_{\{s=r\}}\right)}{R}.$

After such a draw, one tracked rank moves from class \(r\) to class \(r-1\). The new active value becomes \(r-1\), except that when \(r=1\) the rank is exhausted and the active value resets to \(0\).

Each state has at most five outgoing transitions: one for an untracked card and one for each multiplicity class \(r=1,2,3,4\). Memoization makes this recurrence practical because the same compressed state can be reached through many different reveal histories.

Step 5: Recover the exact distribution of \(Y\)

Once \(q_j=\Psi_j(52-4j,0,0,0,j,0)\) is known for every \(j\), we compute

$B_j=\binom{13}{j}q_j.$

The triangular relation

$B_j=\sum_{k=j}^{13}\binom{k}{j}D_k$

can then be inverted from \(k=13\) downward:

$D_k=B_k-\sum_{t=k+1}^{13}\binom{t}{k}D_t.$

Finally, the required probability is

$\Pr\!\left(Y\in\{2,3,5,7,11,13\}\right)=D_2+D_3+D_5+D_7+D_{11}+D_{13}.$

Worked Example: A single fixed rank

For \(j=1\), let us track one specific rank. It is perfect exactly when its four positions in the shuffled deck contain no adjacent pair.

If those positions are

$1\le x_1<x_2<x_3<x_4\le 52,\qquad x_{i+1}\ge x_i+2,$

then setting \(y_i=x_i-(i-1)\) transforms them into

$1\le y_1<y_2<y_3<y_4\le 49.$

So the number of non-adjacent \(4\)-position choices is \(\binom{49}{4}\), while the total number of \(4\)-position choices is \(\binom{52}{4}\). Therefore

$q_1=\frac{\binom{49}{4}}{\binom{52}{4}}=\frac{4324}{5525}.$

By linearity of expectation,

$\mathbb{E}[Y]=13q_1=\frac{4324}{425}.$

This matches the consistency check built into the implementations and confirms that the subset-probability viewpoint is aligned with the combinatorics of the shuffle.

How the Code Works

The C++, Python, and Java implementations all follow the same sequence. First they precompute the binomial coefficients \(\binom{n}{k}\) for \(0\le n,k\le 13\). Then, for each subset size \(j\), they evaluate the memoized probability recurrence starting from the initial state with \(j\) tracked ranks and \(52-4j\) untracked cards.

Those thirteen probabilities produce the values \(q_1,\dots,q_{13}\), and from them the implementations form the binomial moments \(B_j=\binom{13}{j}q_j\). The triangular system is then solved backwards to obtain the full distribution \(D_0,D_1,\dots,D_{13}\).

After the distribution is known, the implementation adds the six probabilities corresponding to prime values of \(Y\) and prints the result to ten decimal places. One implementation also checks the identity \(\mathbb{E}[Y]=4324/425\), which is exactly the worked example above written as a sanity test.

Complexity Analysis

For a fixed subset size \(j\), let \(S_j\) be the number of reachable compressed states \((u,a_1,a_2,a_3,a_4,s)\). Memoization ensures that each such state is evaluated once, and each evaluation considers at most five transitions. Therefore the running time for that \(j\) is \(O(S_j)\) and the memory usage is also \(O(S_j)\).

Across all subset sizes, the total cost is

$O\!\left(\sum_{j=1}^{13} S_j\right)+O(13^2),$

where the \(O(13^2)\) term is the final backward substitution. There is no simple closed form for \(S_j\), but the state space is tightly constrained by

$u+a_1+2a_2+3a_3+4a_4\le 52$

and by the fact that only \(13\) ranks exist. In practice the reachable-state count is small enough that the exact probability can be computed comfortably.

Footnotes and References

1. Project Euler problem page: https://projecteuler.net/problem=687
2. Binomial coefficient: Wikipedia - Binomial coefficient
3. Linearity of expectation: Wikipedia - Expected value
4. Dynamic programming: Wikipedia - Dynamic programming
5. Binomial transform and inversion: Wikipedia - Binomial transform