Problem 298: Selective Amnesia

Project Euler Problem 298 Solution

EulerSolve provides an optimized solution for Project Euler Problem 298, Selective Amnesia, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Larry and Robin observe the same random sequence of calls from an alphabet of size $10$. Both memories have capacity $5$, but they update differently. After $50$ turns we want the exact value of $$\mathbb{E}(|L-R|),$$ where $L$ and $R$ are the total hit counts of Larry and Robin. Mathematical Approach 1. Exact state and scoring rules A state is the ordered pair of memories $$\sigma=(M_L,M_R), \qquad |M_L|,|M_R|\le 5.$$ Both memories keep order from oldest to newest. The update rules are: Larry. On a hit, the called label is removed from its current position and moved to the back. On a miss, the label is appended, and if the memory is full the oldest label is evicted first. Robin. On a hit, the score increases but the memory order does not change. On a miss, Robin also appends and evicts the oldest label if necessary. For one call, the score increment is $$\Delta=(\text{Larry hit})-(\text{Robin hit})\in\{-1,0,+1\}.$$ 2. Why canonical relabeling is valid The actual digit names do not matter. Only equality relations such as "present in Larry only", "present in both", and relative order inside each memory affect future behavior. Therefore two states that differ only by a permutation of labels have exactly the same transition structure....

Detailed mathematical approach

Problem Summary

Larry and Robin play a memory game. In each of $50$ turns a number from $1$ to $10$ is called, the calls being independent and uniform as usual. Each player remembers at most $5$ numbers and starts with an empty memory. A player whose memory contains the called number scores a point; otherwise he adds the number to his memory, first removing another one if the memory is full. Larry removes the number that has not been called for the longest time, Robin the number that has been in his memory for the longest time. With $L$ and $R$ their scores after $50$ turns, we must find the expected value of $|L-R|$, rounded to eight decimal places. The example game of the statement, with the calls $1,2,4,6,1,8,10,2,4,1$, ends with $L=2$ and $R=3$. ^[1] In computer science the two strategies are the cache replacement policies LRU, least recently used, and FIFO, first in first out. ^[2]

The whole solution in one sentence: we write both memories as short ordered lists, rename the numbers so that states that differ only in the names of the numbers become one state, which leaves just $439$ states, and push the probability of every pair of state and score difference through the $50$ turns.

Mathematical Approach

Throughout, $D_t$ denotes $L-R$ after $t$ turns, and $\Delta_t=D_t-D_{t-1}$ is the change caused by the call at turn $t$; it is $+1$, $-1$ or $0$.

1. The memories as ordered lists

We write Larry's memory as a list ordered by the time of the last call, the least recently called number first, and Robin's memory as a list ordered by the time of insertion, the oldest number first. A call of $x$ changes the lists as follows.

Larry. If $x$ is in his list, he scores and $x$ moves to the end, as it is now the most recently called number. Otherwise $x$ is appended, after removing the first entry, the number not called for the longest time, if the list is full.
Robin. If $x$ is in his list, he scores and the list stays as it is, since a call does not change when a number was inserted. Otherwise $x$ is appended, after removing the first entry, the oldest number, if the list is full.

Three facts follow by induction on the turns. Larry's memory always consists of the five most recently called different numbers, or of all numbers called so far while there are fewer: a hit leaves this set unchanged, and a miss adds the new number and drops the least recently called of the old five, which is the number Larry removes. Until the first removal both memories contain exactly the numbers called so far, and from then on both are full, so they always have the same size. And both contain the number called last.

2. The example game as lists

After the calls $1,2,4,6,1,8$ of the example, Larry's list is $(2,4,6,1,8)$, since $1$ was called again at turn 5, while Robin's is $(1,2,4,6,8)$. The new number $10$ at turn 7 removes $2$ from Larry's list but $1$ from Robin's. Then the calls $2$ and $4$ are hits for Robin only, while Larry removes $4$ and $6$, and the call $1$ is a hit for Larry only, while Robin removes $2$. The final lists $(8,10,2,4,1)$ and $(4,6,8,10,1)$ hold the sets of the statement, with $L=2$ and $R=3$.

3. Equal chances: why the absolute value is asked for

Before each call both memories hold the same number $m$ of numbers, and the call is independent of the past, so each player scores with probability $m/10$, whatever happened before. Given the first $t-1$ calls, $\Delta_t$ therefore has conditional mean $0$: the differences $D_t$ form a martingale, ^[5] so $\mathbb E\,D_t=0$ and $\mathbb E\,L=\mathbb E\,R$ after any number of turns: on average the strategies are equally good. If $k$ numbers are in Larry's memory only, then, as the sizes agree, $k$ numbers are in Robin's memory only; a call of one of Larry's $k$ gives $\Delta_t=+1$, a call of one of Robin's $k$ gives $\Delta_t=-1$, and every other call is a hit for both or a miss for both. The distribution of $D_{50}$ is not symmetric, though: Larry is ahead with probability about $0.4331$ and Robin with probability about $0.3889$. Equal expected scores do not determine $\mathbb E|D_{50}|$; Section 10 computes it from the joint distribution of state and difference.

4. The first turn at which the scores can differ

Until the first removal both memories hold the same set, so $\Delta_t=0$. The first removal comes with the sixth different number. If it is called at turn 6, all six calls differ, the least recently called number is also the oldest, both players remove the first number, and the memories stay equal. So $\Delta_t=0$ for $t\le7$, and $\Delta_8\ne0$ needs different memories after turn 7: the first six calls consist of five different numbers $a,b,c,d,e$, in the order of their first calls, and one repetition, and the seventh is a new number $f$. Robin removes $a$, and Larry removes a different number, namely $b$, exactly when the repetition is a call of $a$ at turn 3, 4, 5 or 6. With $4$ turns for the repetition and $10\cdot9\cdot8\cdot7\cdot6\cdot5$ choices of the numbers there are $604800$ such sequences. After each of them the only difference is that Larry holds $a$ and Robin holds $b$, and the eighth call changes the difference when it is $a$ or $b$, so

$$\mathbb E|D_8|=\frac{604800}{10^7}\cdot\frac{2}{10}=\frac{189}{15625},$$

that is, $0.012096$.

5. A Markov chain on pairs of lists

The present pair of lists and the next call, which is independent of the past, determine the next pair and the change, so the pair of lists together with $D_t$ is a Markov chain: the probability of the next state depends only on the present one. ^[3] Its distribution after $50$ turns gives $\mathbb E|D_{50}|$, but $67\,258\,271$ pairs of lists can occur (Section 11), each with up to $101$ values of $D_t$.

6. Renaming the numbers

The rules of Section 1 only ask whether two entries are equal, so the names of the numbers are irrelevant. We rename them $0,1,2,\dots$ in the order in which they first appear when Larry's list and then Robin's list are read from left to right, and call the result the canonical form of the state. In the final state of the example, Larry's list $(8,10,2,4,1)$ becomes $(0,1,2,3,4)$ and Robin's list $(4,6,8,10,1)$ becomes $(3,5,0,1,4)$. Two states have the same canonical form exactly when one is a renaming of the other, a permutation $\pi$ of $\{1,\dots,10\}$ applied to every entry: a renaming does not change where entries first appear, and two states with the same canonical form both arise from it by giving its labels distinct names.

7. Why renamed states are enough

Write $T(s,x)$ for the state after a call of $x$ in the state $s$, and $\Delta(s,x)$ for the change it causes. Since the rules only compare entries, every renaming $\pi$ satisfies

$$\begin{aligned}T(\pi s,\pi x)&=\pi\,T(s,x),\\ \Delta(\pi s,\pi x)&=\Delta(s,x).\end{aligned}$$

For a canonical form $c'$ and a change $\delta$, let $X(s)$ be the set of calls $x$ for which $T(s,x)$ has the canonical form $c'$ and $\Delta(s,x)=\delta$. By the display, $x\mapsto\pi x$ maps $X(s)$ bijectively onto $X(\pi s)$, so the probability $|X(s)|/10$ of this move is the same for all states with the same canonical form. This is the lumpability condition of Kemeny and Snell for the classes of equal canonical form. ^[4] Summing over the states of each class, the probability $P_t(c,d)$ that after $t$ turns the state has canonical form $c$ and $D_t=d$ therefore obeys a recurrence on canonical states only (Section 10).

8. The transitions of a canonical state

Let a canonical state use the $u$ labels $0,\dots,u-1$. For each of them, called with probability $1/10$, we apply the rules, rename the result and record the change: $+1$ for a label in Larry's list only, $-1$ for one in Robin's list only, $0$ for one in both. The other $10-u$ numbers are in neither memory; two of them differ only by a renaming that fixes the state, so by Section 7 they lead to the same canonical form, and both players miss them. One representative, the label $u$, is simulated with weight $10-u$. Calls that lead to the same canonical form with the same change are merged into one transition whose weight counts them; the weights of every state add up to $10$, and a transition of weight $w$ is taken with probability $w/10$.

9. Worked example: three numbers and memories of two

With the numbers $1,2,3$ and memories of two places the renamed chain has only five states, written as (Larry's list, Robin's list):

$A=((),())$, the start: all three calls lead to $B$.
$B=((0),(0))$: the call of $0$ leads back to $B$, the two new numbers to $C$.
$C=((0,1),(0,1))$, both lists in the same order: the call of $1$ and the new number lead back to $C$, and the call of $0$ leads to $D$.
$D=((0,1),(1,0))$, the same numbers in opposite orders: the call of $0$ leads to $C$, the call of $1$ back to $D$, and the new number to $E$, since Larry removes $0$ and Robin removes $1$.
$E=((0,1),(2,1))$, where each player holds one number the other lacks: the call of $0$ leads to $C$ with change $+1$, the call of $2$ to $D$ with change $-1$, the call of $1$ back to $E$.

The difference stays $0$ until the chain reaches $E$. After four turns the chain is in $B$, $C$, $D$ or $E$ with probabilities $\frac1{27}$, $\frac{16}{27}$, $\frac8{27}$ and $\frac2{27}$; the only path to $E$ is $A\to B\to C\to D\to E$, of probability $1\cdot\frac23\cdot\frac13\cdot\frac13$. The fifth call then changes the difference with probability $\frac23$, so

$$\mathbb E|D_5|=\frac{2}{27}\cdot\frac{2}{3}=\frac{4}{81}.$$

Two more turns of the recurrence of Section 10 give $\mathbb E|D_6|=32/243$ and $\mathbb E|D_7|=164/729\approx0.2249657$, and listing all $3^5$, $3^6$ and $3^7$ call sequences gives the same three values.

10. Pushing the distribution through the turns

Let $c_0$ be the empty state. Then $P_0(c_0,0)=1$, and for every canonical state $c'$ and every $d$

$$P_{t+1}(c',d)=\sum\frac{w}{10}\,P_t(c,d-\delta),$$

where the sum runs over the transitions from any state $c$ to $c'$, with change $\delta$ and weight $w$. Since $|\Delta_t|\le1$, after $t$ turns the difference satisfies $-t\le d\le t$, and

$$\mathbb E|D_{50}|=\sum_{c}\sum_{d=-50}^{50}|d|\,P_{50}(c,d).$$

Each of the $10^t$ sequences of $t$ calls has probability $10^{-t}$. With $w$ in place of $w/10$ the recurrence therefore counts call sequences in integers, so $\mathbb E|D_{50}|$ is an integer divided by $10^{50}$, which we use to check the floating-point result.

11. Why there are only 439 states

The renamed chain is small because Larry and Robin never disagree about more than two numbers each. Suppose both memories are full; before that they are equal. Let $x$ be the number in Larry's memory whose last call, at turn $\tau$, is the oldest. Every number called after turn $\tau$ was called more recently than $x$, so by Section 1 it is one of the other four numbers in Larry's memory.

Any five consecutive misses of Robin are calls of five different numbers: a number he inserts is still there at his next four misses, because each miss removes only the oldest entry. Hence Robin misses at most $j\le4$ times after turn $\tau$, and the numbers inserted at these misses are still in his memory and belong to Larry's.
Each of the other $5-j$ numbers in Larry's memory was in Robin's memory right after turn $\tau$: $x$ was called then, and any other one was called later, was in Robin's memory after that call and was not inserted in between, so it had been there all along.
Since turn $\tau$ at most $j$ numbers have left Robin's memory, at most one per miss.

So at most $\min(j,5-j)\le2$ numbers in Larry's memory are missing from Robin's, and a canonical state uses at most $7$ labels. Listing all canonical states reachable from the empty state, as the program does, gives $154$, $240$ and $45$ states with $k=0$, $1$ and $2$. For $k=0$ both lists hold the same $m\le5$ numbers, and each of the $m!$ relative orders occurs: call the numbers in Robin's order and then once more in Larry's order, hits that rearrange only Larry's list; so $0!+1!+\dots+5!=154$. A canonical state with $u$ labels stands for the $10!/(10-u)!$ ways of giving its labels distinct names, all of which occur, and over the $439$ states these numbers add up to the $67\,258\,271$ pairs of lists of Section 5. For memories of $M$ places the argument gives $k\le\lfloor M/2\rfloor$, attained for every $M\le7$.

How the Code Works

The C++ program ^[7] reads the optional parameters --alphabet=, --capacity= and --turns=, by default $10$, $5$ and $50$, with at most five places and $15$ numbers, builds the renamed chain with build_graph, tests it with run_checkpoints (next section) and computes the expectation with solve_exact.

A State holds the lengths l_size, r_size and the arrays l, r of the two lists, first entry first. apply_call_in_place applies Section 1 to one call by shifting array entries and returns the change in a CallResult, and canonicalize_state renames the labels as in Section 6. build_graph explores the chain from the empty state by breadth-first search with a queue. ^[6] intern_state packs a canonical state into a 64-bit key with encode_state, four bits per entry, looks it up in an unordered_map and, if the key is present, confirms with same_state that it encodes the same lists; a new state receives the next number. expand_state queues new states, simulates every present label with weight $1$ and the smallest absent label with weight outside_count, which is $10-u$, and push_or_merge merges calls with the same target and change into one Transition (Section 8). This gives $439$ states and $2886$ transitions, merged from $2913$ simulated calls.

solve_exact is the recurrence of Section 10. It keeps two arrays current and next of long double, one row of $2T+1$ entries per state, with the entry for $d$ at position $d+T$. At turn $t$, for every state, every transition and every $d$ with $|d|\le t-1$, it adds the entry times weight divided by alphabet to the entry of the target state at $d+\delta$, skipping zero entries, and then swaps the arrays. It returns $\sum|d|\,P_T(c,d)$ and the total probability. The program prints the expectation with $15$ decimals and rounded to eight, the total probability, the number of states and the running time.

The Python and Java versions run the same construction and recurrence for the default parameters, without options and checks. Python's solve calls build_graph and solve_exact; there apply_call_in_place returns a new tuple, and the canonical tuples themselves are the dictionary keys. Java's main calls buildGraph, which copies a State before applyCall and canonicalizeState and uses State and Transition, with their own equals and hashCode, as HashMap keys, and then solveExact. Both use 64-bit floating-point numbers and print only the expectation with $15$ decimals.

Complexity and Verification

At turn $t$ the loops pair each of the $2886$ transitions with $2t-1$ values of $d$; as $\sum_{t=1}^{50}(2t-1)=50^2$, that makes $2886\cdot50^2=7\,215\,000$ candidate updates, of which the $2\,779\,770$ with a nonzero source entry cost a multiply-add. All three programs need well under a second, and each array holds $439\cdot101=44\,339$ numbers. In general the method takes $O(E\,T^2)$ time and $O(S\,T+E)$ memory for $S$ states, $E$ transitions and $T$ turns; given enough numbers, memories of $1$ to $7$ places give $2$, $5$, $16$, $73$, $439$, $3229$ and $27754$ states.

For the default parameters run_checkpoints checks the labels and weights of every state, replays the example game turn by turn, compares solve_exact with brute_force_expectation, an enumeration of all call sequences, for the case of Section 9 with seven turns and for four numbers, memories of two and eight turns, and compares the main result with a stored reference value to within $10^{-15}$; --skip-checkpoints turns the checks off. On many systems long double is wider than double; ^[8] on our arm64 macOS test machine the two agree, and all checks pass. The stored value matches such a double-precision run but lies $1.7\cdot10^{-14}$ from the exact value, so with a wider long double the more accurate result misses it by more than the tolerance, as an emulation with 64-bit and 113-bit significands confirms, and the program stops without a result.

For this article we also verified:

Clang 22 and GCC 16 builds pass the checks, with a total probability about $5\cdot10^{-15}$ below $1$, and Python and Java print the same value to all $15$ decimals.
A Python program with call and insertion times instead of lists agrees with the list rules in $20\,000$ random games of various sizes and confirms the facts of Section 1.
An exact integer recurrence with its own construction of the chain finds the same chain, $D_t=0$ for $t\le7$, $\mathbb E|D_8|=189/15625$ and $\mathbb E\,D_t=0$ for all $t\le50$; its $\mathbb E|D_{50}|$ differs from the values printed by all four builds by about $1.7\cdot10^{-14}$ and rounds to the same eight decimals.
Enumerating all call sequences reproduces the recurrence in nine small cases, such as $531/2048$ for four numbers, memories of two and eight turns. A recurrence on pairs of lists without renaming agrees for six numbers with memories of three and for seven and eight numbers with memories of four, whose $2947$, $127\,352$ and $344\,857$ reachable pairs match the counting rule of Section 11.
Four runs of $10^8$ games each, simulated with the literal rules, land between $0.3$ and $1.4$ standard errors from the exact value. Since $|D_{50}|$ has standard deviation $1.4207$, that standard error is $1.4\cdot10^{-4}$, far too coarse for eight decimals.

Footnotes and References

The references below supply background facts; the list rules, the renaming argument, the bound of Section 11 and all numbers are worked out in this article.

Project Euler 298 — Selective Amnesia. The official statement defines the game, its two removal strategies and a ten-turn example, and asks for the expected value of |L − R| after 50 turns to eight decimal places.
Wikipedia — Cache replacement policies. FIFO evicts blocks in the order in which they were added, regardless of how often they were accessed, and LRU discards the least recently used items first.
Wikipedia — Markov chain. In a discrete-time Markov chain the probability of moving to the next state depends only on the present state and not on the previous states.
Wikipedia — Lumpability. A method of Kemeny and Snell for reducing the state space: a transition matrix is lumpable on a partition of the states if all states of one block have the same total probability of moving into each block.
Wikipedia — Martingale (probability theory). Y₁, Y₂, … is a martingale with respect to X₁, X₂, … if E|Yₙ| is finite and the expected value of Yₙ₊₁ given X₁, …, Xₙ equals Yₙ; here Yₙ is L − R after n turns and Xₙ the n-th call.
Wikipedia — Breadth-first search. Breadth-first search explores all nodes at the present depth before moving on to the next level, keeps its frontier in a queue and takes O(|V| + |E|) time.
C++ — Euler298.cpp. The linked immutable revision contains build_graph, solve_exact and run_checkpoints. It is the implementation record for this article.
cppreference — Fundamental types. The extended precision type long double need not follow IEEE-754: many x86 and x86-64 implementations use the x87 80-bit format, some ARM64 implementations binary128, and MSVC the format of double.

Problem 298 source code

C++

#include <algorithm>
#include <array>
#include <chrono>
#include <cmath>
#include <cstdint>
#include <cstdlib>
#include <iomanip>
#include <iostream>
#include <limits>
#include <functional>
#include <queue>
#include <string>
#include <thread>
#include <unordered_map>
#include <vector>

namespace {

using u32 = std::uint32_t;
using u64 = std::uint64_t;

constexpr int kDefaultAlphabet = 10;
constexpr int kDefaultCapacity = 5;
constexpr int kDefaultTurns = 50;
constexpr int kMaxCapacity = 5;
constexpr int kMaxLabelsEncodable = 15;

constexpr long double kMainExpected = 1.7688229423800932L;
constexpr long double kToleranceTight = 1e-15L;
constexpr long double kToleranceLoose = 1e-12L;
constexpr long double kMassTolerance = 5e-14L;

struct Options {
    int alphabet = kDefaultAlphabet;
    int capacity = kDefaultCapacity;
    int turns = kDefaultTurns;
    bool run_checkpoints = true;
    bool allow_multithreading = true;
    unsigned requested_threads = 0U;
};

struct State {
    int l_size = 0;
    int r_size = 0;
    std::array<int, kMaxCapacity> l{};
    std::array<int, kMaxCapacity> r{};
};

struct CallResult {
    bool larry_hit = false;
    bool robin_hit = false;
    int delta = 0;  // Larry hit minus Robin hit.
};

struct Transition {
    int next_state = -1;
    int delta = 0;
    int weight = 0;
};

struct Graph {
    int alphabet = 0;
    int capacity = 0;
    std::vector<State> states;
    std::vector<std::vector<Transition>> transitions;
};

struct ExactSolveResult {
    long double expectation = 0.0L;
    long double probability_mass = 0.0L;
};

bool parse_u32_after_prefix(const std::string& arg, const char* prefix, u32& value) {
    const std::string p(prefix);
    if (arg.rfind(p, 0) != 0U) {
        return false;
    }

    const std::string tail = arg.substr(p.size());
    if (tail.empty()) {
        return false;
    }

    u64 parsed = 0ULL;
    for (const char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10ULL + static_cast<u64>(c - '0');
        if (parsed > static_cast<u64>(std::numeric_limits<u32>::max())) {
            return false;
        }
    }

    value = static_cast<u32>(parsed);
    return true;
}

bool parse_unsigned_after_prefix(const std::string& arg,
                                 const char* prefix,
                                 unsigned& value) {
    u32 parsed = 0U;
    if (!parse_u32_after_prefix(arg, prefix, parsed)) {
        return false;
    }
    value = static_cast<unsigned>(parsed);
    return true;
}

bool parse_int_after_prefix(const std::string& arg, const char* prefix, int& value) {
    u32 parsed = 0U;
    if (!parse_u32_after_prefix(arg, prefix, parsed)) {
        return false;
    }
    if (parsed > static_cast<u32>(std::numeric_limits<int>::max())) {
        return false;
    }
    value = static_cast<int>(parsed);
    return true;
}

bool parse_arguments(const int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);

        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (arg == "--single-thread") {
            options.allow_multithreading = false;
            continue;
        }
        int parsed_int = 0;
        if (parse_int_after_prefix(arg, "--alphabet=", parsed_int)) {
            options.alphabet = parsed_int;
            continue;
        }
        if (parse_int_after_prefix(arg, "--capacity=", parsed_int)) {
            options.capacity = parsed_int;
            continue;
        }
        if (parse_int_after_prefix(arg, "--turns=", parsed_int)) {
            options.turns = parsed_int;
            continue;
        }

        unsigned parsed_unsigned = 0U;
        if (parse_unsigned_after_prefix(arg, "--threads=", parsed_unsigned)) {
            options.requested_threads = parsed_unsigned;
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }

    if (options.alphabet < 1) {
        std::cerr << "--alphabet must be at least 1.\n";
        return false;
    }
    if (options.capacity < 1 || options.capacity > kMaxCapacity) {
        std::cerr << "--capacity must be in [1, " << kMaxCapacity << "].\n";
        return false;
    }
    if (options.capacity > options.alphabet) {
        std::cerr << "--capacity cannot exceed --alphabet.\n";
        return false;
    }
    if (options.turns < 0) {
        std::cerr << "--turns must be non-negative.\n";
        return false;
    }
    if (options.alphabet > kMaxLabelsEncodable) {
        std::cerr << "--alphabet must be <= " << kMaxLabelsEncodable
                  << " (encoding limit).\n";
        return false;
    }

    return true;
}

unsigned choose_thread_count(const bool allow_multithreading,
                             const unsigned requested_threads,
                             const std::size_t workload) {
    if (!allow_multithreading || workload < 50'000ULL) {
        return 1U;
    }

    unsigned threads = requested_threads;
    if (threads == 0U) {
        threads = std::thread::hardware_concurrency();
        if (threads == 0U) {
            threads = 1U;
        }
    }

    return std::max(1U, std::min<unsigned>(threads, 32U));
}

int find_index(const int* arr, const int size, const int value) {
    for (int i = 0; i < size; ++i) {
        if (arr[i] == value) {
            return i;
        }
    }
    return -1;
}

void canonicalize_state(State& state) {
    std::array<int, kMaxLabelsEncodable> remap{};
    remap.fill(-1);

    int next_id = 0;

    for (int i = 0; i < state.l_size; ++i) {
        const int x = state.l[static_cast<std::size_t>(i)];
        if (remap[static_cast<std::size_t>(x)] < 0) {
            remap[static_cast<std::size_t>(x)] = next_id++;
        }
    }
    for (int i = 0; i < state.r_size; ++i) {
        const int x = state.r[static_cast<std::size_t>(i)];
        if (remap[static_cast<std::size_t>(x)] < 0) {
            remap[static_cast<std::size_t>(x)] = next_id++;
        }
    }

    for (int i = 0; i < state.l_size; ++i) {
        const int x = state.l[static_cast<std::size_t>(i)];
        state.l[static_cast<std::size_t>(i)] = remap[static_cast<std::size_t>(x)];
    }
    for (int i = 0; i < state.r_size; ++i) {
        const int x = state.r[static_cast<std::size_t>(i)];
        state.r[static_cast<std::size_t>(i)] = remap[static_cast<std::size_t>(x)];
    }
}

u64 encode_state(const State& state, const int capacity) {
    u64 key = 0ULL;
    key |= static_cast<u64>(state.l_size);
    key |= (static_cast<u64>(state.r_size) << 3U);

    unsigned shift = 6U;
    for (int i = 0; i < capacity; ++i) {
        const u64 nibble = (i < state.l_size)
                               ? static_cast<u64>(state.l[static_cast<std::size_t>(i)] + 1)
                               : 0ULL;
        key |= (nibble << shift);
        shift += 4U;
    }

    for (int i = 0; i < capacity; ++i) {
        const u64 nibble = (i < state.r_size)
                               ? static_cast<u64>(state.r[static_cast<std::size_t>(i)] + 1)
                               : 0ULL;
        key |= (nibble << shift);
        shift += 4U;
    }

    return key;
}

bool same_state(const State& a, const State& b) {
    if (a.l_size != b.l_size || a.r_size != b.r_size) {
        return false;
    }

    for (int i = 0; i < a.l_size; ++i) {
        if (a.l[static_cast<std::size_t>(i)] != b.l[static_cast<std::size_t>(i)]) {
            return false;
        }
    }
    for (int i = 0; i < a.r_size; ++i) {
        if (a.r[static_cast<std::size_t>(i)] != b.r[static_cast<std::size_t>(i)]) {
            return false;
        }
    }
    return true;
}

bool has_contiguous_labels(const State& state) {
    std::array<unsigned char, kMaxLabelsEncodable> seen{};
    seen.fill(0U);

    int count = 0;
    int mx = -1;
    for (int i = 0; i < state.l_size; ++i) {
        const int x = state.l[static_cast<std::size_t>(i)];
        if (seen[static_cast<std::size_t>(x)] == 0U) {
            seen[static_cast<std::size_t>(x)] = 1U;
            ++count;
            mx = std::max(mx, x);
        }
    }
    for (int i = 0; i < state.r_size; ++i) {
        const int x = state.r[static_cast<std::size_t>(i)];
        if (seen[static_cast<std::size_t>(x)] == 0U) {
            seen[static_cast<std::size_t>(x)] = 1U;
            ++count;
            mx = std::max(mx, x);
        }
    }

    return (mx + 1) == count;
}

CallResult apply_call_in_place(State& state, const int called_label, const int capacity) {
    CallResult result;

    {
        const int idx = find_index(state.l.data(), state.l_size, called_label);
        if (idx >= 0) {
            result.larry_hit = true;
            ++result.delta;

            for (int i = idx; i + 1 < state.l_size; ++i) {
                state.l[static_cast<std::size_t>(i)] = state.l[static_cast<std::size_t>(i + 1)];
            }
            state.l[static_cast<std::size_t>(state.l_size - 1)] = called_label;
        } else {
            if (state.l_size < capacity) {
                state.l[static_cast<std::size_t>(state.l_size)] = called_label;
                ++state.l_size;
            } else {
                for (int i = 0; i + 1 < capacity; ++i) {
                    state.l[static_cast<std::size_t>(i)] = state.l[static_cast<std::size_t>(i + 1)];
                }
                state.l[static_cast<std::size_t>(capacity - 1)] = called_label;
            }
        }
    }

    {
        const int idx = find_index(state.r.data(), state.r_size, called_label);
        if (idx >= 0) {
            result.robin_hit = true;
            --result.delta;
        } else {
            if (state.r_size < capacity) {
                state.r[static_cast<std::size_t>(state.r_size)] = called_label;
                ++state.r_size;
            } else {
                for (int i = 0; i + 1 < capacity; ++i) {
                    state.r[static_cast<std::size_t>(i)] = state.r[static_cast<std::size_t>(i + 1)];
                }
                state.r[static_cast<std::size_t>(capacity - 1)] = called_label;
            }
        }
    }

    return result;
}

int intern_state(const State& state,
                 const int capacity,
                 std::unordered_map<u64, int>& id_by_key,
                 std::vector<State>& states,
                 std::vector<std::vector<Transition>>& transitions,
                 std::vector<unsigned char>& expanded,
                 bool* is_new = nullptr) {
    State canonical = state;
    canonicalize_state(canonical);
    const u64 key = encode_state(canonical, capacity);

    const auto it = id_by_key.find(key);
    if (it != id_by_key.end()) {
        const State& existing = states[static_cast<std::size_t>(it->second)];
        if (!same_state(existing, canonical)) {
            std::cerr << "Internal encoding collision detected.\n";
            std::abort();
        }
        if (is_new != nullptr) {
            *is_new = false;
        }
        return it->second;
    }

    const int id = static_cast<int>(states.size());
    id_by_key.emplace(key, id);
    states.push_back(canonical);
    transitions.emplace_back();
    expanded.push_back(0U);

    if (is_new != nullptr) {
        *is_new = true;
    }
    return id;
}

void expand_state(const int state_id,
                  const int alphabet,
                  const int capacity,
                  std::unordered_map<u64, int>& id_by_key,
                  std::vector<State>& states,
                  std::vector<std::vector<Transition>>& transitions,
                  std::vector<unsigned char>& expanded,
                  std::queue<int>& bfs_queue) {
    if (expanded[static_cast<std::size_t>(state_id)] != 0U) {
        return;
    }
    expanded[static_cast<std::size_t>(state_id)] = 1U;

    const State state = states[static_cast<std::size_t>(state_id)];

    std::array<unsigned char, kMaxLabelsEncodable> present{};
    present.fill(0U);

    for (int i = 0; i < state.l_size; ++i) {
        const int x = state.l[static_cast<std::size_t>(i)];
        present[static_cast<std::size_t>(x)] = 1U;
    }
    for (int i = 0; i < state.r_size; ++i) {
        const int x = state.r[static_cast<std::size_t>(i)];
        present[static_cast<std::size_t>(x)] = 1U;
    }

    std::vector<int> known_labels;
    known_labels.reserve(static_cast<std::size_t>(state.l_size + state.r_size));
    for (int x = 0; x < alphabet; ++x) {
        if (present[static_cast<std::size_t>(x)] != 0U) {
            known_labels.push_back(x);
        }
    }

    const int outside_count = alphabet - static_cast<int>(known_labels.size());

    std::vector<Transition> out;
    out.reserve(known_labels.size() + 1ULL);

    auto push_or_merge = [&out](const int next_id, const int delta, const int weight) {
        for (Transition& tr : out) {
            if (tr.next_state == next_id && tr.delta == delta) {
                tr.weight += weight;
                return;
            }
        }
        out.push_back(Transition{next_id, delta, weight});
    };

    for (const int called : known_labels) {
        State next = state;
        const CallResult result = apply_call_in_place(next, called, capacity);
        canonicalize_state(next);

        bool is_new = false;
        const int next_id =
            intern_state(next, capacity, id_by_key, states, transitions, expanded, &is_new);
        if (is_new) {
            bfs_queue.push(next_id);
        }

        push_or_merge(next_id, result.delta, 1);
    }

    if (outside_count > 0) {
        int representative = -1;
        for (int x = 0; x < alphabet; ++x) {
            if (present[static_cast<std::size_t>(x)] == 0U) {
                representative = x;
                break;
            }
        }

        State next = state;
        const CallResult result = apply_call_in_place(next, representative, capacity);
        canonicalize_state(next);

        bool is_new = false;
        const int next_id =
            intern_state(next, capacity, id_by_key, states, transitions, expanded, &is_new);
        if (is_new) {
            bfs_queue.push(next_id);
        }

        push_or_merge(next_id, result.delta, outside_count);
    }

    transitions[static_cast<std::size_t>(state_id)] = std::move(out);
}

Graph build_graph(const int alphabet, const int capacity) {
    Graph graph;
    graph.alphabet = alphabet;
    graph.capacity = capacity;

    std::unordered_map<u64, int> id_by_key;
    id_by_key.reserve(1024U);

    std::vector<unsigned char> expanded;
    expanded.reserve(1024U);

    State initial;
    (void)intern_state(initial,
                       capacity,
                       id_by_key,
                       graph.states,
                       graph.transitions,
                       expanded,
                       nullptr);

    std::queue<int> bfs_queue;
    bfs_queue.push(0);

    while (!bfs_queue.empty()) {
        const int id = bfs_queue.front();
        bfs_queue.pop();

        expand_state(id,
                     alphabet,
                     capacity,
                     id_by_key,
                     graph.states,
                     graph.transitions,
                     expanded,
                     bfs_queue);
    }

    return graph;
}

ExactSolveResult solve_exact(const Graph& graph,
                             const int turns,
                             const bool allow_multithreading,
                             const unsigned requested_threads) {
    const int state_count = static_cast<int>(graph.states.size());
    const int width = 2 * turns + 1;
    const int offset = turns;

    std::vector<long double> current(static_cast<std::size_t>(state_count) *
                                         static_cast<std::size_t>(width),
                                     0.0L);
    std::vector<long double> next(static_cast<std::size_t>(state_count) *
                                      static_cast<std::size_t>(width),
                                  0.0L);

    current[static_cast<std::size_t>(offset)] = 1.0L;

    for (int turn = 1; turn <= turns; ++turn) {
        std::fill(next.begin(), next.end(), 0.0L);

        const int d_prev_min = offset - (turn - 1);
        const int d_prev_max = offset + (turn - 1);

        const std::size_t workload = static_cast<std::size_t>(state_count) *
                                     static_cast<std::size_t>(2 * turn - 1) * 11ULL;
        const unsigned threads = choose_thread_count(allow_multithreading,
                                                     requested_threads,
                                                     workload);

        if (threads == 1U) {
            for (int sid = 0; sid < state_count; ++sid) {
                const std::size_t src_base =
                    static_cast<std::size_t>(sid) * static_cast<std::size_t>(width);

                for (const Transition& tr : graph.transitions[static_cast<std::size_t>(sid)]) {
                    const long double p = static_cast<long double>(tr.weight) /
                                          static_cast<long double>(graph.alphabet);
                    const std::size_t dst_base = static_cast<std::size_t>(tr.next_state) *
                                                 static_cast<std::size_t>(width);

                    for (int d = d_prev_min; d <= d_prev_max; ++d) {
                        const long double prob = current[src_base + static_cast<std::size_t>(d)];
                        if (prob == 0.0L) {
                            continue;
                        }

                        const int nd = d + tr.delta;
                        next[dst_base + static_cast<std::size_t>(nd)] += prob * p;
                    }
                }
            }
        } else {
            std::vector<std::vector<long double>> local_buffers(
                threads,
                std::vector<long double>(static_cast<std::size_t>(state_count) *
                                             static_cast<std::size_t>(width),
                                         0.0L));

            std::vector<std::thread> workers;
            workers.reserve(threads);

            for (unsigned t = 0U; t < threads; ++t) {
                const int begin = (state_count * static_cast<int>(t)) /
                                  static_cast<int>(threads);
                const int end = (state_count * static_cast<int>(t + 1U)) /
                                static_cast<int>(threads);

                workers.emplace_back([&, begin, end, t]() {
                    std::vector<long double>& local = local_buffers[t];

                    for (int sid = begin; sid < end; ++sid) {
                        const std::size_t src_base =
                            static_cast<std::size_t>(sid) * static_cast<std::size_t>(width);

                        for (const Transition& tr :
                             graph.transitions[static_cast<std::size_t>(sid)]) {
                            const long double p = static_cast<long double>(tr.weight) /
                                                  static_cast<long double>(graph.alphabet);
                            const std::size_t dst_base =
                                static_cast<std::size_t>(tr.next_state) *
                                static_cast<std::size_t>(width);

                            for (int d = d_prev_min; d <= d_prev_max; ++d) {
                                const long double prob =
                                    current[src_base + static_cast<std::size_t>(d)];
                                if (prob == 0.0L) {
                                    continue;
                                }

                                const int nd = d + tr.delta;
                                local[dst_base + static_cast<std::size_t>(nd)] += prob * p;
                            }
                        }
                    }
                });
            }

            for (std::thread& worker : workers) {
                worker.join();
            }

            for (std::size_t idx = 0ULL; idx < next.size(); ++idx) {
                long double acc = 0.0L;
                for (unsigned t = 0U; t < threads; ++t) {
                    acc += local_buffers[t][idx];
                }
                next[idx] = acc;
            }
        }

        current.swap(next);
    }

    long double expectation = 0.0L;
    long double mass = 0.0L;

    for (int sid = 0; sid < state_count; ++sid) {
        const std::size_t base = static_cast<std::size_t>(sid) * static_cast<std::size_t>(width);
        for (int i = 0; i < width; ++i) {
            const long double p = current[base + static_cast<std::size_t>(i)];
            if (p == 0.0L) {
                continue;
            }

            const int diff = i - offset;
            expectation += p * static_cast<long double>(std::abs(diff));
            mass += p;
        }
    }

    return ExactSolveResult{expectation, mass};
}

long double brute_force_expectation(const int alphabet, const int capacity, const int turns) {
    long double total_abs_diff = 0.0L;

    std::function<void(int, State&, int)> dfs =
        [&](const int step, State& state, const int diff) {
            if (step == turns) {
                total_abs_diff += static_cast<long double>(std::abs(diff));
                return;
            }

            for (int x = 0; x < alphabet; ++x) {
                State next = state;
                const CallResult result = apply_call_in_place(next, x, capacity);
                dfs(step + 1, next, diff + result.delta);
            }
        };

    State initial;
    dfs(0, initial, 0);

    long double total_sequences = 1.0L;
    for (int i = 0; i < turns; ++i) {
        total_sequences *= static_cast<long double>(alphabet);
    }

    return total_abs_diff / total_sequences;
}

bool compare_sets(const std::array<int, kMaxCapacity>& memory,
                  const int size,
                  std::initializer_list<int> expected_values) {
    std::vector<int> got;
    got.reserve(static_cast<std::size_t>(size));
    for (int i = 0; i < size; ++i) {
        got.push_back(memory[static_cast<std::size_t>(i)]);
    }
    std::sort(got.begin(), got.end());

    std::vector<int> expected(expected_values.begin(), expected_values.end());
    std::sort(expected.begin(), expected.end());

    return got == expected;
}

bool run_example_checkpoint() {
    State game;
    int cumulative_diff = 0;

    const std::vector<int> sequence = {1, 2, 4, 6, 1, 8, 10, 2, 4, 1};
    const std::vector<int> expected_larry = {0, 0, 0, 0, 1, 1, 1, 1, 1, 2};
    const std::vector<int> expected_robin = {0, 0, 0, 0, 1, 1, 1, 2, 3, 3};

    int larry_score = 0;
    int robin_score = 0;

    for (std::size_t i = 0ULL; i < sequence.size(); ++i) {
        const CallResult result = apply_call_in_place(game, sequence[i], kDefaultCapacity);
        cumulative_diff += result.delta;

        if (result.larry_hit) {
            ++larry_score;
        }
        if (result.robin_hit) {
            ++robin_score;
        }

        if (larry_score != expected_larry[i] || robin_score != expected_robin[i]) {
            std::cerr << "Example checkpoint mismatch on turn " << (i + 1ULL)
                      << ": got (L=" << larry_score << ", R=" << robin_score
                      << "), expected (L=" << expected_larry[i]
                      << ", R=" << expected_robin[i] << ").\n";
            return false;
        }
    }

    if (cumulative_diff != -1) {
        std::cerr << "Example checkpoint mismatch in total difference: got "
                  << cumulative_diff << ", expected -1.\n";
        return false;
    }

    if (!compare_sets(game.l, game.l_size, {1, 2, 4, 8, 10})) {
        std::cerr << "Example checkpoint mismatch in Larry's final memory set.\n";
        return false;
    }

    if (!compare_sets(game.r, game.r_size, {1, 4, 6, 8, 10})) {
        std::cerr << "Example checkpoint mismatch in Robin's final memory set.\n";
        return false;
    }

    return true;
}

bool run_bruteforce_cross_checkpoint() {
    const struct {
        int alphabet;
        int capacity;
        int turns;
    } cases[] = {
        {3, 2, 7},
        {4, 2, 8},
    };

    for (const auto& test_case : cases) {
        const Graph graph = build_graph(test_case.alphabet, test_case.capacity);
        const ExactSolveResult exact = solve_exact(graph, test_case.turns, false, 1U);
        const long double brute =
            brute_force_expectation(test_case.alphabet, test_case.capacity, test_case.turns);

        if (std::fabs(exact.expectation - brute) > 1e-15L) {
            std::cerr << "Brute-force cross-check failed for (alphabet="
                      << test_case.alphabet << ", capacity=" << test_case.capacity
                      << ", turns=" << test_case.turns << "): exact="
                      << std::setprecision(18) << exact.expectation
                      << ", brute=" << brute << '\n';
            return false;
        }

        if (std::fabs(exact.probability_mass - 1.0L) > kMassTolerance) {
            std::cerr << "Probability mass drift in brute-force cross-check case.\n";
            return false;
        }
    }

    return true;
}

bool run_thread_consistency_checkpoint(const Graph& graph) {
    unsigned hw = std::thread::hardware_concurrency();
    if (hw < 2U) {
        return true;
    }

    const ExactSolveResult single = solve_exact(graph, kDefaultTurns, false, 1U);
    const ExactSolveResult threaded = solve_exact(graph, kDefaultTurns, true, 2U);

    if (std::fabs(single.expectation - threaded.expectation) > kToleranceLoose) {
        std::cerr << "Thread consistency mismatch: single=" << std::setprecision(18)
                  << single.expectation << ", threaded=" << threaded.expectation << '\n';
        return false;
    }

    return true;
}

bool run_graph_structure_checkpoint(const Graph& graph) {
    for (const State& state : graph.states) {
        if (!has_contiguous_labels(state)) {
            std::cerr << "Non-contiguous labels detected in canonical state.\n";
            return false;
        }
    }

    for (std::size_t sid = 0ULL; sid < graph.transitions.size(); ++sid) {
        int sum_w = 0;
        for (const Transition& tr : graph.transitions[sid]) {
            sum_w += tr.weight;
        }
        if (sum_w != graph.alphabet) {
            std::cerr << "Transition weight sum mismatch at state " << sid
                      << ": got " << sum_w << ", expected " << graph.alphabet << ".\n";
            return false;
        }
    }

    return true;
}

bool run_checkpoints(const Graph& graph) {
    if (!run_graph_structure_checkpoint(graph)) {
        std::cerr << "Graph-structure checkpoint failed.\n";
        return false;
    }

    if (!run_example_checkpoint()) {
        std::cerr << "Example checkpoint failed.\n";
        return false;
    }

    if (!run_bruteforce_cross_checkpoint()) {
        std::cerr << "Brute-force checkpoint failed.\n";
        return false;
    }

    if (!run_thread_consistency_checkpoint(graph)) {
        std::cerr << "Thread consistency checkpoint failed.\n";
        return false;
    }

    const ExactSolveResult main_exact = solve_exact(graph, kDefaultTurns, false, 1U);
    if (std::fabs(main_exact.expectation - kMainExpected) > kToleranceTight) {
        std::cerr << "Main-value checkpoint mismatch: got " << std::setprecision(18)
                  << main_exact.expectation << ", expected " << kMainExpected << '\n';
        return false;
    }

    if (std::fabs(main_exact.probability_mass - 1.0L) > kMassTolerance) {
        std::cerr << "Probability mass checkpoint failed for main case.\n";
        return false;
    }

    return true;
}

}  // namespace

int main(int argc, char** argv) {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);

    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }

    const Graph graph = build_graph(options.alphabet, options.capacity);

    if (options.run_checkpoints &&
        options.alphabet == kDefaultAlphabet &&
        options.capacity == kDefaultCapacity &&
        options.turns == kDefaultTurns) {
        if (!run_checkpoints(graph)) {
            return 1;
        }
    }

    const auto start = std::chrono::steady_clock::now();
    const ExactSolveResult result =
        solve_exact(graph, options.turns, options.allow_multithreading, options.requested_threads);
    const auto end = std::chrono::steady_clock::now();

    const std::chrono::duration<long double> elapsed = end - start;

    std::cout << std::fixed << std::setprecision(15)
              << "E(|L-R|) = " << result.expectation << '\n';
    std::cout << std::setprecision(8)
              << "Answer (8 d.p.): " << result.expectation << '\n';
    std::cout << std::setprecision(15)
              << "Probability mass: " << result.probability_mass << '\n';
    std::cout << "States in canonical graph: " << graph.states.size() << '\n';
    std::cout << "Elapsed: " << elapsed.count() << " seconds\n";

    return 0;
}

Python

from collections import deque

def apply_call_in_place(state, called_label, capacity):
    l_size, r_size, l_arr, r_arr = state
    l_arr = list(l_arr)
    r_arr = list(r_arr)
    
    larry_hit = False
    robin_hit = False
    delta = 0
    
    if called_label in l_arr:
        larry_hit = True
        delta += 1
        l_arr.remove(called_label)
        l_arr.append(called_label)
    else:
        if len(l_arr) < capacity:
            l_arr.append(called_label)
        else:
            l_arr.pop(0)
            l_arr.append(called_label)
            
    if called_label in r_arr:
        robin_hit = True
        delta -= 1
    else:
        if len(r_arr) < capacity:
            r_arr.append(called_label)
        else:
            r_arr.pop(0)
            r_arr.append(called_label)
            
    return (len(l_arr), len(r_arr), tuple(l_arr), tuple(r_arr)), delta

def canonicalize_state(state):
    l_size, r_size, l_arr, r_arr = state
    remap = {}
    next_id = 0
    
    new_l = []
    for x in l_arr:
        if x not in remap:
            remap[x] = next_id
            next_id += 1
        new_l.append(remap[x])
        
    new_r = []
    for x in r_arr:
        if x not in remap:
            remap[x] = next_id
            next_id += 1
        new_r.append(remap[x])
        
    return (l_size, r_size, tuple(new_l), tuple(new_r))

def build_graph(alphabet, capacity):
    states = []
    transitions = []
    id_by_state = {}
    
    initial = (0, 0, (), ())
    canonical_initial = canonicalize_state(initial)
    id_by_state[canonical_initial] = 0
    states.append(canonical_initial)
    transitions.append([])
    
    queue = deque([0])
    
    while queue:
        u = queue.popleft()
        state = states[u]
        l_size, r_size, l_arr, r_arr = state
        
        present = set(l_arr) | set(r_arr)
        known_labels = list(present)
        outside_count = alphabet - len(known_labels)
        
        out_transitions = {}
        
        for called in known_labels:
            next_state, delta = apply_call_in_place(state, called, capacity)
            can = canonicalize_state(next_state)
            
            if can not in id_by_state:
                id_by_state[can] = len(states)
                states.append(can)
                transitions.append([])
                queue.append(id_by_state[can])
                
            nxt_id = id_by_state[can]
            key = (nxt_id, delta)
            out_transitions[key] = out_transitions.get(key, 0) + 1
            
        if outside_count > 0:
            representative = -1
            for x in range(alphabet):
                if x not in present:
                    representative = x
                    break
                    
            next_state, delta = apply_call_in_place(state, representative, capacity)
            can = canonicalize_state(next_state)
            
            if can not in id_by_state:
                id_by_state[can] = len(states)
                states.append(can)
                transitions.append([])
                queue.append(id_by_state[can])
                
            nxt_id = id_by_state[can]
            key = (nxt_id, delta)
            out_transitions[key] = out_transitions.get(key, 0) + outside_count
            
        transitions[u] = [(nxt, delta, weight) for ((nxt, delta), weight) in out_transitions.items()]
        
    return states, transitions

def solve_exact(states, transitions, alphabet, turns):
    state_count = len(states)
    width = 2 * turns + 1
    offset = turns
    
    current = [0.0] * (state_count * width)
    next_probs = [0.0] * (state_count * width)
    
    current[offset] = 1.0
    
    for turn in range(1, turns + 1):
        for i in range(len(next_probs)):
            next_probs[i] = 0.0
            
        d_prev_min = offset - (turn - 1)
        d_prev_max = offset + (turn - 1)
        
        for sid in range(state_count):
            src_base = sid * width
            
            for nxt, delta, weight in transitions[sid]:
                p = weight / float(alphabet)
                dst_base = nxt * width
                
                for d in range(d_prev_min, d_prev_max + 1):
                    prob = current[src_base + d]
                    if prob == 0.0:
                        continue
                        
                    nd = d + delta
                    next_probs[dst_base + nd] += prob * p
                    
        current, next_probs = next_probs, current
        
    expectation = 0.0
    for sid in range(state_count):
        base = sid * width
        for i in range(width):
            p = current[base + i]
            if p > 0.0:
                diff = i - offset
                expectation += p * abs(diff)
                
    return expectation

def solve():
    alphabet = 10
    capacity = 5
    turns = 50
    states, transitions = build_graph(alphabet, capacity)
    ans = solve_exact(states, transitions, alphabet, turns)
    return "{:.15f}".format(ans)

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

Java

import java.util.*;

public class Euler298 {

    static final int kDefaultAlphabet = 10;
    static final int kDefaultCapacity = 5;
    static final int kDefaultTurns = 50;

    static class State {
        int lSize, rSize;
        int[] l = new int[kDefaultCapacity];
        int[] r = new int[kDefaultCapacity];

        @Override
        public boolean equals(Object o) {
            if (this == o)
                return true;
            if (o == null || getClass() != o.getClass())
                return false;
            State state = (State) o;
            if (lSize != state.lSize || rSize != state.rSize)
                return false;
            for (int i = 0; i < lSize; i++)
                if (l[i] != state.l[i])
                    return false;
            for (int i = 0; i < rSize; i++)
                if (r[i] != state.r[i])
                    return false;
            return true;
        }

        @Override
        public int hashCode() {
            int result = Objects.hash(lSize, rSize);
            for (int i = 0; i < lSize; i++)
                result = 31 * result + l[i];
            for (int i = 0; i < rSize; i++)
                result = 31 * result + r[i];
            return result;
        }

        State copy() {
            State c = new State();
            c.lSize = lSize;
            c.rSize = rSize;
            System.arraycopy(l, 0, c.l, 0, lSize);
            System.arraycopy(r, 0, c.r, 0, rSize);
            return c;
        }
    }

    static class CallResult {
        int delta;
    }

    static int findIndex(int[] arr, int size, int val) {
        for (int i = 0; i < size; i++)
            if (arr[i] == val)
                return i;
        return -1;
    }

    static CallResult applyCall(State state, int calledLabel, int capacity) {
        CallResult result = new CallResult();
        int idxL = findIndex(state.l, state.lSize, calledLabel);
        if (idxL >= 0) {
            result.delta++;
            for (int i = idxL; i + 1 < state.lSize; i++)
                state.l[i] = state.l[i + 1];
            state.l[state.lSize - 1] = calledLabel;
        } else {
            if (state.lSize < capacity) {
                state.l[state.lSize++] = calledLabel;
            } else {
                for (int i = 0; i + 1 < capacity; i++)
                    state.l[i] = state.l[i + 1];
                state.l[capacity - 1] = calledLabel;
            }
        }

        int idxR = findIndex(state.r, state.rSize, calledLabel);
        if (idxR >= 0) {
            result.delta--;
        } else {
            if (state.rSize < capacity) {
                state.r[state.rSize++] = calledLabel;
            } else {
                for (int i = 0; i + 1 < capacity; i++)
                    state.r[i] = state.r[i + 1];
                state.r[capacity - 1] = calledLabel;
            }
        }
        return result;
    }

    static void canonicalizeState(State state) {
        int[] remap = new int[20];
        Arrays.fill(remap, -1);
        int nextId = 0;

        for (int i = 0; i < state.lSize; i++) {
            if (remap[state.l[i]] < 0)
                remap[state.l[i]] = nextId++;
        }
        for (int i = 0; i < state.rSize; i++) {
            if (remap[state.r[i]] < 0)
                remap[state.r[i]] = nextId++;
        }
        for (int i = 0; i < state.lSize; i++)
            state.l[i] = remap[state.l[i]];
        for (int i = 0; i < state.rSize; i++)
            state.r[i] = remap[state.r[i]];
    }

    static class Transition {
        int nextState;
        int delta;
        int weight;

        Transition(int n, int d, int w) {
            nextState = n;
            delta = d;
            weight = w;
        }

        @Override
        public boolean equals(Object o) {
            if (this == o)
                return true;
            if (o == null || getClass() != o.getClass())
                return false;
            Transition that = (Transition) o;
            return nextState == that.nextState && delta == that.delta;
        }

        @Override
        public int hashCode() {
            return Objects.hash(nextState, delta);
        }
    }

    static class Graph {
        List<State> states = new ArrayList<>();
        List<List<Transition>> transitions = new ArrayList<>();
    }

    static Graph buildGraph(int alphabet, int capacity) {
        Graph graph = new Graph();
        Map<State, Integer> idMap = new HashMap<>();

        State initial = new State();
        canonicalizeState(initial);
        idMap.put(initial, 0);
        graph.states.add(initial);
        graph.transitions.add(new ArrayList<>());

        Queue<Integer> queue = new ArrayDeque<>();
        queue.add(0);

        while (!queue.isEmpty()) {
            int u = queue.poll();
            State state = graph.states.get(u);

            boolean[] present = new boolean[20];
            for (int i = 0; i < state.lSize; i++)
                present[state.l[i]] = true;
            for (int i = 0; i < state.rSize; i++)
                present[state.r[i]] = true;

            List<Integer> knownLabels = new ArrayList<>();
            for (int i = 0; i < alphabet; i++)
                if (present[i])
                    knownLabels.add(i);

            int outsideCount = alphabet - knownLabels.size();
            Map<Transition, Integer> outTransitions = new HashMap<>();

            for (int called : knownLabels) {
                State next = state.copy();
                CallResult res = applyCall(next, called, capacity);
                canonicalizeState(next);

                if (!idMap.containsKey(next)) {
                    int nxtId = graph.states.size();
                    idMap.put(next, nxtId);
                    graph.states.add(next);
                    graph.transitions.add(new ArrayList<>());
                    queue.add(nxtId);
                }
                int nxtId = idMap.get(next);
                Transition key = new Transition(nxtId, res.delta, 0);
                outTransitions.put(key, outTransitions.getOrDefault(key, 0) + 1);
            }

            if (outsideCount > 0) {
                int rep = -1;
                for (int x = 0; x < alphabet; x++) {
                    if (!present[x]) {
                        rep = x;
                        break;
                    }
                }
                State next = state.copy();
                CallResult res = applyCall(next, rep, capacity);
                canonicalizeState(next);

                if (!idMap.containsKey(next)) {
                    int nxtId = graph.states.size();
                    idMap.put(next, nxtId);
                    graph.states.add(next);
                    graph.transitions.add(new ArrayList<>());
                    queue.add(nxtId);
                }
                int nxtId = idMap.get(next);
                Transition key = new Transition(nxtId, res.delta, 0);
                outTransitions.put(key, outTransitions.getOrDefault(key, 0) + outsideCount);
            }

            List<Transition> finalTrans = new ArrayList<>();
            for (Map.Entry<Transition, Integer> e : outTransitions.entrySet()) {
                Transition tr = e.getKey();
                tr.weight = e.getValue();
                finalTrans.add(tr);
            }
            graph.transitions.set(u, finalTrans);
        }

        return graph;
    }

    static double solveExact(Graph graph, int alphabet, int turns) {
        int stateCount = graph.states.size();
        int width = 2 * turns + 1;
        int offset = turns;

        double[] current = new double[stateCount * width];
        double[] next = new double[stateCount * width];

        current[offset] = 1.0;

        for (int turn = 1; turn <= turns; turn++) {
            Arrays.fill(next, 0.0);
            int dPrevMin = offset - (turn - 1);
            int dPrevMax = offset + (turn - 1);

            for (int sid = 0; sid < stateCount; sid++) {
                int srcBase = sid * width;
                for (Transition tr : graph.transitions.get(sid)) {
                    double p = (double) tr.weight / alphabet;
                    int dstBase = tr.nextState * width;

                    for (int d = dPrevMin; d <= dPrevMax; d++) {
                        double prob = current[srcBase + d];
                        if (prob == 0.0)
                            continue;
                        next[dstBase + d + tr.delta] += prob * p;
                    }
                }
            }
            double[] temp = current;
            current = next;
            next = temp;
        }

        double expectation = 0.0;
        for (int sid = 0; sid < stateCount; sid++) {
            int base = sid * width;
            for (int i = 0; i < width; i++) {
                double p = current[base + i];
                if (p > 0.0) {
                    expectation += p * Math.abs(i - offset);
                }
            }
        }
        return expectation;
    }

    public static void main(String[] args) {
        Graph graph = buildGraph(kDefaultAlphabet, kDefaultCapacity);
        double ans = solveExact(graph, kDefaultAlphabet, kDefaultTurns);
        System.out.printf(Locale.US, "%.15f\n", ans);
    }
}