Detailed mathematical approach

Problem Summary

Start with a uniformly random permutation of \(\{1,2,\dots,n\}\). Whenever adjacent cards form an increasing consecutive run, they are glued into one bundle. The bundles are then shuffled uniformly at random, the gluing rule is applied again, and the process continues until the whole deck becomes one fully sorted bundle.

If \(S(n)\) denotes the expected number of shuffles still needed from that random starting arrangement, the task is to evaluate \(S(52)\) and print it to eight decimal places. The core difficulty is that the deck contains many different visible arrangements, but the implementations show that all of them collapse to a very small state space.

Mathematical Approach

The key observation is that after every gluing phase, each bundle is already internally correct: it contains a consecutive interval of labels in increasing order. That makes it possible to ignore the exact card values inside each bundle and track only how many bundles remain.

Step 1: Reduce the State to the Number of Bundles

Suppose a gluing phase has just finished and there are \(m\) bundles. Read those bundles in increasing label order and rename them \(1,2,\dots,m\). Their lengths may differ, but that no longer matters for the next shuffle.

After a random shuffle, what matters is only the permutation \(\pi\) of these \(m\) renamed bundles. Two neighboring bundles merge exactly when the second one is the immediate successor of the first one in label order. In the relabeled picture, that condition is simply

$\pi(j+1)=\pi(j)+1.$

Therefore every state with \(m\) bundles has the same transition law, independent of the bundle sizes. The process is a Markov chain on bundle counts \(1,2,\dots,n\).

Step 2: Count Arrangements That End with Exactly \(b\) Bundles

Fix a current state with \(m\) bundles and ask for the probability of having exactly \(b\) bundles after the next shuffle and gluing phase.

To finish with \(b\) bundles, the labels \(1,2,\dots,m\) must first be split into \(b\) consecutive intervals. Choosing the \(b-1\) cut positions among the \(m-1\) gaps gives

$\binom{m-1}{b-1}$

possible interval decompositions.

Once these \(b\) intervals are fixed, treat each interval as one macro-block. They may be permuted, but the resulting order must avoid any place where macro-block \(i\) is immediately followed by macro-block \(i+1\), because such a pair would glue again and the decomposition would not be maximal.

Let \(A_b\) be the number of permutations of \(1,2,\dots,b\) with no adjacency of the form \(i\) immediately followed by \(i+1\). Then the number of permutations of the original \(m\) bundles that collapse to exactly \(b\) bundles is

$\binom{m-1}{b-1}A_b.$

Step 3: Evaluate \(A_b\) with Inclusion-Exclusion

For \(1\le i\le b-1\), let \(E_i\) be the event that block \(i\) is immediately followed by block \(i+1\). We want permutations avoiding all these events.

If we force a chosen set of \(r\) such adjacencies, the path \(1,2,\dots,b\) is compressed into \(b-r\) super-blocks. Each super-block has internally fixed order, and the super-blocks themselves can be permuted freely, so the number of permutations satisfying those \(r\) forced adjacencies is

$(b-r)!.$

There are \(\binom{b-1}{r}\) ways to choose which adjacencies are forced. Inclusion-exclusion therefore gives

$A_b=\sum_{r=0}^{b-1}(-1)^r\binom{b-1}{r}(b-r)!.$

Since every permutation of the \(m\) current bundles is equally likely, the transition probability is

$p_{m,b}=\frac{\binom{m-1}{b-1}A_b}{m!},\qquad 1\le b\le m.$

Step 4: Derive the Expected-Value Recurrence

Let \(T_m\) be the expected number of future shuffles needed when a gluing phase has just produced \(m\) bundles. Clearly

$T_1=0,$

because one bundle means the deck is already sorted.

For \(m>1\), we perform one shuffle immediately, then land in some new bundle count \(b\). Hence

$T_m=1+\sum_{b=1}^{m}p_{m,b}T_b.$

The term with \(b=m\) represents the possibility that a shuffle creates no new merge at all. Moving that self-loop term to the left yields the usable recurrence

$T_m=\frac{1+\sum_{b=1}^{m-1}p_{m,b}T_b}{1-p_{m,m}}.$

This allows \(T_2,T_3,\dots,T_n\) to be computed in increasing order.

Step 5: Convert the Random Start into the Same State Model

The problem does not start from a post-shuffle state with a fixed bundle count. It starts from a uniformly random permutation of \(n\) single cards, before any counted shuffle has happened.

If we apply the gluing rule once to that initial random permutation, the resulting distribution of bundle counts is exactly the same as the transition distribution from \(n\) single-card bundles. Therefore

$S(n)=\sum_{b=1}^{n}p_{n,b}T_b.$

This is why the implementations first solve the recurrence for the post-compression expectation \(T_m\), and only then average over the initial distribution to obtain \(S(n)\).

Worked Example: \(m=4\) and \(n=2\)

For \(b=1,2,3,4\), inclusion-exclusion gives

$A_1=1,\qquad A_2=1,\qquad A_3=3,\qquad A_4=11.$

Thus the transition probabilities from \(m=4\) bundles are

$p_{4,1}=\frac{\binom{3}{0}A_1}{4!}=\frac{1}{24},\qquad p_{4,2}=\frac{\binom{3}{1}A_2}{4!}=\frac{3}{24},$

$p_{4,3}=\frac{\binom{3}{2}A_3}{4!}=\frac{9}{24},\qquad p_{4,4}=\frac{\binom{3}{3}A_4}{4!}=\frac{11}{24}.$

These probabilities sum to \(1\), as they must.

For the smallest nontrivial deck, \(n=2\), we have \(p_{2,1}=p_{2,2}=1/2\). Hence

$T_2=\frac{1+\frac{1}{2}T_1}{1-\frac{1}{2}}=2.$

The initial random permutation is already sorted with probability \(1/2\), so

$S(2)=\frac{1}{2}T_1+\frac{1}{2}T_2=1.$

This matches the exact small-case checkpoint used by the method.

How the Code Works

The C++, Python, and Java implementations all follow the same pipeline. First they precompute factorials exactly as integers up to \(52!\), because all transition counts are rational combinations of factorials and binomial coefficients.

Next they evaluate the inclusion-exclusion count \(A_b\) for every \(1\le b\le 52\). With those counts available, they fill a triangular table of transition probabilities \(p_{m,b}\) for all \(1\le b\le m\le 52\).

The counting stage is done exactly, while the divisions are performed in high-precision decimal arithmetic with 100 digits of precision. That is more than enough to stabilize the final rounding to eight decimal places.

After the probability table is built, the implementations compute \(T_2,T_3,\dots,T_{52}\) using the rearranged recurrence. Finally they evaluate

$S(52)=\sum_{b=1}^{52}p_{52,b}T_b$

and print the result with exactly eight digits after the decimal point. The same exact recurrence also reproduces useful checkpoints such as \(S(1)=0\), \(S(2)=1\), and \(S(5)=4213/871\).

Complexity Analysis

Let \(N=52\), or more generally let \(N\) be the maximum deck size to be evaluated. Precomputing factorials costs \(O(N)\) big-integer operations. Computing all values \(A_b\) by inclusion-exclusion costs \(O(N^2)\). Filling the transition table \(p_{m,b}\) also costs \(O(N^2)\), and the dynamic-programming recurrence for \(T_m\) is another \(O(N^2)\).

So, ignoring the internal cost of high-precision arithmetic, the overall algorithm uses \(O(N^2)\) arithmetic steps and \(O(N^2)\) memory because the full triangular probability table is stored. For \(N=52\), this is easily fast enough.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=595
2. Inclusion-exclusion principle: Wikipedia — Inclusion-exclusion principle
3. Composition of an integer: Wikipedia — Composition (combinatorics)
4. Permutation: Wikipedia — Permutation
5. Markov chain: Wikipedia — Markov chain