Problem 371: Licence Plates

Project Euler Problem 371 Solution

EulerSolve provides an optimized solution for Project Euler Problem 371, Licence Plates, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary In the original problem we observe random three-digit licence plates from 000 to 999 . We stop as soon as two observed plates sum to $1000$. The local solvers generalize the setting to any odd maximum $M$, so the available values are $\{0,1,\dots,M\}$ and the target sum is $M+1$. For the Project Euler input $M=999$, the program computes the expected stopping time and prints it to eight decimal places. Mathematical Approach Step 1: Partition the Numbers Let $$N=M+1,\qquad m=\frac{M-1}{2},\qquad s=\frac{M+1}{2}.$$ Then the numbers split into three types. For each $1\le i\le m$ we have one complementary pair $$P_i=\{i,N-i\}.$$ There are $m$ such pairs, one self-complementary value $s$ because $s+s=N$, and the singleton $0$, whose complement $N$ lies outside the range. When $M=999$, this means $499$ ordinary pairs, the special value $500$, and the harmless value $000$. Step 2: Compress the State Space Before the game ends, an ordinary pair can only be in two relevant conditions: unopened, meaning neither side has appeared yet, or half-open, meaning exactly one side has appeared. A fully completed pair would already end the process, so it never appears inside the recurrence. Therefore the detailed identity of the opened pairs does not matter....

Detailed mathematical approach

Problem Summary

An Oregon licence plate carries three letters and then a three digit number, so its number is one of $1000$ strings, 000 to 999. Seth wins as soon as two plates seen on one trip carry numbers adding to $1000$: the statement offers MIC-012 with HAN-988, and RYU-500 with SET-500. Each plate carries any of the $1000$ numbers with equal probability, independently of the others, and the letters are irrelevant. We must find the expected number of plates he has to see, counting the one that completes the win, to eight decimal places. ^[1]

The whole solution in one sentence: the $1000$ numbers split into one that can never help, one that matches itself and $499$ complementary pairs, so a trip without a win is described by how many pairs are half used and whether $500$ has turned up, and the waiting times of those few positions satisfy linear equations that one backward pass evaluates in order.

Mathematical Approach

1. Three kinds of plate number

Write $T=1000$ for how many plate numbers exist; the target sum is the same $T$, so two numbers $x$ and $y$ win when $x+y=T$. Sorting them by their partners gives three kinds.

The dead number. $x=0$ would need the partner $1000$, which is not a three digit number, so a plate numbered 000 never contributes to a win.
The self-matching number. $x=500$ is its own partner, so two plates numbered 500 win, as the second example of the statement shows.
The complementary pairs. Every other $x$ has exactly one partner $T-x\ne x$, and these $998$ numbers form $m=499$ pairs $\{x,1000-x\}$ for $x=1,\dots,499$.

The count confirms the split: $1+1+2\cdot499=1000$. Among the $T^2=10^6$ ordered pairs (first plate, second plate) a win needs the two to be partners, which holds for the $2\cdot499=998$ pairs made of a number and its partner and for $(500,500)$. So $999$ of them win, and the chance that the second plate already ends the trip is $999/10^6=0.000999$.

2. What a trip without a win looks like

Let $S$ be the list of numbers seen so far, repetitions included. No win has occurred precisely when no two entries of $S$ add up to $T$, and by Section 1 that one condition says two things: $500$ occurs at most once in $S$, since two copies would already win, and for each of the $m$ pairs at most one member occurs, possibly many times. The number $0$ is unrestricted. Both directions are immediate, since a winning couple is either two copies of $500$ or the two members of one pair.

3. Only two pieces of information matter

Call a pair half used when exactly one of its two members has appeared. Neither the order of the plates nor the repetitions change which numbers can still win, so by Section 2 a win-free position is recorded completely by the seen member, if any, of each pair and a flag for $500$: $3^{m}$ choices times $2$. Even the identity of those members is unnecessary; only the number $k$ of half-used pairs and the flag $s\in\{0,1\}$ matter. To see it, count the $T$ equally likely numbers the next plate can carry.

$1$ number is $0$: nothing changes.
$k$ numbers are the seen members of the half-used pairs: nothing changes.
$k$ numbers are their missing partners: the trip ends in a win.
$2(m-k)$ numbers belong to untouched pairs: whichever appears, that pair becomes half used and $k$ grows by $1$.
$1$ number is $500$: if $s=0$ it sets $s=1$; if $s=1$ the trip ends in a win.

The five counts add up to $2+2m=T$, and none depends on which members were seen, on the order of the plates or on how often $0$ appeared. In a Markov chain, whose next state depends only on the current one, states may be grouped into blocks whenever the total probability of moving from a state into each block is the same for all states of its own block; the chain on the blocks then serves just as well. This is the lumping of Kemeny and Snell. ^[4] Here the blocks are the pairs $(k,s)$, so the $2\cdot3^{m}$ detailed positions collapse to $2(m+1)$.

4. The expectation equations

Let $X$ be the number of plates seen up to and including the first win, and $E_s[k]$ the expected number of further plates from the position $(k,s)$, so the quantity asked for is $E_0[0]$. Conditioning on the next plate and weighting each outcome by its probability $1/T$ is the law of total expectation; ^[2] a winning draw ends the trip and contributes no further wait, so it adds no term. For $0\le k\lt m$ the five counts of Section 3 give

$$\begin{aligned}E_1[k]&=1+\frac{k+1}{T}E_1[k]\\ &\quad+\frac{2(m-k)}{T}E_1[k+1],\\ E_0[k]&=1+\frac{k+1}{T}E_0[k]+\frac{1}{T}E_1[k]\\ &\quad+\frac{2(m-k)}{T}E_0[k+1].\end{aligned}$$

The leading $1$ is the plate we are about to see, and the factor $(k+1)/T$ collects the $1+k$ numbers that change nothing. Moving that term to the left multiplies $E_s[k]$ by $1-(k+1)/T$; multiplying through by $T$ and writing $d_k=T-k-1$ for the denominator gives the recurrences the programs use.

$$\begin{aligned}E_1[k]&=\frac{T+2(m-k)E_1[k+1]}{d_k},\\ E_0[k]&=\frac{T+2(m-k)E_0[k+1]}{d_k}\\ &\quad+\frac{E_1[k]}{d_k}.\end{aligned}$$

Since $d_k=2m+1-k$ is at least $m+1$ for every $k\le m$, no division is degenerate. Both lines read a value at the larger index $k+1$, and the second also reads $E_1[k]$; hence the sweep runs downwards and computes $E_1[k]$ first. It starts at $k=m$, where there is no $k+1$ term, as Section 5 shows.

5. The two values the sweep starts from

At $k=m$ every pair is half used, so $m-k=0$ and the recurrences close.

$$\begin{aligned}E_1[m]&=\frac{T}{d_m}=\frac{2m+2}{m+1}=2,\\ E_0[m]&=\frac{T+2}{m+1}.\end{aligned}$$

The first value needs no algebra: there the $m$ missing partners and the second $500$ all win, which is $m+1$ of the $T=2m+2$ numbers, exactly half, so the waiting time is geometric with success probability $\tfrac12$ and mean $2$. The second is $2+2/(m+1)$, which for $m=499$ equals $2.004$: four thousandths of a plate more, because a first $500$ only sets the flag instead of winning. A missing partner can of course still end the trip before any $500$ appears.

6. A worked example: one-digit plates

Shrink the puzzle to plates with a single digit, $0$ to $9$, and target sum $10$. Now $T=10$, the dead number is $0$, the self-matching number is $5$, and the $m=4$ pairs are $\{1,9\}$, $\{2,8\}$, $\{3,7\}$ and $\{4,6\}$. The sweep of Section 4, with $d_k=9-k$ and multiplier $2(4-k)$, gives these exact values.

$k=4$: $E_1[4]=10/5=2$ and $E_0[4]=12/5$.
$k=3$: $E_1[3]=(10+2\cdot2)/6=7/3$, then $E_0[3]=(10+\tfrac{24}{5}+\tfrac73)/6=\tfrac{257}{90}$.
$k=2$: $E_1[2]=58/21$ and $E_0[2]=7618/2205$.
$k=1$: $E_1[1]=93/28$ and $E_0[1]=100109/23520$.
$k=0$: $E_1[0]=256/63$ and $E_0[0]=424367/79380$.

A player watching one-digit plates therefore waits $424367/79380\approx5.34601915$ plates on average. The same sweep with $m=1$, that is plates $0$ to $3$ with target sum $4$, gives $E_1[1]=2$, $E_0[1]=3$, $E_1[0]=8/3$ and $E_0[0]=38/9=4.2222\dots$, the first checkpoint built into the C++ program.

7. Why the expectation is finite, and how large it is

The recurrences are legitimate only if the expectations are finite, and they are, for a crude reason. From any position, if the next plate is not 000, which has probability $999/1000$, its number has a partner among the plate numbers, and the plate after it wins whenever it carries that partner, with probability $1/1000$. Every block of two plates therefore ends the trip with probability at least $999/10^6$, so the remaining wait is dominated by twice a geometric variable with that success probability and $E_0[0]\le2\cdot10^6/999\lt2003$. Adjoining one winning state to the $2(m+1)$ win-free positions gives an absorbing chain, whose expected absorption times are the unique solution of the system above. ^[3]

The truth lies far below that bound, and the classical birthday problem says where. Each couple of plates wins with probability $(T-1)/T^{2}$, so the expected number of winning couples among $n$ plates is $\lambda_n=\binom n2(T-1)/T^{2}$. The couples are dependent, but treating their number as Poisson is the classical heuristic ^[7] and is accurate here: $e^{-\lambda_n}$ is $0.45876371$ at $n=40$ against the exact survival probability $0.45872783$. Replacing $\lambda_n$ by $n^{2}/(2T)$ and integrating over $n$ leaves the birthday scale $\sqrt{\pi T/2}$, ^[7] about $39.6333$ plates at $T=1000$ and $8.1$ per cent below the expectation at $T=100$, which rounds to $13.63434282$.

8. An independent closed form

A second route reaches the same number without recursion, and we use it as a check. For $X$ with values in $\{0,1,2,\dots\}$ the tail-sum formula ^[5] gives $\mathbb E[X]=\sum_{n\ge0}\Pr(X\gt n)$, and here $\Pr(X\gt n)=V(n)/T^{\,n}$, where $V(n)$ counts the sequences of $n$ plate numbers free of a winning couple. To count those, fix a set $A$ of $j$ pairs with one chosen member of each, and count the sequences over those $j$ members, the number $0$ and at most one $500$: there are $(j+1)^{n}$ with no $500$ and $n(j+1)^{n-1}$ with exactly one, whose position is free. Running over all $\binom{m}{j}2^{j}$ such choices counts a valid sequence once for every $A$ that contains all $t$ pairs the sequence touches and whose chosen members agree with it there, that is $\binom{m-t}{j-t}2^{\,j-t}$ times. When $m=2$, for instance, a sequence touching exactly one pair belongs to one choice with $j=1$, of weight $-1$, and to two choices with $j=2$, each of weight $+1$, so its total weight is $-1+2=1$. In general, writing $u=m-t$, the identity

$$\sum_{i=0}^{u}(-1)^{u-i}\binom{u}{i}2^{i}=(2-1)^{u}=1,$$

shows that the alternating combination of these counts, the pattern of inclusion and exclusion, ^[6] keeps every valid sequence exactly once. With $a_j=(-1)^{m-j}\binom{m}{j}2^{j}$ and $c=j+1$,

$$V(n)=\sum_{j=0}^{m}a_j\bigl(c^{\,n}+nc^{\,n-1}\bigr).$$

Since $c\le m+1\lt T$, the geometric series gives $\sum_n c^{n}/T^{n}=T/(T-c)$ and its derivative $\sum_n nc^{n-1}/T^{n}=T/(T-c)^{2}$, so summing over $n$ collapses the tail sum to a finite closed form.

$$\mathbb E[X]=\sum_{j=0}^{m}a_j\left(\frac{T}{T-c}+\frac{T}{(T-c)^{2}}\right)$$

For $m=1$ and $T=4$ the terms are $-\bigl(\tfrac43+\tfrac49\bigr)=-\tfrac{16}{9}$ and $2\bigl(\tfrac42+\tfrac44\bigr)=6$, summing to $38/9$, the value found by hand in Section 6. For $m=499$ the formula has $500$ terms whose coefficients $|a_j|=2^{j}\binom{499}{j}$ reach a $237$-digit maximum and alternate in sign, so summing it in double precision is hopeless; in exact rational arithmetic it reproduces the value of the sweep.

How the Code Works

All three files compute the expectation with the backward sweep of Section 4, and the C++ file runs its checkpoints first unless they are switched off. ^[10] From the largest plate number $M$ they set pair_count, or pairCount in Java, to $m=(M-1)/2$ and total_numbers to $2m+2=T$, allocate e0 and e1 of length $m+1$ and run $k$ downwards to $0$. In the loop den is $d_k=T-k-1$ and grow is $2(m-k)$, the count of numbers from untouched pairs; the branch for $k=m$ uses the closed values of Section 5, the other the two recurrences, e1[k] before e0[k] because the second formula needs it. The printed result is e0[0] rounded to eight decimal places.

The C++ file adds two things. parse_arguments accepts --max-number= followed by digits and the flag --skip-checkpoints, rejecting an even $M$ or an $M\lt3$, since the three kinds of Section 1 need an odd $M\ge3$. Those checks do not by themselves fix a safe range: the digits are accumulated into an int without an overflow test, and $2m+2$ overflows for the largest odd int. run_checkpoints compares expected_steps_compressed(1) with $38/9$ to within $10^{-15}$, and the sweep at $m=2$ and $m=3$ with expected_steps_explicit_states to within $10^{-12}$.

That last function is the honest, unlumped chain of Section 3. It encodes a position as a base-$3$ numeral with one digit per pair — $0$ untouched, $1$ lower member seen, $2$ upper member seen — shifted up by one bit for the flag of the self-matching number, giving $2\cdot3^{m}$ states, and writes one row of $(I-Q)t=\mathbf 1$ per state: the diagonal entry starts at $1$ and loses $1/T$ for the number $0$ and $1/T$ per already seen member, each untouched pair adds two off-diagonal entries, the self-matching number adds one more while $s=0$, and a winning draw adds nothing, since after absorption no time remains. solve_linear_system applies Gauss–Jordan elimination with partial pivoting. ^[8] A default run therefore calls this dense solver at $m=2$ and $m=3$ before the sweep; it checks the lumping of Section 3 rather than contributing to the result.

The languages differ only in arithmetic and formatting. The precision of the C++ long double is left to the platform: on Apple ARM64 it is exactly double, elsewhere it may carry $64$ or $113$ significand bits. Python's float and Java's double are both normally binary64, and Java formats through Locale.US so the separator is a point. One asymmetry: C++ and Python set the largest plate number to $999$, while the Java solve() passes $1000$, the target sum; integer division sends both to $m=499$, so the output agrees, but the arguments mean different things and $1000$ would fail the C++ check.

Complexity and Verification

The sweep does a constant amount of arithmetic at each of the $m+1$ indices, so the method costs $O(m)$ time and $O(m)$ memory: at $m=499$ that is $500$ iterations and two sequences of $500$ values, an $8000$-byte payload in the Java double[] arrays and more in the other two, whose elements are wider or are separate objects. The whole C++ program, checkpoints included, takes a couple of milliseconds. The unlumped verifier is the opposite, with $n=2\cdot3^{m}$ unknowns and an augmented matrix of $n(n+1)$ entries; since Gauss–Jordan elimination costs $O(n^{3})$ operations, ^[8] that is about $54^{3}\approx1.6\cdot10^{5}$ operations at $m=3$. At $m=499$ the state count would be a $239$-digit number, so the lumping of Section 3 is not a refinement but the reason the problem is solvable.

Precision deserves a word, since eight decimals must survive a chain of $1000$ divisions. Repeating the sweep in exact rational arithmetic gives a fraction whose denominator has $727$ digits; the double result differs from it by $2.8\cdot10^{-15}$, less than half the gap of about $7.1\cdot10^{-15}$ between neighbouring doubles there. The sweep is therefore correctly rounded, and a double, whose $53$-bit significand is worth $15$ to $17$ decimal digits, ^[9] keeps some seven digits in reserve.

The checks below were run for this article; all but the C++ checkpoints are additional validation, not packaged in the published files.

The three programs print the same text, and it is the exact rational value of $E_0[0]$ rounded to eight decimal places. With largest plate number $3$, $5$, $7$, $9$, $11$ the C++ binary prints $4.22222222$, $4.60666667$, $4.98721088$, $5.34601915$ and $5.68263089$, the eight-place roundings of $38/9$, $691/150$, $18328/3675$, $424367/79380$ and $4548463/800415$.
The unlumped chain agrees with the sweep for $m=1,\dots,5$: exactly, in rational arithmetic, for $m\le3$, and to $1.8\cdot10^{-14}$ at $m=5$, where it has $486$ unknowns.
In exact rational arithmetic the closed form of Section 8 equals the sweep for $m=1,\dots,6$ and for $m=499$; its counts $V(n)$ match an exhaustive enumeration of all $6^{n}$ sequences when $T=6$, giving $1,6,31,140,569,2150,7731$ for $n=0,\dots,6$; and the tail sum $\sum_{n\ge0}\Pr(X\gt n)$, with survival probabilities pushed forward plate by plate, reproduces the sweep when truncated at $n=324$, where the surviving probability is $9.4\cdot10^{-23}$ and the bound of Section 7 keeps the discarded remainder below $2\cdot10^{-19}$.
A simulation of $400\,000$ trips agrees to within one standard error. The sample standard deviation was about $20.7$ plates, and the standard error of a mean falls only as $1/\sqrt N$, so eight decimal places are out of reach of any practical sample.

Footnotes and References

The references supply the background facts used above; the classification of the numbers, the lumped state $(k,s)$, the recurrences and the closed form of Section 8 are derived here.

Project Euler 371 — Licence Plates. The official statement describes plates of three letters and a three digit number, declares a win when two plates seen on one trip have numbers adding to 1000, gives MIC-012 with HAN-988 and RYU-500 with SET-500 as examples, and asks for the expected number of plates needed, to eight decimal places, each plate being equally likely to carry any such number.
Wikipedia — Law of total expectation. For a countable partition {Aᵢ}, E[X] = Σᵢ E[X | Aᵢ] Pr(Aᵢ); Section 4 splits the waiting time this way over the 1000 numbers the next plate can carry.
Wikipedia — Absorbing Markov chain. The expected steps before absorption form t = N1 with N = (I − Q)⁻¹, Q the transition matrix among transient states, so the system of Section 4 has exactly one solution.
Wikipedia — Lumpability. After Kemeny and Snell, a chain is lumpable for a partition exactly when any two states of one block have the same summed transition probability into each block; Section 3 verifies this for (k, s).
Wikipedia — Expected value. The article states the tail-sum formula E[X] = Σ Pr(X > k), summed from k = 0, for values in {0, 1, 2, …}, the form Section 8 uses.
Wikipedia — Inclusion–exclusion principle. Over-generous inclusion followed by compensating exclusion counts each element once, with the sign (−1)^(|J|+1) on each nonempty index set J; Section 8 uses that pattern with the weights 2^j.
Wikipedia — Birthday problem. The article approximates the probability of no shared birthday among n people from d days by exp(−n(n − 1)/(2d)), and gives √(πd/2) as the leading term of the expected number of people needed for the first repeat.
Wikipedia — Gaussian elimination. Row reduction solves a linear system in O(n³) arithmetic operations, and taking the largest available absolute value as the pivot improves stability.
Wikipedia — Double-precision floating-point format. The binary64 format has a 53-bit significand, 52 bits of it stored, worth 15 to 17 significant decimal digits.
C++ — Euler371.cpp. The linked immutable revision contains expected_steps_compressed, expected_steps_explicit_states, solve_linear_system, run_checkpoints and the --max-number= option; it is the implementation record for this article.

Problem 371 source code

C++

#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <string>
#include <vector>

namespace {

struct Options {
    int max_number = 999;
    bool run_checkpoints = true;
};

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }
    int parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<int>(c - '0');
    }
    value = parsed;
    return true;
}

bool parse_arguments(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 (parse_int_after_prefix(arg, "--max-number=", options.max_number)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.max_number >= 3 && (options.max_number % 2 == 1);
}

long double expected_steps_compressed(const int pair_count) {
    const int total_numbers = 2 * pair_count + 2;
    std::vector<long double> e0(static_cast<std::size_t>(pair_count + 1), 0.0L);
    std::vector<long double> e1(static_cast<std::size_t>(pair_count + 1), 0.0L);

    for (int k = pair_count; k >= 0; --k) {
        const long double den = static_cast<long double>(total_numbers - (k + 1));
        const long double grow = static_cast<long double>(2 * (pair_count - k));
        if (k == pair_count) {
            e1[static_cast<std::size_t>(k)] = static_cast<long double>(total_numbers) / den;
            e0[static_cast<std::size_t>(k)] =
                (static_cast<long double>(total_numbers) + e1[static_cast<std::size_t>(k)]) / den;
        } else {
            e1[static_cast<std::size_t>(k)] =
                (static_cast<long double>(total_numbers) +
                 grow * e1[static_cast<std::size_t>(k + 1)]) /
                den;
            e0[static_cast<std::size_t>(k)] =
                (static_cast<long double>(total_numbers) +
                 grow * e0[static_cast<std::size_t>(k + 1)] +
                 e1[static_cast<std::size_t>(k)]) /
                den;
        }
    }
    return e0[0];
}

std::vector<long double> solve_linear_system(std::vector<std::vector<long double>> matrix) {
    const int n = static_cast<int>(matrix.size());
    for (int col = 0; col < n; ++col) {
        int pivot = col;
        for (int row = col + 1; row < n; ++row) {
            if (std::fabsl(matrix[static_cast<std::size_t>(row)][static_cast<std::size_t>(col)]) >
                std::fabsl(matrix[static_cast<std::size_t>(pivot)][static_cast<std::size_t>(col)])) {
                pivot = row;
            }
        }
        std::swap(matrix[static_cast<std::size_t>(col)], matrix[static_cast<std::size_t>(pivot)]);
        const long double diag = matrix[static_cast<std::size_t>(col)][static_cast<std::size_t>(col)];
        for (int j = col; j <= n; ++j) {
            matrix[static_cast<std::size_t>(col)][static_cast<std::size_t>(j)] /= diag;
        }
        for (int row = 0; row < n; ++row) {
            if (row == col) {
                continue;
            }
            const long double factor = matrix[static_cast<std::size_t>(row)][static_cast<std::size_t>(col)];
            if (factor == 0.0L) {
                continue;
            }
            for (int j = col; j <= n; ++j) {
                matrix[static_cast<std::size_t>(row)][static_cast<std::size_t>(j)] -=
                    factor * matrix[static_cast<std::size_t>(col)][static_cast<std::size_t>(j)];
            }
        }
    }

    std::vector<long double> solution(static_cast<std::size_t>(n), 0.0L);
    for (int i = 0; i < n; ++i) {
        solution[static_cast<std::size_t>(i)] = matrix[static_cast<std::size_t>(i)][static_cast<std::size_t>(n)];
    }
    return solution;
}

long double expected_steps_explicit_states(const int pair_count) {
    const int total_numbers = 2 * pair_count + 2;
    std::vector<int> pow3(static_cast<std::size_t>(pair_count + 1), 1);
    for (int i = 1; i <= pair_count; ++i) {
        pow3[static_cast<std::size_t>(i)] = 3 * pow3[static_cast<std::size_t>(i - 1)];
    }

    const int states = 2 * pow3[static_cast<std::size_t>(pair_count)];
    std::vector<std::vector<long double>> matrix(
        static_cast<std::size_t>(states),
        std::vector<long double>(static_cast<std::size_t>(states + 1), 0.0L));

    const long double prob = 1.0L / static_cast<long double>(total_numbers);

    for (int idx = 0; idx < states; ++idx) {
        const int self_seen = idx & 1;
        const int code = idx >> 1;

        matrix[static_cast<std::size_t>(idx)][static_cast<std::size_t>(idx)] = 1.0L;
        matrix[static_cast<std::size_t>(idx)][static_cast<std::size_t>(states)] = 1.0L;

        matrix[static_cast<std::size_t>(idx)][static_cast<std::size_t>(idx)] -= prob;  // Draw 0.

        if (self_seen == 0) {
            const int next_idx = (code << 1) | 1;
            matrix[static_cast<std::size_t>(idx)][static_cast<std::size_t>(next_idx)] -= prob;
        }
        // If self_seen == 1 and we draw self, the game ends (absorbing win), so no matrix term.

        for (int p = 0; p < pair_count; ++p) {
            const int place = pow3[static_cast<std::size_t>(p)];
            const int state = (code / place) % 3;

            if (state == 0) {
                const int left_code = code + place;      // 0 -> 1
                const int right_code = code + 2 * place; // 0 -> 2
                matrix[static_cast<std::size_t>(idx)][static_cast<std::size_t>((left_code << 1) | self_seen)] -=
                    prob;
                matrix[static_cast<std::size_t>(idx)][static_cast<std::size_t>((right_code << 1) | self_seen)] -=
                    prob;
            } else if (state == 1) {
                // Draw left: no change. Draw right: win.
                matrix[static_cast<std::size_t>(idx)][static_cast<std::size_t>(idx)] -= prob;
            } else {
                // Draw right: no change. Draw left: win.
                matrix[static_cast<std::size_t>(idx)][static_cast<std::size_t>(idx)] -= prob;
            }
        }
    }

    const std::vector<long double> solution = solve_linear_system(std::move(matrix));
    const int initial_state = 0;  // No pair-side seen, self not seen.
    return solution[static_cast<std::size_t>(initial_state)];
}

long double solve(const int max_number) {
    const int pair_count = (max_number - 1) / 2;
    return expected_steps_compressed(pair_count);
}

bool run_checkpoints() {
    const long double expected_m1 = 38.0L / 9.0L;
    if (std::fabsl(expected_steps_compressed(1) - expected_m1) > 1e-15L) {
        std::cerr << "Checkpoint failed for pair_count=1 closed form" << '\n';
        return false;
    }

    const long double m2_fast = expected_steps_compressed(2);
    const long double m2_exact = expected_steps_explicit_states(2);
    if (std::fabsl(m2_fast - m2_exact) > 1e-12L) {
        std::cerr << "Checkpoint failed for explicit-state cross-check at pair_count=2" << '\n';
        return false;
    }

    const long double m3_fast = expected_steps_compressed(3);
    const long double m3_exact = expected_steps_explicit_states(3);
    if (std::fabsl(m3_fast - m3_exact) > 1e-12L) {
        std::cerr << "Checkpoint failed for explicit-state cross-check at pair_count=3" << '\n';
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }
    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }
    std::cout << std::fixed << std::setprecision(8) << solve(options.max_number) << '\n';
    return 0;
}

Python

def solve():
    max_number = 999
    pair_count = (max_number - 1) // 2
    total_numbers = 2 * pair_count + 2

    e0 = [0.0] * (pair_count + 1)
    e1 = [0.0] * (pair_count + 1)

    for k in range(pair_count, -1, -1):
        den = float(total_numbers - (k + 1))
        grow = float(2 * (pair_count - k))
        if k == pair_count:
            e1[k] = float(total_numbers) / den
            e0[k] = (float(total_numbers) + e1[k]) / den
        else:
            e1[k] = (float(total_numbers) + grow * e1[k + 1]) / den
            e0[k] = (float(total_numbers) + grow * e0[k + 1] + e1[k]) / den

    result = e0[0]
    return f"{result:.8f}"

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

Java

public class Euler371 {
    public static String solve(int maxNumber) {
        int pairCount = (maxNumber - 1) / 2;
        int totalNumbers = 2 * pairCount + 2;

        double[] e0 = new double[pairCount + 1];
        double[] e1 = new double[pairCount + 1];

        for (int k = pairCount; k >= 0; --k) {
            double den = totalNumbers - (k + 1);
            double grow = 2 * (pairCount - k);

            if (k == pairCount) {
                e1[k] = totalNumbers / den;
                e0[k] = (totalNumbers + e1[k]) / den;
            } else {
                e1[k] = (totalNumbers + grow * e1[k + 1]) / den;
                e0[k] = (totalNumbers + grow * e0[k + 1] + e1[k]) / den;
            }
        }

        return String.format(java.util.Locale.US, "%.8f", e0[0]);
    }

    public static String solve() {
        return solve(1000);
    }

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