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

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. We only need:

$k$: the number of half-open pairs, and a binary flag telling us whether $s$ has already appeared once.

Define $E_0[k]$ as the expected remaining number of draws when $k$ ordinary pairs are half-open and $s$ has not been seen. Define $E_1[k]$ as the same expectation when $s$ has been seen exactly once. This is the exact meaning of the arrays e0 and e1 in the C++, Python, and Java solutions.

Step 3: First-Step Equation for $E_1[k]$

Assume we are in state $E_1[k]$. Then:

The draw $0$ keeps the state unchanged. For each of the $k$ half-open pairs, drawing the already seen side also keeps the state unchanged. So there are $k+1$ same-state outcomes.

Drawing $s$ again wins immediately. For each of the $k$ half-open pairs, drawing the missing complement also wins immediately. So there are $k+1$ winning outcomes.

The remaining $m-k$ pairs are unopened. Any one of their $2(m-k)$ numbers opens a new pair, so the process moves to $E_1[k+1]$.

First-step analysis therefore gives

$$E_1[k]=1+\frac{k+1}{N}E_1[k]+\frac{2(m-k)}{N}E_1[k+1].$$

Moving the same-state term to the left yields the recurrence implemented in the code:

$$\boxed{E_1[k]=\frac{N+2(m-k)E_1[k+1]}{N-(k+1)}}.$$

The denominator $N-(k+1)$ is exactly the variable den in the source files. It appears because we subtract the probability of staying in the same state from the left-hand side.

Step 4: First-Step Equation for $E_0[k]$

Now assume $s$ has not appeared yet. The same-state outcomes are still $0$ plus the already seen sides of the $k$ half-open pairs, so again there are $k+1$ of them.

Drawing the missing complement of one of the $k$ half-open pairs wins immediately. Drawing $s$ does not win yet; it changes the state from $E_0[k]$ to $E_1[k]$. Finally, any number from one of the $m-k$ unopened pairs moves the process to $E_0[k+1]$.

Hence

$$E_0[k]=1+\frac{k+1}{N}E_0[k]+\frac{1}{N}E_1[k]+\frac{2(m-k)}{N}E_0[k+1],$$

$$\boxed{E_0[k]=\frac{N+2(m-k)E_0[k+1]+E_1[k]}{N-(k+1)}}.$$

Step 5: Boundary Condition and Final Answer

When $k=m$, every ordinary pair is already half-open, so there is no transition to $k+1$. The recurrence becomes

$$E_1[m]=\frac{N}{N-(m+1)}=\frac{2m+2}{m+1}=2,$$

$$E_0[m]=\frac{N+E_1[m]}{N-(m+1)}=\frac{2m+4}{m+1}.$$

From there we sweep backward for $k=m-1,m-2,\dots,0$. The required expectation is

$$\boxed{E_0[0]}. $$

Checkpoint Example: $M=3$

This is the exact case hard-coded in the C++ verification. Here $N=4$, there is one ordinary pair $\{1,3\}$, the special value $2$, and the inert value $0$. The backward recurrence gives

$$E_1[1]=2,\qquad E_0[1]=3,$$

$$E_1[0]=\frac{4+2E_1[1]}{3}=\frac{8}{3},$$

$$E_0[0]=\frac{4+2E_0[1]+E_1[0]}{3}=\frac{38}{9}.$$

This matches the checkpoint in Euler371.cpp. For the real input $M=999$, the same recurrence with $m=499$ gives

$$E_0[0]\approx 40.66368097.$$

How the Code Works

All three solution files implement the same compressed dynamic program. They compute pairCount = (maxNumber - 1) / 2, set totalNumbers = 2 * pairCount + 2, allocate arrays e0 and e1, and iterate $k$ downward from pairCount to $0$.

The helper quantities are

$$\texttt{grow}=2(\texttt{pairCount}-k),\qquad \texttt{den}=\texttt{totalNumbers}-(k+1).$$

Here grow counts the numbers that belong to still unopened pairs and therefore increase $k$ by one. Then the code applies the two boxed recurrences above. The final formatted result is e0[0].

The C++ file also contains an explicit-state verifier for tiny instances. In that checker, every pair has three microstates: unseen, left side seen, or right side seen. Together with the flag for whether $s$ has been seen, this creates $2\cdot 3^m$ states. The program builds the linear system for those exact Markov expectations, solves it by Gaussian elimination, and confirms that the compressed recurrence agrees for $m=2$ and $m=3$.

Complexity Analysis

Let $m=(M-1)/2$. The implemented recurrence evaluates each $k$ once, so the running time is $O(m)$. The current implementations store two arrays of length $m+1$, hence $O(m)$ memory. The explicit-state verifier is exponential in $m$ and is only used for very small checkpoint cases.

Footnotes and References

Problem page: https://projecteuler.net/problem=371
First-step analysis: Wikipedia — First-step analysis
Markov chains and expected hitting times: Wikipedia — Markov chain
Dynamic programming overview: Wikipedia — Dynamic programming

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());
    }
}