Solution

Answer: 0.3285320869

Let $X$ be the number of perfect ranks in a uniformly shuffled deck.

We work with the $52!$ permutations of the distinct cards.

A rank is perfect if no two cards of that rank are adjacent.

1. Mathematical analysis

For each rank, there are $4$ cards.

Define:

$$A_r = {\text{rank } r \text{ is perfect}}.$$

We need

$$\Pr(X \text{ is prime}).$$

The possible prime values are

$$2,3,5,7,11,13.$$

Step 1: Probability that a fixed set of ranks are perfect

Suppose a fixed set of $t$ ranks must all be perfect.

We count permutations satisfying these conditions by inclusion–exclusion.

Inclusion–exclusion for one rank

Take one rank with cards

$$a_1,a_2,a_3,a_4.$$

To force adjacency, we merge adjacent cards into blocks.

If exactly $j$ adjacencies are enforced ($j=0,1,2,3$), then the four cards become

$$4-j$$

blocks.

The number of ways to choose those adjacencies is

$$\binom{3}{j}.$$

After merging, the total number of objects becomes

$$52-j.$$

Hence the contribution is

$$(-1)^j \binom{3}{j}\frac{(52-j)!}{(4-j)!}.$$

For $t$ specified ranks

Let rank $i$ contribute $j_i\in{0,1,2,3}$.

Then:

total merges:

$$s=j_1+\cdots+j_t,$$

total objects:

$$52-s,$$

each constrained rank contributes a factor

$$(-1)^{j_i}\binom{3}{j_i}\frac1{(4-j_i)!}.$$

The remaining $13-t$ ranks contribute $(4!)^{-(13-t)}$.

Thus the number of favorable permutations is

$$N_t = \sum_{j_1,\dots,j_t} (52-s)! \prod_{i=1}^t \left( (-1)^{j_i} \binom{3}{j_i} \frac1{(4-j_i)!} \right) \frac1{(4!)^{13-t}}.$$

The total number of rank-pattern permutations is

$$\frac{52!}{(4!)^{13}}.$$

Therefore

$$p_t = \Pr(\text{a fixed set of }t\text{ ranks are perfect}) = \frac{N_t}{52!/(4!)^{13}}.$$

Step 2: Recover the distribution of $X$

By symmetry,

$$\binom{13}{t}p_t = \mathbb E!\left[\binom Xt\right].$$

Hence

$$\binom{13}{t}p_t = \sum_{k=t}^{13}\binom{k}{t}\Pr(X=k).$$

This is a binomial transform. Inverting it gives

$$\Pr(X=k) = \sum_{t=k}^{13} (-1)^{t-k} \binom{t}{k} \binom{13}{t} p_t.$$

Then we sum over prime $k$.

2. Python implementation

from math import factorial, comb
from collections import defaultdict

# Total number of rank-pattern permutations
TOTAL = factorial(52) / (factorial(4) ** 13)

# p[t] = probability that a fixed set of t ranks are perfect
p = []

for t in range(14):

    # DP over total number of merges s
    dp = {0: 1.0}

    for _ in range(t):
        ndp = defaultdict(float)

        for s, val in dp.items():
            for j in range(4):
                ndp[s + j] += (
                    val
                    * ((-1) ** j)
                    * comb(3, j)
                    / factorial(4 - j)
                )

        dp = ndp

    favorable = 0.0

    for s, coeff in dp.items():
        favorable += (
            factorial(52 - s)
            * coeff
            / (factorial(4) ** (13 - t))
        )

    p.append(favorable / TOTAL)

# Recover exact distribution of X
P = [0.0] * 14

for k in range(14):
    total = 0.0

    for t in range(k, 14):
        total += (
            ((-1) ** (t - k))
            * comb(t, k)
            * comb(13, t)
            * p[t]
        )

    P[k] = total

# Prime values
primes = [2, 3, 5, 7, 11, 13]

answer = sum(P[k] for k in primes)

print(f"{answer:.10f}")

3. Code walkthrough

Computing $p_t$

The dictionary dp[s] stores the total inclusion–exclusion coefficient for exactly $s$ merges.

For each constrained rank:

for j in range(4):

corresponds to enforcing $j$ adjacencies.

The factor

((-1) ** j) * comb(3, j) / factorial(4 - j)

is exactly the inclusion–exclusion contribution.

After processing all $t$ ranks, each state s contributes

factorial(52 - s)

because the merged configuration has $52-s$ objects.

Recovering the distribution

We use the inversion formula

$$\Pr(X=k) = \sum_{t=k}^{13} (-1)^{t-k} \binom{t}{k} \binom{13}{t} p_t.$$

This produces the exact probability for every possible value of $X$.

Finally we sum the prime cases.

4. Final verification

We are told

$$\mathbb E[X] = \frac{4324}{425} \approx 10.1741,$$

so values around $10$ are most likely.

The resulting probability

$$0.3285\ldots$$

is reasonable because among the nearby values only $11$ is prime.

Answer: 0.3285320868

Project Euler Problem 687