Project Euler Problem 595

A deck of cards numbered from 1 to n is shuffled randomly such that each permutation is equally likely.

Project Euler Problem 595

Solution

Answer: 54.17529329

Let $S(n)$ denote the expected number of shuffles needed when starting from a uniformly random permutation of $n$ cards.

1. Mathematical analysis

Step 1: Reformulating the process

At any stage, cards that are already in consecutive ascending order are attached into bundles.

Example:

$$4123756 \to (4)(123)(7)(56)$$

The crucial observation:

Once cards are bundled, only the relative order of bundles matters.

If there are $k$ bundles, then after throwing them in the air and picking them up, we obtain a uniformly random permutation of the $k$ bundles.

Thus the state depends only on:

$$k = \text{number of bundles}$$

Let:

$$E_k = \text{expected additional shuffles needed from a random arrangement of } k \text{ bundles}.$$

We want:

$$S(n)=E_n.$$


Step 2: Transition probabilities

Suppose we currently have $k$ bundles.

A bundle boundary disappears when two consecutive bundles appear adjacent in the correct order.

Define:

$$r = \text{number of successful adjacencies}$$

Then after inspection:

$$j = k-r$$

bundles remain.

We need the probability $p_{k,j}$ that a random permutation of $k$ objects produces exactly $j$ bundles.

Equivalently, we count permutations with exactly $r=k-j$ adjacent ordered pairs:

$$(i,i+1)$$

for $i=1,\dots,k-1$.

Using inclusion–exclusion:

$$N_{k,r} = \sum_{s=r}^{k-1} (-1)^{s-r} \binom{s}{r} \binom{k-1}{s} (k-s)!$$

where $N_{k,r}$ is the number of permutations with exactly $r$ successful adjacencies.

Hence

$$p_{k,j} = \frac{N_{k,k-j}}{k!}.$$


Step 3: Recurrence for the expectation

If the permutation is already sorted, no shuffle is needed.

That occurs only when $j=1$, with probability $1/k!$.

Otherwise, exactly one shuffle occurs and we continue from the new bundle count.

Thus:

$$E_k = \left(1-\frac1{k!}\right) + \sum_{j=2}^{k} p_{k,j}E_j.$$

Since the term $j=k$ contains $E_k$ itself:

$$E_k(1-p_{k,k}) = \left(1-\frac1{k!}\right) + \sum_{j=2}^{k-1} p_{k,j}E_j.$$

Therefore:

$$E_k = \frac{ 1-\frac1{k!} + \sum_{j=2}^{k-1} p_{k,j}E_j }{ 1-p_{k,k} }.$$

Base case:

$$E_1=0.$$


2. Python implementation

from math import comb, factorial

def expected_shuffles(n):
    # E[k] = expected shuffles for k bundles
    E = [0.0] * (n + 1)

    # Precompute factorials
    fact = [1] * (n + 1)
    for i in range(1, n + 1):
        fact[i] = fact[i - 1] * i

    # Build E[k] iteratively
    for k in range(2, n + 1):

        rhs = 1.0 - 1.0 / fact[k]

        # Contributions from smaller states
        for j in range(2, k):
            r = k - j

            # Inclusion-exclusion count N_{k,r}
            count = 0
            for s in range(r, k):
                count += (
                    (-1) ** (s - r)
                    * comb(s, r)
                    * comb(k - 1, s)
                    * fact[k - s]
                )

            p_kj = count / fact[k]
            rhs += p_kj * E[j]

        # Probability of staying at k bundles
        count_same = 0
        for s in range(k):
            count_same += (
                (-1) ** s
                * comb(k - 1, s)
                * fact[k - s]
            )

        p_kk = count_same / fact[k]

        E[k] = rhs / (1.0 - p_kk)

    return E[n]

# Check given examples
print(expected_shuffles(2))  # 1
print(expected_shuffles(5))  # 4213/871 ≈ 4.836969001148105

# Required answer
ans = expected_shuffles(52)
print(f"{ans:.8f}")

3. Code walkthrough

Factorials

We precompute:

fact[i]

to avoid repeated factorial computation.


Inclusion–exclusion count

For fixed $k$ and $r$:

count += (-1) ** (s - r) * comb(s, r) * comb(k - 1, s) * fact[k - s]

This computes:

$$N_{k,r}$$

the number of permutations with exactly $r$ successful adjacencies.


Dynamic programming

We compute $E_k$ in increasing order:

for k in range(2, n + 1):

because $E_k$ only depends on smaller states.

The self-loop probability $p_{k,k}$ is moved to the denominator:

E[k] = rhs / (1.0 - p_kk)

Verification with examples

The code reproduces:

$$S(2)=1$$

and

$$S(5)=\frac{4213}{871} \approx 4.836969001148105$$

matching the problem statement exactly.


4. Final verification

For $n=52$, the expectation should be much larger than 5 but still modest, because each shuffle tends to create larger bundles and accelerate convergence.

The computed value is:

$$S(52) \approx 54.175293292644746$$

Rounded to 8 decimal places:

Answer: 54.17529329