Problem 232: The Race

Project Euler Problem 232 Solution

EulerSolve provides an optimized solution for Project Euler Problem 232, The Race, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Player 1 and Player 2 both start at 0 and race to 100 points. Player 1 always plays the simple strategy: on each turn, flip one fair coin and score 1 point on heads, 0 on tails. Player 2 has a much richer move: before each turn, choose a positive integer $T$, then succeed with probability $2^{-T}$; on success, score $2^{T-1}$ points, otherwise score nothing. Player 1 moves first, and the question asks for Player 2's maximum possible winning probability under optimal choices of $T$. The solutions do not simulate individual games. Instead they solve a finite optimal-control problem on the remaining distances to the target. For the actual target 100, the computed starting probability is $0.83648556$ to eight decimal places. Mathematical Approach The right viewpoint is to forget absolute scores and track only how many points each player still needs. That turns the game into a dynamic program with an explicit Bellman recurrence. Two state functions for the two turn types Let $F(a,b)$ be the optimal probability that Player 2 eventually wins when Player 1 still needs $a$ points, Player 2 still needs $b$ points, and it is Player 2's turn to move. It is also convenient to introduce $G(a,b)$, the corresponding probability at the start of Player 1's turn with the same remaining distances....

Detailed mathematical approach

Problem Summary

Player 1 and Player 2 both start at 0 and race to 100 points. Player 1 always plays the simple strategy: on each turn, flip one fair coin and score 1 point on heads, 0 on tails. Player 2 has a much richer move: before each turn, choose a positive integer $T$, then succeed with probability $2^{-T}$; on success, score $2^{T-1}$ points, otherwise score nothing. Player 1 moves first, and the question asks for Player 2's maximum possible winning probability under optimal choices of $T$.

The solutions do not simulate individual games. Instead they solve a finite optimal-control problem on the remaining distances to the target. For the actual target 100, the computed starting probability is $0.83648556$ to eight decimal places.

Mathematical Approach

The right viewpoint is to forget absolute scores and track only how many points each player still needs. That turns the game into a dynamic program with an explicit Bellman recurrence.

Two state functions for the two turn types

Let $F(a,b)$ be the optimal probability that Player 2 eventually wins when Player 1 still needs $a$ points, Player 2 still needs $b$ points, and it is Player 2's turn to move.

It is also convenient to introduce $G(a,b)$, the corresponding probability at the start of Player 1's turn with the same remaining distances. The initial position of the real problem is therefore $G(100,100)$, not $F(100,100)$, because Player 1 moves first.

The boundary conditions are immediate:

$F(a,0)=1$ and $G(a,0)=1$, because Player 2 has already reached the target; $G(0,b)=0$, because if Player 1 has already arrived before Player 2 gets another move, Player 2 has lost.

The easy transition: Player 1's turn

Player 1 either scores 1 point with probability $1/2$ or scores nothing with probability $1/2$. Therefore, for $a\ge 1$ and $b\ge 1$,

$$G(a,b)=\frac{F(a,b)+F(a-1,b)}{2}.$$

This formula is the first invariant used by the code: every time the game passes through Player 1's turn, the continuation value is just the average of two neighboring $F$-states.

Player 2's move and the raw Bellman equation

Suppose Player 2 chooses an integer $t\ge 1$. Then

$$q_t=2^{-t},\qquad d_t=2^{t-1}$$

are, respectively, the success probability and the number of points gained on success.

Define the success continuation value by

$S_t(a,b)=1$ if $d_t\ge b$, because Player 2 wins immediately, and otherwise $S_t(a,b)=G(a,b-d_t)$, because after scoring $d_t$ points the turn passes to Player 1.

If $t$ is fixed, the one-step Bellman equation is therefore

$$F(a,b)=q_t\,S_t(a,b)+(1-q_t)\,G(a,b).$$

The subtle point is the failure branch: when Player 2 misses, the scores do not change, so the game returns to the start of Player 1's turn at the same pair $(a,b)$. That is why $F(a,b)$ appears on both sides after we expand $G(a,b)$.

Removing the self-reference

Substitute $G(a,b)=\frac{F(a,b)+F(a-1,b)}{2}$ into the previous equation:

$$F(a,b)=q_t\,S_t(a,b)+(1-q_t)\frac{F(a,b)+F(a-1,b)}{2}.$$

Now solve this linear equation for $F(a,b)$. After collecting terms,

$$F(a,b)=\frac{2q_t\,S_t(a,b)+(1-q_t)F(a-1,b)}{1+q_t}.$$

Since Player 2 chooses the best $t$, the true recurrence is

$$\boxed{F(a,b)=\max_{t\ge 1}\frac{2q_t\,S_t(a,b)+(1-q_t)F(a-1,b)}{1+q_t}.}$$

This closed form is exactly what the implementations evaluate for every state $(a,b)$.

Why the dynamic program can be filled bottom-up

For a fixed state $(a,b)$, the formula only refers to two kinds of previously known values:

$F(a-1,b)$, which lies in the previous row, and $S_t(a,b)$, which either equals 1 or involves $G(a,b-d_t)$ with $b-d_t<b$.

That means a row-by-row, left-to-right fill order works. When the code reaches $(a,b)$, every state on which it depends has already been computed.

The action set is also finite. Once $d_t\ge b$, success already means immediate victory, and any larger $t$ only lowers the success probability $q_t$. So only $O(\log b)$ candidates can matter, and a small global cap is enough for the whole target 100 table.

Worked example: why a bigger gamble can be optimal

Consider the state $(a,b)=(1,2)$: Player 1 needs 1 more point, Player 2 needs 2 more points, and it is Player 2's turn.

First compute the simpler state $(1,1)$. If Player 2 chooses $t=1$, then $q_1=1/2$ and success ends the game immediately. The formula gives

$$F(1,1)=\frac{2\cdot(1/2)\cdot 1+(1-1/2)\cdot 0}{1+1/2}=\frac{2}{3},$$

$$G(1,1)=\frac{F(1,1)+F(0,1)}{2}=\frac{2/3+0}{2}=\frac{1}{3}.$$

Now compare two actions at $(1,2)$.

If Player 2 chooses $t=1$, a success scores only 1 point, so the success branch leads to $G(1,1)=1/3$. Since $F(0,2)=0$,

$$F_{t=1}(1,2)=\frac{2\cdot(1/2)\cdot(1/3)+(1-1/2)\cdot 0}{1+1/2}=\frac{2}{9}\approx 0.2222.$$

If Player 2 chooses $t=2$, then $q_2=1/4$ but the gain is 2 points, which wins immediately:

$$F_{t=2}(1,2)=\frac{2\cdot(1/4)\cdot 1+(1-1/4)\cdot 0}{1+1/4}=\frac{2}{5}=0.4.$$

So the optimal move at $(1,2)$ is the riskier $t=2$. This is exactly the kind of local tradeoff the recurrence captures across the whole $100\times 100$ state space.

How the Code Works

Table setup and precomputation

The C++, Python, and Java implementations store a two-dimensional table for $F(a,b)$. The column $b=0$ is initialized to 1, expressing that Player 2 has already won whenever no further points are needed.

They also precompute the candidate actions: for each tested $t$, the code stores $q_t=2^{-t}$ and $d_t=2^{t-1}$. Because the gains double each time, only a small number of actions must be examined for target 100.

The bottom-up sweep

The main loop visits states in increasing order of $a$ and $b$. At each cell, the implementation reads the already known value $F(a-1,b)$, evaluates every admissible action $t$, computes the corresponding success value $S_t(a,b)$, plugs those quantities into the closed-form recurrence, and keeps the maximum.

The C++ implementation additionally stores the maximizing choice of $t$ for each state, while the Python and Java implementations keep only the probability table because the final answer needs only the optimal value.

Extracting the start state and cross-checking

After the table is complete, the required answer is obtained by one final Player 1 transition:

$$G(100,100)=\frac{F(100,100)+F(99,100)}{2}.$$

The C++ program also includes an independent Bellman value-iteration verifier on smaller targets and checks that its threaded verifier agrees with the single-threaded version. Those checks are not part of the mathematical core, but they confirm that the closed-form recurrence has been implemented correctly.

Complexity Analysis

For a general target $N$, there are $N^2$ states $(a,b)$. Each state tests only $L$ candidate actions, where $L=O(\log N)$ because the possible gains are powers of two. The running time is therefore $O(N^2L)=O(N^2\log N)$.

The memory usage is $O(N^2)$ for the probability table. One implementation also stores an auxiliary policy table and a verification routine, but the main solver itself is still a straightforward $O(N^2)$-space dynamic program. For the actual Project Euler input $N=100$, this is easily fast enough.

Footnotes and References

Problem page: https://projecteuler.net/problem=232
Dynamic programming: Wikipedia - Dynamic programming
Bellman equation: Wikipedia - Bellman equation
Markov decision process: Wikipedia - Markov decision process
Bernoulli trial: Wikipedia - Bernoulli trial

Problem 232 source code

C++

#include <algorithm>
#include <atomic>
#include <cmath>
#include <cstdlib>
#include <iomanip>
#include <iostream>
#include <thread>
#include <vector>

namespace {

struct SolveResult {
    std::vector<std::vector<double>> f;      // P2 win prob at start of P2 turn
    std::vector<std::vector<int>> best_t;    // Optimal T at each (a,b)
    double start_probability = 0.0;          // P2 win prob from initial state (P1 starts)
};

unsigned compute_t_limit(int target) {
    unsigned t = 1;
    long long score = 1;
    while (score < target && t < 60) {
        ++t;
        score <<= 1;
    }
    return t + 4;
}

inline double p1_turn_value(const std::vector<std::vector<double>>& f, int a, int b) {
    if (a <= 0) return 0.0;
    if (b <= 0) return 1.0;
    const double prev_row = (a == 1) ? 0.0 : f[static_cast<std::size_t>(a - 1)][static_cast<std::size_t>(b)];
    return 0.5 * (f[static_cast<std::size_t>(a)][static_cast<std::size_t>(b)] + prev_row);
}

SolveResult solve_closed_form(int target) {
    SolveResult out;
    out.f.assign(static_cast<std::size_t>(target + 1), std::vector<double>(static_cast<std::size_t>(target + 1), 0.0));
    out.best_t.assign(static_cast<std::size_t>(target + 1), std::vector<int>(static_cast<std::size_t>(target + 1), 1));

    for (int a = 1; a <= target; ++a) {
        out.f[static_cast<std::size_t>(a)][0] = 1.0;
    }

    const unsigned t_limit = compute_t_limit(target);
    std::vector<double> q(static_cast<std::size_t>(t_limit + 1), 0.0);
    std::vector<int> score(static_cast<std::size_t>(t_limit + 1), 0);
    for (unsigned t = 1; t <= t_limit; ++t) {
        q[static_cast<std::size_t>(t)] = std::ldexp(1.0, -static_cast<int>(t));
        score[static_cast<std::size_t>(t)] = 1 << (t - 1);
    }

    for (int a = 1; a <= target; ++a) {
        for (int b = 1; b <= target; ++b) {
            const double prev_same_b = (a == 1) ? 0.0 : out.f[static_cast<std::size_t>(a - 1)][static_cast<std::size_t>(b)];

            double best_value = 0.0;
            int best_action = 1;

            for (unsigned t = 1; t <= t_limit; ++t) {
                const int gain = score[static_cast<std::size_t>(t)];
                const double success_value = (gain >= b)
                    ? 1.0
                    : p1_turn_value(out.f, a, b - gain);

                const double qt = q[static_cast<std::size_t>(t)];
                const double candidate = (2.0 * qt * success_value + (1.0 - qt) * prev_same_b) / (1.0 + qt);

                if (candidate > best_value) {
                    best_value = candidate;
                    best_action = static_cast<int>(t);
                }
            }

            out.f[static_cast<std::size_t>(a)][static_cast<std::size_t>(b)] = best_value;
            out.best_t[static_cast<std::size_t>(a)][static_cast<std::size_t>(b)] = best_action;
        }
    }

    out.start_probability = p1_turn_value(out.f, target, target);
    return out;
}

template <typename RowFn>
void parallel_rows(int row_begin, int row_end, unsigned threads, RowFn&& fn) {
    if (row_begin > row_end) return;
    if (threads <= 1) {
        for (int row = row_begin; row <= row_end; ++row) fn(row);
        return;
    }

    std::atomic<int> next_row{row_begin};
    std::vector<std::thread> workers;
    workers.reserve(static_cast<std::size_t>(threads));

    for (unsigned tid = 0; tid < threads; ++tid) {
        workers.emplace_back([&]() {
            while (true) {
                const int row = next_row.fetch_add(1, std::memory_order_relaxed);
                if (row > row_end) break;
                fn(row);
            }
        });
    }

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

// Independent Bellman value-iteration verifier.
double solve_value_iteration(int target, unsigned threads, double tolerance = 1e-14, int max_iterations = 5000) {
    if (target <= 0) return 0.0;

    const unsigned t_limit = compute_t_limit(target);
    std::vector<double> q(static_cast<std::size_t>(t_limit + 1), 0.0);
    std::vector<int> score(static_cast<std::size_t>(t_limit + 1), 0);
    for (unsigned t = 1; t <= t_limit; ++t) {
        q[static_cast<std::size_t>(t)] = std::ldexp(1.0, -static_cast<int>(t));
        score[static_cast<std::size_t>(t)] = 1 << (t - 1);
    }

    std::vector<std::vector<double>> p1(static_cast<std::size_t>(target + 1), std::vector<double>(static_cast<std::size_t>(target + 1), 0.0));
    std::vector<std::vector<double>> p2(static_cast<std::size_t>(target + 1), std::vector<double>(static_cast<std::size_t>(target + 1), 0.0));
    std::vector<std::vector<double>> next_p1(static_cast<std::size_t>(target + 1), std::vector<double>(static_cast<std::size_t>(target + 1), 0.0));
    std::vector<std::vector<double>> next_p2(static_cast<std::size_t>(target + 1), std::vector<double>(static_cast<std::size_t>(target + 1), 0.0));

    for (int a = 0; a <= target; ++a) {
        p1[static_cast<std::size_t>(a)][0] = 1.0;
        p2[static_cast<std::size_t>(a)][0] = 1.0;
        next_p1[static_cast<std::size_t>(a)][0] = 1.0;
        next_p2[static_cast<std::size_t>(a)][0] = 1.0;
    }

    std::vector<double> row_delta(static_cast<std::size_t>(target + 1), 0.0);

    for (int iteration = 0; iteration < max_iterations; ++iteration) {
        parallel_rows(1, target, threads, [&](int a) {
            for (int b = 1; b <= target; ++b) {
                const double p1_if_head = (a == 1) ? 0.0 : p2[static_cast<std::size_t>(a - 1)][static_cast<std::size_t>(b)];
                const double p1_if_tail = p2[static_cast<std::size_t>(a)][static_cast<std::size_t>(b)];
                next_p1[static_cast<std::size_t>(a)][static_cast<std::size_t>(b)] = 0.5 * (p1_if_head + p1_if_tail);
            }
        });

        parallel_rows(1, target, threads, [&](int a) {
            double local_max = 0.0;
            for (int b = 1; b <= target; ++b) {
                const double stay = next_p1[static_cast<std::size_t>(a)][static_cast<std::size_t>(b)];
                double best_value = 0.0;

                for (unsigned t = 1; t <= t_limit; ++t) {
                    const int gain = score[static_cast<std::size_t>(t)];
                    const double success_value = (gain >= b)
                        ? 1.0
                        : next_p1[static_cast<std::size_t>(a)][static_cast<std::size_t>(b - gain)];
                    const double qt = q[static_cast<std::size_t>(t)];
                    const double candidate = qt * success_value + (1.0 - qt) * stay;
                    if (candidate > best_value) best_value = candidate;
                }

                next_p2[static_cast<std::size_t>(a)][static_cast<std::size_t>(b)] = best_value;

                const double d1 = std::fabs(next_p1[static_cast<std::size_t>(a)][static_cast<std::size_t>(b)] -
                                            p1[static_cast<std::size_t>(a)][static_cast<std::size_t>(b)]);
                const double d2 = std::fabs(next_p2[static_cast<std::size_t>(a)][static_cast<std::size_t>(b)] -
                                            p2[static_cast<std::size_t>(a)][static_cast<std::size_t>(b)]);
                local_max = std::max(local_max, std::max(d1, d2));
            }
            row_delta[static_cast<std::size_t>(a)] = local_max;
        });

        double max_delta = 0.0;
        for (int a = 1; a <= target; ++a) {
            max_delta = std::max(max_delta, row_delta[static_cast<std::size_t>(a)]);
        }

        p1.swap(next_p1);
        p2.swap(next_p2);

        if (max_delta < tolerance) {
            break;
        }
    }

    return p1[static_cast<std::size_t>(target)][static_cast<std::size_t>(target)];
}

bool validate() {
    {
        const SolveResult tiny = solve_closed_form(1);
        const double expected = 1.0 / 3.0;
        if (std::fabs(tiny.start_probability - expected) > 1e-13) {
            std::cerr << "Validation failed for target=1: expected=" << expected
                      << ", got=" << tiny.start_probability << "\n";
            return false;
        }
    }

    const std::vector<int> deterministic_targets = {5, 10, 20};
    for (int target : deterministic_targets) {
        const double exact = solve_closed_form(target).start_probability;
        const double iter = solve_value_iteration(target, 1);
        if (std::fabs(exact - iter) > 2e-12) {
            std::cerr << "Validation failed for target=" << target
                      << ": exact=" << exact << ", value-iteration=" << iter << "\n";
            return false;
        }
    }

    unsigned hw = std::thread::hardware_concurrency();
    if (hw == 0) hw = 2;
    const unsigned threads = std::min<unsigned>(hw, 8);
    if (threads > 1) {
        const int thread_check_target = 35;
        const double single = solve_value_iteration(thread_check_target, 1);
        const double multi = solve_value_iteration(thread_check_target, threads);
        if (std::fabs(single - multi) > 1e-12) {
            std::cerr << "Thread consistency failed for target=" << thread_check_target
                      << ": single=" << single << ", multi=" << multi << "\n";
            return false;
        }
    }

    return true;
}

}  // namespace

int main(int argc, char** argv) {
    if (!validate()) {
        return 1;
    }

    int target = 100;
    if (argc > 1) {
        target = std::max(1, std::atoi(argv[1]));
    }

    const SolveResult result = solve_closed_form(target);
    std::cout << std::fixed << std::setprecision(8) << result.start_probability << '\n';
    return 0;
}

Python

def compute_t_limit(target):
    t = 1
    score = 1
    while score < target and t < 60:
        t += 1
        score <<= 1
    return t + 4

def p1_turn_value(f, a, b):
    if a <= 0: return 0.0
    if b <= 0: return 1.0
    prev_row = 0.0 if a == 1 else f[a - 1][b]
    return 0.5 * (f[a][b] + prev_row)

def solve(target=100):
    f = [[0.0] * (target + 1) for _ in range(target + 1)]
    
    for a in range(1, target + 1):
        f[a][0] = 1.0
        
    t_limit = compute_t_limit(target)
    q = [0.0] * (t_limit + 1)
    score = [0] * (t_limit + 1)
    
    for t in range(1, t_limit + 1):
        q[t] = 2.0 ** (-t)
        score[t] = 1 << (t - 1)
        
    for a in range(1, target + 1):
        for b in range(1, target + 1):
            prev_same_b = 0.0 if a == 1 else f[a - 1][b]
            best_value = 0.0
            
            for t in range(1, t_limit + 1):
                gain = score[t]
                if gain >= b:
                    success_value = 1.0
                else:
                    success_value = p1_turn_value(f, a, b - gain)
                    
                qt = q[t]
                candidate = (2.0 * qt * success_value + (1.0 - qt) * prev_same_b) / (1.0 + qt)
                
                if candidate > best_value:
                    best_value = candidate
                    
            f[a][b] = best_value
            
    ans = p1_turn_value(f, target, target)
    return f"{ans:.8f}"

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

Java

import java.util.Locale;

public class Euler232 {
    static int computeTLimit(int target) {
        int t = 1;
        long score = 1;
        while (score < target && t < 60) {
            t++;
            score <<= 1;
        }
        return t + 4;
    }

    static double p1TurnValue(double[][] f, int a, int b) {
        if (a <= 0)
            return 0.0;
        if (b <= 0)
            return 1.0;
        double prevRow = (a == 1) ? 0.0 : f[a - 1][b];
        return 0.5 * (f[a][b] + prevRow);
    }

    public static String solve() {
        int target = 100;
        double[][] f = new double[target + 1][target + 1];

        for (int a = 1; a <= target; ++a) {
            f[a][0] = 1.0;
        }

        int tLimit = computeTLimit(target);
        double[] q = new double[tLimit + 1];
        int[] score = new int[tLimit + 1];

        for (int t = 1; t <= tLimit; ++t) {
            q[t] = Math.pow(2.0, -t);
            score[t] = 1 << (t - 1);
        }

        for (int a = 1; a <= target; ++a) {
            for (int b = 1; b <= target; ++b) {
                double prevSameB = (a == 1) ? 0.0 : f[a - 1][b];
                double bestValue = 0.0;

                for (int t = 1; t <= tLimit; ++t) {
                    int gain = score[t];
                    double successValue = (gain >= b) ? 1.0 : p1TurnValue(f, a, b - gain);

                    double qt = q[t];
                    double candidate = (2.0 * qt * successValue + (1.0 - qt) * prevSameB) / (1.0 + qt);

                    if (candidate > bestValue) {
                        bestValue = candidate;
                    }
                }

                f[a][b] = bestValue;
            }
        }

        double startProbability = p1TurnValue(f, target, target);
        return String.format(Locale.US, "%.8f", startProbability);
    }

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