Problem 253: Tidying Up A

Project Euler Problem 253 Solution

EulerSolve provides an optimized solution for Project Euler Problem 253, Tidying Up A, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary We start with 40 occupied cells in a row. At each step, one occupied cell is chosen uniformly at random and removed. The occupied cells always form a disjoint union of contiguous segments. If $C_t$ denotes the number of occupied segments after the $t$-th removal, the goal is to compute $$\mathbb{E}\!\left[\max_t C_t\right].$$ The key difficulty is that the process has $40!$ possible removal orders, so a direct permutation average is infeasible. The code solves the problem exactly by compressing configurations to segment lengths and then applying a memoized expectation recursion. Mathematical Approach Reverse-Time Viewpoint Every deletion history is just a permutation of the 40 positions. If we reverse that permutation and insert cells instead of deleting them, we obtain the same sequence of occupied sets in reverse order. Therefore the maximum number of segments is identical in the forward deletion process and in the reverse insertion process. This observation is not needed for the main DP, but it explains why the brute-force validator in the C++ file can enumerate insertion permutations for small $n$. State Compression by Segment Lengths A configuration does not need absolute coordinates. Only the lengths of the current occupied segments matter....

Detailed mathematical approach

Problem Summary

A number caterpillar consists of $40$ jigsaw pieces numbered $1$ to $40$. The father picks the pieces up in a random order, all $40!$ orders being equally likely, and puts each piece at its correct place in the line. The pieces already placed form maximal runs of consecutive numbers, the segments; their number starts at $0$ and ends at $1$. Let $M$ be the largest number of segments seen during the tidy-up. We must find the expected value of $M$ for $40$ pieces, rounded to six decimal places. For ten pieces the statement gives the numbers of orders with $M=1,2,3,4,5$, out of $10!=3628800$: $512$, $250912$, $1815264$, $1418112$ and $144000$. So the most likely value is $3$, and the average is $385643/113400\approx3.400732$. For $40$ pieces the most likely value is $11$. ^[1]

The whole solution in one sentence: if we read the order of the pieces backwards, the segments become the runs of pieces that are still missing, and the number of these runs changes in a way that depends only on the list of their lengths, so a recursion over these lists, which remembers the largest number of runs seen so far and stores every value it has computed, gives the exact expected value.

Mathematical Approach

We start with small cases by hand, turn the process around in time, and then set up a recursion over the lengths of the missing runs.

1. What one piece does to the count

Write $n$ for the number of pieces and $s_k$ for the number of segments after $k$ placements, so that $s_0=0$, $s_n=1$ and $M=\max_k s_k$. When piece $x$ is placed, only its neighbours $x-1$ and $x+1$ matter; a neighbour beyond the ends of the line counts as missing. There are three cases:

both neighbours are missing: $x$ starts a new segment, and the count goes up by $1$;
exactly one neighbour is placed: $x$ extends that segment, and the count stays the same;
both neighbours are placed: $x$ joins two segments into one, and the count goes down by $1$.

In the example of the statement, the first pieces $12, 4, 29, 6, 34$ all have both neighbours missing, so the count climbs to $5$; piece $5$ joins $\{4\}$ and $\{6\}$, so it drops to $4$; piece $35$ extends $\{34\}$, so it stays at $4$.

2. Three and four pieces by hand

For $n=3$ the orders $(1,3,2)$ and $(3,1,2)$ give the counts $1,2,1$, because the two end pieces come first and the middle piece joins them; the other four orders give $1,1,1$. So the expected value of $M$ is $(2\cdot2+4\cdot1)/6=4/3$. For $n=4$ the $24$ orders give $M=1$ eight times and $M=2$ sixteen times, so the expected value is $40/24=5/3$.

Two values of $M$ are clear for every $n$. First, $M=1$ happens for exactly $2^{n-1}$ orders: then the placed pieces always form a single segment, so going backwards from the full line, every step takes away its left or its right end, two choices at each of the $n-1$ steps down to a single piece, and different choices give different orders. For $n=10$ this is the $512$ of the statement. Second, $s$ segments need at least $s-1$ missing pieces between them, so $2s-1\le n$ and $M\le\lceil n/2\rceil$, which is $5$ for ten pieces and $20$ for forty; placing the odd pieces first reaches this bound.

3. Why simple counting does not reach forty pieces

There are $40!\approx8.16\cdot10^{47}$ orders, and even the $2^{40}\approx1.1\cdot10^{12}$ sets of placed pieces are too many to store with a running maximum. Averages at one moment are easy, though: after $k$ placements every set of $k$ pieces is equally likely, and a placed piece $x$ starts a segment when $x-1$ is missing or does not exist, which has probability $k/n$ for $x=1$ and $\frac kn\cdot\frac{n-k}{n-1}$ for $x\ge2$. Expected values add up even for dependent events, ^[2] so

$$\begin{aligned}\mathbb E[s_k]&=\frac kn+(n-1)\,\frac{k(n-k)}{n(n-1)}\\&=\frac{k(n-k+1)}{n}.\end{aligned}$$

For $n=10$ this is at most $3$, at $k=5$ and $k=6$, while $M$ averages about $3.4$: the maximum picks the most crowded moment of each tidy-up. Fixed-time averages say nothing about how the counts at different moments go together, so they cannot give the expected maximum.

4. Reading the order backwards

The pieces that are still missing also form maximal runs of consecutive numbers; we call them gaps and write $g_j$ for their number after $j$ placements, so $g_0=1$ and $g_n=0$. Take an order $\sigma=(\sigma_1,\dots,\sigma_n)$ and its reverse $\tau=(\sigma_n,\dots,\sigma_1)$. After $k$ placements of $\tau$ the placed pieces are $\sigma_n,\dots,\sigma_{n-k+1}$, exactly the pieces still missing after $n-k$ placements of $\sigma$. The segments of the one are the gaps of the other:

$$s_k(\tau)=g_{n-k}(\sigma)\qquad(1\le k\le n).$$

Taking the largest value over all moments gives $M(\tau)=G(\sigma)$, where $G=\max_j g_j$ is the largest number of gaps ever seen. Reversing twice gives back the original order, so reversal pairs the $n!$ orders off one to one, and for every $m$ exactly as many orders have $M=m$ as have $G=m$. In particular $\mathbb E[M]=\mathbb E[G]$. For example, $\sigma=(2,1,3)$ leaves the gaps $\{1,2,3\}$, then $\{1\}$ and $\{3\}$, then $\{3\}$, then none, so $g_j=1,2,1,0$ and $G=2$; its reverse $\tau=(3,1,2)$ has the counts $1,2,1$ and $M=2$. From now on we follow the gaps forward in time.

5. Gaps only care about their lengths

Take a gap of length $L$, the missing pieces $u,u+1,\dots,u+L-1$, and place one of them:

if $L=1$, the gap disappears, and the number of gaps goes down by $1$;
if $L\ge2$ and the piece is one of the two end pieces $u$ and $u+L-1$, the gap shrinks to length $L-1$, and the number stays the same;
if the piece is one of the $L-2$ inner pieces, it leaves $a$ missing pieces on its left and $L-1-a$ on its right, with $1\le a\le L-2$: the gap splits into two, and the number goes up by $1$.

What lies beyond the gap, a placed piece or the end of the line, plays no role. For segments it does matter: with pieces $1$ and $2$ missing and piece $3$ placed, placing piece $1$ starts a new segment, although it is an end piece of its gap. A recursion over the segments must keep the two end gaps apart from the others, which for $n=40$ gives about twelve times as many states as there are lists of gap lengths, even with mirror images identified.

The next piece is equally likely to be any missing piece: after $j$ placements each of the $(n-j)!$ ways to finish the order is equally likely, and each missing piece comes next in $(n-j-1)!$ of them. So if the gap lengths form the list $\lambda$, with $r$ missing pieces in total, the next list is $\mu$ with probability $w(\lambda\to\mu)/r$, where the weight $w(\lambda\to\mu)$ counts the missing pieces whose placement turns $\lambda$ into $\mu$. This depends only on $\lambda$, not on the positions of the gaps or on the history: the list of gap lengths is a Markov chain. ^[3] For example, from $(5,2)$ the two pieces of the short gap give $(5,1)$, the two end pieces of the long gap give $(4,2)$, its second and fourth pieces give $(3,2,1)$, and its middle piece gives $(2,2,2)$; the weights $2,2,2,1$ add up to $r=7$.

6. How many lists can occur

A list of gap lengths is a partition of $r$: a way of writing $r$ as a sum of positive integers in which the order of the terms does not matter. ^[4] Not every partition occurs: $g$ gaps need at least $g-1$ placed pieces between them, so $g\le n-r+1$. To count the lists that do occur, imagine an extra piece $n+1$ at the right end that is always placed, and cut the line just after every placed piece. Each part consists of one placed piece and the gap just before it, possibly empty, so the part sizes are the gap lengths plus $1$, together with a $1$ for every placed piece without a gap before it, and they add up to $n+1$. For $n=10$ with the pieces $3$, $4$ and $8$ placed, the gaps $(3,2,2)$ become the parts $4,3,3,1$ of $11$. Conversely, subtracting $1$ from every part larger than $1$ gives the gap lengths back, and every partition of $n+1$ arises: write its parts side by side, each as its missing pieces followed by one placed piece. So the possible lists correspond one to one to the partitions of $n+1$, and there are exactly $p(n+1)$ of them, the partition number: $p(11)=56$ for ten pieces and $p(41)=44583$ for forty, the empty list included. ^[5]

7. A recursion for the expected maximum

The largest number of gaps in a tidy-up is the larger of the largest number $m$ seen so far and the largest number still to come, so its expected value depends only on the current list $\lambda$ and on $m$. Let $E(\lambda,m)$ be this expected value, and write $\ell(\mu)$ for the number of parts of a list $\mu$, its number of gaps. When nothing is missing the maximum is final; otherwise we sum over the possible next lists $\mu$, with $m_\mu=\max(m,\ell(\mu))$:

$$E(\varnothing,m)=m,$$

$$E(\lambda,m)=\sum_{\mu}\frac{w(\lambda\to\mu)}{r}\,E(\mu,m_\mu).$$

The second line is the law of total expectation over the $r$ possible next pieces. ^[6] It holds because of Section 5: once the next list is $\mu$, its future does not depend on how it was reached, and the maximum seen so far becomes $m_\mu$. Every step removes one missing piece, so induction on $r$ shows that the recursion gives the true expected values. At the start the whole line is one gap, so

$$\mathbb E[M]=\mathbb E[G]=E\bigl((n),1\bigr).$$

Many orders lead to the same pair $(\lambda,m)$. Computing each pair only once and storing its value in a table is memoization, the top-down form of dynamic programming. ^[7]

8. A worked example: four pieces

For $n=4$ we start from $\bigl((4),1\bigr)$. The two end pieces of the gap give $(3)$, and the two inner pieces give $(2,1)$, so

$$\begin{aligned}E\bigl((4),1\bigr)&=\tfrac24\,E\bigl((3),1\bigr)\\&\quad+\tfrac24\,E\bigl((2,1),2\bigr).\end{aligned}$$

From $(3)$ the end pieces give $(2)$ and the middle piece gives $(1,1)$. The list $(2)$ can only shrink to $(1)$ and then disappear, so $E((2),1)=1$; from $(1,1)$ the number of gaps can only fall, so $E((1,1),2)=2$. Hence $E((3),1)=\frac23\cdot1+\frac13\cdot2=\frac43$. After $(2,1)$ at most two pieces are missing, so there are never more than two gaps, and $E((2,1),2)=2$. Altogether

$$E\bigl((4),1\bigr)=\tfrac12\cdot\tfrac43+\tfrac12\cdot2=\tfrac53,$$

the value found by hand in Section 2.

9. The checkpoint for ten pieces

The same weights also count orders exactly. Let $C_k(\lambda)$ be the number of ways to place the missing pieces of $\lambda$ one by one so that no later list has more than $k$ gaps. Then $C_k(\varnothing)=1$ and

$$C_k(\lambda)=\sum_{\ell(\mu)\le k}w(\lambda\to\mu)\,C_k(\mu),$$

and for $k\ge1$ the number of orders with $G\le k$ is $C_k((n))$. For $n=10$ this gives $512$, $251424$, $2066688$, $3484800$ and $3628800$ for $k=1,\dots,5$. The differences $512$, $250912$, $1815264$, $1418112$ and $144000$ are exactly the numbers of the statement, which by Section 4 count the orders by $M$ as well as by $G$, and the recursion of Section 7 with exact fractions gives $E((10),1)=385643/113400$.

10. Rounding errors stay tiny

In double precision every arithmetic operation returns its exact result times a factor $1+\delta$ with $|\delta|\le u=2^{-53}\approx1.1\cdot10^{-16}$; here all numbers lie between $1/40$ and $20$, far from underflow and overflow. ^[8] A gap of length $L$ has at most $\lceil L/2\rceil$ different next lists, since the splits for $a$ and $L-1-a$ coincide, so a list with $g$ gaps and $r$ missing pieces has at most $(r+g)/2\le(n+1)/2$ next lists by Section 6, at most $20$ for $n=40$. Each value of the recursion is thus a sum of at most $20$ nonnegative terms, and each term passes through one division, one multiplication and at most $19$ additions: at most $21$ rounding factors per level. Expanding the at most $40$ levels writes the computed value as a sum of positive contributions, each multiplied by at most $840$ such factors, plus one for the final conversion to double in C++. Nothing cancels, so the relative error is at most $(1+u)^{841}-1\lt10^{-13}$, and the absolute error is below $2\cdot10^{-12}$, since the value is at most $20$. This cannot change the sixth decimal unless the exact value lies that close to a rounding boundary, which the checks below exclude.

How the Code Works

The C++ program ^[9] stores a list of gap lengths as a State, a std::vector<int> in increasing order, so that every list has exactly one representation; VecHash hashes it with a multiply-and-xor loop. A Solver holds two hash maps: trans keeps the transitions of every list, shared by all maxima $m$, and memo keeps for every list a row of $41$ values indexed by $m$, with NaN for a value not computed yet.

make_next copies a list, removes one entry of a given length, found with std::lower_bound, and inserts up to two new lengths at their sorted positions. get_transitions builds the weights of Section 5: if the length $L$ occurs cnt times, $L=1$ gives the list without it, with weight cnt, and a longer gap gives the shrunken list, with weight 2 * cnt, and for each $a=1,\dots,L-2$ the split into $a$ and $L-1-a$, with weight cnt. An unordered_map adds up the weights of coinciding splits; for $(40)$ the result has $20$ entries of weight $2$.

expected_max is the recursion of Section 7. It returns the current maximum for the empty list and otherwise the stored value, unless that is NaN, which the test row[current_max] == row[current_max] detects; then it adds up $w/r$ times the values of the next lists, with the maximum raised where needed, and stores the sum. It keeps references into both maps while its recursive calls insert new entries, which is safe because such references survive insertions and rehashing. ^[10]

First, main runs two checks: brute_expectation averages max_segments_of_perm, which follows the segments with the three cases of Section 1, over all $6!$ and all $7!$ orders from std::next_permutation, and two assert statements require the recursion to match these averages, $203/90$ and $229/90$, to within $10^{-12}$. One side counts segments and the other gaps, so the comparisons also test Section 4. With NDEBUG defined, the comparisons are skipped, although both values are still computed. Then expected_max is evaluated for $(40)$ with maximum $1$, and the result is printed with six decimals.

The Python and Java versions implement the same recursion without the checks. Python uses sorted tuples, bisect.insort and one dictionary keyed by the pair of list and maximum; its recursion limit of $2000$ is far more than the $41$ nested calls needed. Java wraps the lengths in a State class with a cached Arrays.hashCode, sorts every new list with Arrays.sort, keeps NaN-marked rows of $41$ doubles and prints with Locale.US. Python and Java use 64-bit doubles; C++ computes in long double, whose precision depends on the compiler and platform (on Apple silicon it equals double), and converts the result to double before printing. Section 10 covers all these cases.

Complexity and Verification

For $n=40$ the recursion visits all $44582$ nonempty lists and evaluates $142008$ pairs $(\lambda,m)$; the loop over next lists runs $1516196$ times in total. The $530383$ stored transitions, each holding its own copy of a list, dominate the memory. On a desktop computer the C++ and Java programs take about $0.2$ seconds and the Python program less than one second. In general there are $p(n+1)$ lists, and by the asymptotic formula of Hardy and Ramanujan ^[5]

$$p(N)\sim\frac{1}{4\sqrt3\,N}\,e^{\pi\sqrt{2N/3}},$$

which grows faster than any power of $N$ but far more slowly than $2^N$. Each list is evaluated for at most $\lceil n/2\rceil$ maxima and has at most $(n+1)/2$ next lists, and each lookup hashes or compares a list of up to $n$ entries. With the usual average-case cost of hash tables the time is $O(n^3p(n+1))$, and the cached transitions, $O(np(n+1))$ lists of length $O(n)$, need $O(n^2p(n+1))$ memory. A simulation would be far weaker: $M$ has a standard deviation of about $1.37$ for $40$ pieces, so a standard error of $5\cdot10^{-7}$ needs about $10^{13}$ random tidy-ups.

Further checks, not part of the published solutions, were run for this article:

For every $n\le9$ all $n!$ orders were enumerated: the identity $s_k(\tau)=g_{n-k}(\sigma)$ holds for every order and every $k$, $M$ and $G$ have the same distribution, and the expected values are $1$, $1$, $4/3$, $5/3$, $59/30$, $203/90$, $229/90$, $3569/1260$ and $14143/4536$.
The formula for $\mathbb E[s_k]$ agrees with a count over all subsets for every $n\le12$, and $p(n+1)$ agrees with a direct count of the partitions with at most $n-r+1$ parts for every $n\le15$ and for $n=40$.
A separately written exact recursion follows the segments forward in time and keeps the two end gaps apart. With whole numbers and no rounding, it gives the same number of orders for every value of $M$ as the gap recursion for every $n\le40$, with the most likely value $11$ for $n=40$.
The C++ program, built with Clang and with GCC, also with NDEBUG, and the Python and Java programs print the same six decimals, which agree with the exact fraction from that distribution. This fraction lies more than $10^{-9}$ from the nearest rounding boundary, and a double-precision run of the recursion differs from it by about $6\cdot10^{-15}$.

Footnotes and References

The references below supply the background facts used above. The reversal argument, the states given by lists of gap lengths, their count and the error estimate are derived explicitly in this article.

Project Euler 253 — Tidying Up A. The official statement: 40 numbered pieces are placed in random order, and M is the largest number of segments seen. For ten pieces it gives the counts of M = 1, …, 5 and their average 385643/113400 ≈ 3.400732; it says that 11 is the most likely value for forty pieces and asks for the average of M for forty pieces, rounded to six decimal places.
Wikipedia — Expected value. Expectation is linear, E[X + Y] = E[X] + E[Y] whenever the right-hand side is defined, also for dependent variables, and the expected value of the indicator of an event is its probability.
Wikipedia — Markov chain. In a Markov chain the probability of moving to the next state depends only on the present state, not on the previous states.
Wikipedia — Integer partition. A partition of n writes n as a sum of positive integers, sums that differ only in the order of their terms being the same partition; p(n) counts them.
OEIS A000041 — The partition numbers. The entry lists p(11) = 56, p(40) = 37338 and p(41) = 44583, and gives the asymptotic formula p(n) ~ e^(π√(2n/3)) / (4n√3) of Hardy and Ramanujan.
Wikipedia — Law of total expectation. For a finite partition of the sample space into events Aᵢ, E[X] = Σ E[X | Aᵢ] P(Aᵢ).
Wikipedia — Dynamic programming. When a naive recursion solves the same subproblems over and over, the top-down approach stores their solutions in a table and checks it before solving a new subproblem.
David Goldberg — What Every Computer Scientist Should Know About Floating-Point Arithmetic. IEEE 754 requires exactly rounded basic operations and p = 53 binary digits in double precision, and rounding to nearest has a relative error of at most ε = (β/2)β^(−p) = 2^(−53).
C++ — Euler253.cpp. The linked immutable revision contains get_transitions, expected_max and the brute-force checks for six and seven pieces; it is the implementation record for this article.
cppreference — std::unordered_map. Insertions invalidate iterators only when they cause a rehash, and references and pointers to stored keys and values are invalidated only by erasing that element.

Problem 253 source code

C++

#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <limits>
#include <numeric>
#include <unordered_map>
#include <vector>

struct VecHash {
    std::size_t operator()(const std::vector<int>& v) const noexcept {
        std::size_t h = 1469598103934665603ULL;
        for (int x : v) {
            h ^= static_cast<std::size_t>(x + 0x9e3779b97f4a7c15ULL);
            h *= 1099511628211ULL;
        }
        return h;
    }
};

using State = std::vector<int>;

struct Solver {
    std::unordered_map<State, std::vector<long double>, VecHash> memo;
    std::unordered_map<State, std::vector<std::pair<State, int>>, VecHash> trans;

    static State make_next(const State& s, int remove_len, int add1, int add2 = 0) {
        State t = s;
        auto it = std::lower_bound(t.begin(), t.end(), remove_len);
        t.erase(it);
        if (add1 > 0) {
            t.insert(std::lower_bound(t.begin(), t.end(), add1), add1);
        }
        if (add2 > 0) {
            t.insert(std::lower_bound(t.begin(), t.end(), add2), add2);
        }
        return t;
    }

    const std::vector<std::pair<State, int>>& get_transitions(const State& s) {
        auto it = trans.find(s);
        if (it != trans.end()) return it->second;

        std::unordered_map<State, int, VecHash> acc;
        const int blocks = static_cast<int>(s.size());
        int i = 0;
        while (i < blocks) {
            int j = i;
            while (j < blocks && s[j] == s[i]) ++j;
            const int len = s[i];
            const int cnt = j - i;

            if (len == 1) {
                State t = make_next(s, 1, 0, 0);
                acc[t] += cnt;
            } else {
                State endpoint = make_next(s, len, len - 1, 0);
                acc[endpoint] += 2 * cnt;
                for (int a = 1; a <= len - 2; ++a) {
                    int b = len - 1 - a;
                    State split = make_next(s, len, a, b);
                    acc[split] += cnt;
                }
            }
            i = j;
        }

        std::vector<std::pair<State, int>> out;
        out.reserve(acc.size());
        for (auto& kv : acc) out.push_back(std::move(kv));
        auto inserted = trans.emplace(s, std::move(out));
        return inserted.first->second;
    }

    long double expected_max(const State& s, int current_max) {
        if (s.empty()) return static_cast<long double>(current_max);
        auto& row = memo[s];
        if (row.empty()) {
            row.assign(41, std::numeric_limits<long double>::quiet_NaN());
        }
        if (row[current_max] == row[current_max]) return row[current_max];

        int remaining = 0;
        for (int x : s) remaining += x;
        const auto& tr = get_transitions(s);

        long double ans = 0;
        for (const auto& [ns, w] : tr) {
            int next_max = std::max(current_max, static_cast<int>(ns.size()));
            ans += (static_cast<long double>(w) / remaining) * expected_max(ns, next_max);
        }
        row[current_max] = ans;
        return ans;
    }
};

static int max_segments_of_perm(const std::vector<int>& p) {
    const int n = static_cast<int>(p.size());
    std::vector<char> on(n + 2, 0);
    int seg = 0;
    int best = 0;
    for (int pos : p) {
        on[pos] = 1;
        const bool l = on[pos - 1];
        const bool r = on[pos + 1];
        if (l && r) --seg;
        else if (!l && !r) ++seg;
        if (seg > best) best = seg;
    }
    return best;
}

static long double brute_expectation(int n) {
    std::vector<int> p(n);
    std::iota(p.begin(), p.end(), 1);
    long double total = 0;
    long double cnt = 0;
    do {
        total += max_segments_of_perm(p);
        cnt += 1;
    } while (std::next_permutation(p.begin(), p.end()));
    return total / cnt;
}

int main() {
    {
        Solver t1;
        long double dp6 = t1.expected_max(State{6}, 1);
        long double br6 = brute_expectation(6);
        assert(std::fabsl(dp6 - br6) < 1e-12L);

        Solver t2;
        long double dp7 = t2.expected_max(State{7}, 1);
        long double br7 = brute_expectation(7);
        assert(std::fabsl(dp7 - br7) < 1e-12L);
    }

    Solver solver;
    long double ans = solver.expected_max(State{40}, 1);
    std::cout << std::fixed << std::setprecision(6) << static_cast<double>(ans) << '\n';
    return 0;
}

Python

import bisect
import sys

# Increase recursion depth just in case, though max depth is 40.
sys.setrecursionlimit(2000)

def make_next(s, remove_len, add1, add2=0):
    t = list(s)
    t.remove(remove_len) # we know it's there
    if add1 > 0:
        bisect.insort(t, add1)
    if add2 > 0:
        bisect.insort(t, add2)
    return tuple(t)

trans_memo = {}

def get_transitions(s):
    if s in trans_memo:
        return trans_memo[s]
        
    acc = {}
    blocks = len(s)
    i = 0
    while i < blocks:
        j = i
        while j < blocks and s[j] == s[i]:
            j += 1
        cnt = j - i
        length = s[i]
        
        if length == 1:
            t = make_next(s, 1, 0, 0)
            acc[t] = acc.get(t, 0) + cnt
        else:
            endpoint = make_next(s, length, length - 1, 0)
            acc[endpoint] = acc.get(endpoint, 0) + 2 * cnt
            for a in range(1, length - 1):
                b = length - 1 - a
                split = make_next(s, length, a, b)
                acc[split] = acc.get(split, 0) + cnt
        i = j
        
    out = list(acc.items())
    trans_memo[s] = out
    return out

expected_memo = {}

def expected_max(s, current_max):
    if not s:
        return float(current_max)
        
    state_key = (s, current_max)
    if state_key in expected_memo:
        return expected_memo[state_key]
        
    remaining = sum(s)
    tr = get_transitions(s)
    
    ans = 0.0
    for ns, w in tr:
        next_max = max(current_max, len(ns))
        ans += (w / float(remaining)) * expected_max(ns, next_max)
        
    expected_memo[state_key] = ans
    return ans

def solve():
    ans = expected_max((40,), 1)
    return "{:.6f}".format(ans)

if __name__ == '__main__':
    print(solve())

Java

import java.util.*;

public class Euler253 {
    static class State implements Comparable<State> {
        int[] pieces;
        int hash;

        State(int[] pieces) {
            this.pieces = pieces;
            this.hash = Arrays.hashCode(pieces);
        }

        @Override
        public int hashCode() {
            return hash;
        }

        @Override
        public boolean equals(Object obj) {
            if (this == obj)
                return true;
            if (obj == null || getClass() != obj.getClass())
                return false;
            State other = (State) obj;
            return Arrays.equals(this.pieces, other.pieces);
        }

        @Override
        public int compareTo(State o) {
            if (this.pieces.length != o.pieces.length)
                return Integer.compare(this.pieces.length, o.pieces.length);
            for (int i = 0; i < this.pieces.length; ++i) {
                if (this.pieces[i] != o.pieces[i])
                    return Integer.compare(this.pieces[i], o.pieces[i]);
            }
            return 0;
        }
    }

    static class Transition {
        State nextState;
        int weight;

        Transition(State nextState, int weight) {
            this.nextState = nextState;
            this.weight = weight;
        }
    }

    static State makeNext(State s, int removeLen, int add1, int add2) {
        int[] old = s.pieces;
        int[] next = new int[old.length + (add1 > 0 ? 1 : 0) + (add2 > 0 ? 1 : 0) - 1];

        boolean removed = false;
        int idx = 0;
        for (int x : old) {
            if (!removed && x == removeLen) {
                removed = true;
            } else {
                next[idx++] = x;
            }
        }
        if (add1 > 0)
            next[idx++] = add1;
        if (add2 > 0)
            next[idx++] = add2;

        Arrays.sort(next);
        return new State(next);
    }

    static Map<State, List<Transition>> transMemo = new HashMap<>();

    static List<Transition> getTransitions(State s) {
        if (transMemo.containsKey(s))
            return transMemo.get(s);

        Map<State, Integer> acc = new HashMap<>();
        int blocks = s.pieces.length;
        int i = 0;

        while (i < blocks) {
            int j = i;
            while (j < blocks && s.pieces[j] == s.pieces[i])
                ++j;
            int cnt = j - i;
            int length = s.pieces[i];

            if (length == 1) {
                State t = makeNext(s, 1, 0, 0);
                acc.put(t, acc.getOrDefault(t, 0) + cnt);
            } else {
                State endpoint = makeNext(s, length, length - 1, 0);
                acc.put(endpoint, acc.getOrDefault(endpoint, 0) + 2 * cnt);
                for (int a = 1; a <= length - 2; ++a) {
                    int b = length - 1 - a;
                    State split = makeNext(s, length, a, b);
                    acc.put(split, acc.getOrDefault(split, 0) + cnt);
                }
            }
            i = j;
        }

        List<Transition> out = new ArrayList<>(acc.size());
        for (Map.Entry<State, Integer> entry : acc.entrySet()) {
            out.add(new Transition(entry.getKey(), entry.getValue()));
        }
        transMemo.put(s, out);
        return out;
    }

    static Map<State, double[]> memo = new HashMap<>();

    static double expectedMax(State s, int currentMax) {
        if (s.pieces.length == 0)
            return (double) currentMax;

        double[] row = memo.get(s);
        if (row == null) {
            row = new double[41];
            Arrays.fill(row, Double.NaN);
            memo.put(s, row);
        }

        if (!Double.isNaN(row[currentMax]))
            return row[currentMax];

        int remaining = 0;
        for (int x : s.pieces)
            remaining += x;

        List<Transition> tr = getTransitions(s);
        double ans = 0.0;

        for (Transition t : tr) {
            int nextMax = Math.max(currentMax, t.nextState.pieces.length);
            ans += ((double) t.weight / remaining) * expectedMax(t.nextState, nextMax);
        }

        row[currentMax] = ans;
        return ans;
    }

    public static String solve() {
        State start = new State(new int[] { 40 });
        double ans = expectedMax(start, 1);
        return String.format(Locale.US, "%.6f", ans);
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}