Problem 1012: Rock Paper Scissors

Project Euler Problem 1012 Solution

Certify the unique mixed equilibrium of every RPS(n) on a sliding odd interval of options, solving the support equations with a parity-separated recurrence for consecutive differences in 50-digit decimal arithmetic. Implementations are available in C++, Python and Java.

Detailed mathematical approach

Problem Summary

In $RPS(n)$ two players simultaneously choose a number from $1,\dots,n$. Equal choices are a draw. If the two numbers differ by an odd amount, the smaller number wins; if they differ by a nonzero even amount, the larger number wins. A win with the number $m$ is paid $2m-1$ dollars by the loser. The game has exactly one Nash equilibrium, a mixed strategy $p=(p_1,\dots,p_n)$, and $P(n)=p_n$ is the probability that it assigns to the largest number. We must compute

$$S(N)=\sum_{n=3}^{N}P(n),\qquad N=10^5,$$

rounded to ten places after the decimal point. The statement gives $P(3)=\tfrac19$, $P(4)=\tfrac15$, $P(10)\approx0.0479638009$, $S(10)\approx1.1546112276$ and $S(100)\approx4.8779925686$. For $n=3$ the numbers $1,2,3$ play the roles of Scissors, Paper and Rock with the payments $1,3,5$ of the statement's weighted Rock-Paper-Scissors, whose equilibrium is $(\tfrac13,\tfrac59,\tfrac19)$ for the numbers $1,2,3$. ^[1]

Mathematical Approach

An equilibrium is described by linear equations on the numbers it uses and by inequalities on the others. The parity rule of the game turns these equations into a short recurrence for consecutive differences of probabilities, so each new equilibrium reduces to a $3\times3$ linear system on a short interval of numbers ending at $n$. As $n$ grows, this interval slides upwards and now and then widens by two. Every candidate is checked against all equilibrium conditions before it is accepted, so the structural guesses only steer the search; they never decide the answer.

1. The payoff matrix and the equilibrium conditions

Let $a_{ij}$ be the payment received by a player who chooses $i$ when the opponent chooses $j$. It is $0$ for $i=j$. Otherwise the winning number $w$ is $\min(i,j)$ when $i-j$ is odd and $\max(i,j)$ when $i-j$ is even, and

$$a_{ij}=\begin{cases}\phantom{-}(2w-1), & w=i,\\ -(2w-1), & w=j.\end{cases}$$

Swapping the players swaps winner and loser, so $a_{ji}=-a_{ij}$: the matrix $A=(a_{ij})$ is skew-symmetric. Against a mixed strategy $p$, the pure counter-strategy $i$ earns

$$E_i=\sum_{j=1}^{n}a_{ij}\,p_j .$$

A mixed counter-strategy $q$ earns $\sum_i q_iE_i\le\max_iE_i$, so $p$ is an equilibrium in the sense of the statement exactly when $E_i\le0$ for every $i$. Moreover $\sum_i p_iE_i=p^{\mathsf T}Ap=0$ for every vector $p$, because $x^{\mathsf T}Ax=0$ for every skew-symmetric matrix. ^[2] A sum of nonpositive terms $p_iE_i$ vanishes only if every term does, so the equilibrium conditions are

$$p_i\ge0,\qquad \sum_{i}p_i=1,\qquad p_i\gt0\ \Rightarrow\ E_i=0,\qquad p_i=0\ \Rightarrow\ E_i\le0.$$

For $n=3$ the payoffs are $E_1=p_2-5p_3$, $E_2=-p_1+3p_3$ and $E_3=5p_1-3p_2$. Setting all three to zero gives $p_1=3p_3$ and $p_2=5p_3$, and the normalisation yields $p=(\tfrac13,\tfrac59,\tfrac19)$, so $P(3)=\tfrac19$.

2. Payoffs split by parity

Write $j\equiv i$ when $j$ and $i$ have the same parity. Against $j\ne i$, the option $i$ wins exactly when $j\equiv i$ and $j\lt i$ (even difference, the larger number wins), or when $j\not\equiv i$ and $j\gt i$ (odd difference, the smaller number wins). In the two remaining cases it loses and pays $2j-1$. Hence

$$E_i=(2i-1)\,X_i-\sum_{\substack{j\lt i\\ j\not\equiv i}}(2j-1)\,p_j-\sum_{\substack{j\gt i\\ j\equiv i}}(2j-1)\,p_j,\qquad X_i=\sum_{\substack{j\lt i\\ j\equiv i}}p_j+\sum_{\substack{j\gt i\\ j\not\equiv i}}p_j .$$

Here $X_i$ is the probability that the option $i$ wins. Let $Y_i$ be the probability that it loses, so that $X_i+Y_i+p_i=1$. The option $i+1$ wins against exactly the choices that beat $i$: a smaller number with the parity of $i+1$ is a number below $i$ of the opposite parity, and a larger number of the opposite parity to $i+1$ is a number above $i$ of the same parity as $i$. Therefore $X_{i+1}=Y_i$, and

$$X_i+X_{i+1}=1-p_i .$$

3. Options below and above the support

Let $S=\{i:p_i\gt0\}$ be the support, with smallest element $f$ and largest element $\ell$. For a parity class $t$ put $m_t=\sum_{j\in S,\,j\equiv t}p_j$ and $W_t=\sum_{j\in S,\,j\equiv t}(2j-1)\,p_j$, and let $r$ be the parity of $i$ and $\bar r$ the other parity. If $i\lt f$, every chosen number lies above $i$, and the formula of Section 2 becomes

$$E_i=(2i-1)\,m_{\bar r}-W_r\qquad(i\lt f).$$

Within one parity class this increases with $i$. Consequently $E_i\le0$ for all $i\lt f$ follows from $E_{f-1}\le0$ and $E_{f-2}\le0$, the largest options of each parity below the support. Symmetrically, $E_i=(2i-1)\,m_r-W_{\bar r}$ for $i\gt\ell$, which also increases with $i$, so above the support the largest options are the dangerous ones.

The payoffs between the numbers $1,\dots,n-1$ do not depend on $n$. Hence the equilibrium of $RPS(n-1)$ is still an equilibrium of $RPS(n)$ exactly when the new option earns $E_n\le0$. Because the equilibrium is unique, this means $P(n)=0$. It happens for $n=8,21,47,88,150,\dots$, for $52$ values of $n$ up to $10^5$ in total.

4. A recurrence for consecutive differences

Moving from $i$ to $i+2$ keeps the parity and changes each of the four sums of Section 2 by one term: $p_i$ joins the same-parity numbers below, $p_{i+1}$ leaves the opposite-parity numbers above, $(2i+1)\,p_{i+1}$ joins the weighted opposite-parity sum below, and $(2i+3)\,p_{i+2}$ leaves the weighted same-parity sum above. Together with $(2i+3)-(2i-1)=4$, this gives, for every vector $p$,

$$E_{i+2}-E_i=4X_i+(2i+3)\,(p_i-p_{i+1}+p_{i+2})-(2i+1)\,p_{i+1}.$$

The sum $X_i$ reaches over the whole range, but $X_i+X_{i+1}=1-p_i$ is local. Adding the identity for $i$ and for $i+1$ and writing $d_i=p_{i+1}-p_i$, all other terms collapse:

$$\bigl(E_{i+2}+E_{i+3}\bigr)-\bigl(E_i+E_{i+1}\bigr)=4-(2i-1)\,d_i+(2i+5)\,d_{i+2}.$$

If $i$, $i+1$, $i+2$ and $i+3$ all belong to the support, the left side vanishes and

$$(2i+5)\,d_{i+2}=(2i-1)\,d_i-4 .$$

The recurrence links $d_i$ only with $d_{i+2}$: the differences that start at even positions and those that start at odd positions form two separate first-order chains. This is the parity separation used by the code.

5. Odd supports and three boundary equations

On the support the conditions read $A_{SS}\,p_S=0$ with $p_S\ne0$, where $A_{SS}$ is the skew-symmetric submatrix of the chosen numbers, so $A_{SS}$ must be singular. Every skew-symmetric matrix of odd order is singular. The orthogonal normal form of a real skew-symmetric matrix consists of $2\times2$ blocks $\left(\begin{smallmatrix}0&\lambda\\-\lambda&0\end{smallmatrix}\right)$ and zeros, so its rank is even. ^[2] For a support of even size, the kernel of $A_{SS}$ would therefore have dimension at least $2$ and contain a vector $q\ne0$ with $\sum_j q_j=0$. If every option outside $S$ loses strictly, $p+\varepsilon q$ would then be a second equilibrium for small $|\varepsilon|$, contradicting uniqueness. The search therefore considers supports of odd size only. As elsewhere, the program does not rely on this argument, because it verifies every candidate.

For an interval support $\{f,\dots,n\}$ of odd size $s\ge3$, the recurrence determines everything from three numbers $p_f$, $d_f$ and $d_{f+1}$:

$$p_{f+1}=p_f+d_f,\qquad p_{f+2}=p_{f+1}+d_{f+1},\qquad d_{i+2}=\frac{(2i-1)\,d_i-4}{2i+5}\quad(f\le i\le n-3).$$

Each $p_j$ is therefore an affine function of the three unknowns, and three linear equations fix them:

$$\sum_{j=f}^{n}p_j=1,\qquad E_f=0,\qquad E_{f+1}=0.$$

The remaining support equations follow, except in a degenerate case. Put $F_j=E_j+E_{j+1}$. By the identity of Section 4 and the normalisation, the recurrence states $F_{j+2}=F_j$ for $f\le j\le n-3$, and $E_f=E_{f+1}=0$ gives $F_f=0$. Hence $F_{f+2t}=0$ and $F_{f+2t+1}=F_{f+1}=E_{f+2}=:c$ for all admissible $t$, which forces $E_{f+2t}=tc$ and $E_{f+2t+1}=-tc$. Finally $\sum_j p_jE_j=0$ gives

$$c\sum_{t\ge0}t\,\bigl(p_{f+2t}-p_{f+2t+1}\bigr)=0\qquad(p_{n+1}=0),$$

so $c=0$ unless this weighted sum vanishes. The program does not depend on this step either: it recomputes every $E_j$ on the support.

6. Worked example: $RPS(10)$

For $n=10$ the support is $\{6,\dots,10\}$, of size $5$, and the recurrence is used twice:

$$17\,d_8=11\,d_6-4,\qquad 19\,d_9=13\,d_7-4.$$

The three boundary equations give $p_6=\tfrac{93}{1105}$, $d_6=\tfrac{14}{65}$ and $d_7=\tfrac{36}{1105}$. The recurrence then gives $d_8=-\tfrac{106}{1105}$ and $d_9=-\tfrac{16}{85}$, so

$$\bigl(p_6,p_7,p_8,p_9,p_{10}\bigr)=\Bigl(\frac{93}{1105},\frac{331}{1105},\frac{367}{1105},\frac{261}{1105},\frac{53}{1105}\Bigr),$$

and $P(10)=\tfrac{53}{1105}\approx0.0479638009$, as in the statement. The two nearest options below the support lose, with $E_5=-\tfrac{4123}{1105}$ and $E_4=-\tfrac{3391}{1105}$, so by Section 3 every smaller option loses as well. The table lists the first equilibria. At $n=8$ the new option earns $E_8=-\tfrac7{11}$ against the equilibrium of $RPS(7)$, so $P(8)=0$, and at $n=9$ the support widens to five numbers.

$n$	support	$P(n)$
$3$	$\{1,2,3\}$	$1/9$
$4$	$\{2,3,4\}$	$1/5$
$5$	$\{3,4,5\}$	$5/21$
$6$	$\{4,5,6\}$	$7/27$
$7$	$\{5,6,7\}$	$3/11$
$8$	$\{5,6,7\}$	$0$
$9$	$\{5,\dots,9\}$	$7/275$
$10$	$\{6,\dots,10\}$	$53/1105$
$11$	$\{7,\dots,11\}$	$31/475$
$12$	$\{8,\dots,12\}$	$47/595$
$13$	$\{9,\dots,13\}$	$197/2185$

The first eight values add up to exactly $S(10)=\tfrac{13262413}{11486475}\approx1.1546112276$.

7. Following the support as $n$ grows

Start from the equilibrium of $RPS(n-1)$, whose support is the interval $\{f,\dots,\ell\}$. If $E_n\le0$, it remains the equilibrium and $P(n)=0$. Otherwise the new support must contain $n$, and the first candidate is $\{f',\dots,n\}$ with $f'=\min\bigl(n-2,\ f+((n-f)\bmod 2)\bigr)$. When $\ell=n-1$, this slides the interval up by one; after a step with $P(n-1)=0$ we have $\ell=n-2$, and it keeps $f$ and widens the interval by two. A candidate with a negative probability is narrowed ($f'$ increases by $2$); a candidate with a profitable option below it is widened ($f'$ decreases by $2$).

In a single pass over $n\le10^5$, the first candidate was accepted every time: the $99\,946$ positive values of $P(n)$ needed exactly $99\,946$ interval solutions. The width of the support grows roughly like $(12n)^{1/3}$: it is $23$ at $n=10^3$, $49$ at $n=10^4$ and $107$ at $n=10^5$. The steps with $P(n)=0$, namely $8,21,47,88,150,235,\dots$, have second differences $13,15,21,23,29,31,\dots$, that is $13+4m-2\,(m\bmod 2)$ for $m=0,1,2,\dots$, for all $52$ of them up to $10^5$. Third differences of average $4$ make the $m$-th such step grow like $\tfrac23m^3$; since the support has width about $2m$ there, this agrees with the $(12n)^{1/3}$ growth. These are observations, and the algorithm does not use them.

How the Code Works

payoff evaluates $a_{ij}$ directly. interval_strategy(first, last) stores the coefficients of every $p_j$ as a vector of four numbers: one for each unknown $p_f$, $d_f$, $d_{f+1}$ and one for the constant term. The two current differences sit in a two-element array indexed by the parity of the offset, so one loop advances both chains of Section 4. The function then accumulates the three boundary equations: the normalisation, and $E_f$ and $E_{f+1}$ divided by $2f-1$ and $2f+1$ to keep the coefficients moderate. It solves the $3\times3$ system by Gauss–Jordan elimination with partial pivoting and substitutes the solution.

check_strategy certifies a candidate. It checks nonnegative probabilities, total mass $1$, and $E_j=0$ at every support position, computed in one pass with running sums per parity class, exactly as in Section 2. It also checks the two nearest options below (Section 3) and every option above the support up to $n$. advance implements Section 7. initial_strategy finds an equilibrium from scratch by trying the odd widths $3,5,7,\dots$ for intervals that end at $n$ or $n-1$.

solve splits the range $3\le n\le N$ into contiguous blocks, one per hardware thread, and runs them as POSIX threads. Each block starts from initial_strategy(first - 1), so the blocks are independent, and their sums are added in block order. The arithmetic type is Boost.Multiprecision's cpp_dec_float_50, a radix-10 type with $50$ decimal digits and extra internal guard digits. ^[3] Every sign test uses the tolerance $\varepsilon=10^{-30}$. The program prints the sum with ten decimals; with --self-test it runs the checks listed below instead.

Precision and the three implementations

The probabilities are rational numbers, but their denominators grow quickly: the exact denominators already have $13$ digits at $n=100$ and $23$ digits at $n=400$. All three versions therefore use $50$-digit decimal arithmetic instead of exact fractions. Measured over the whole run, the support equations hold to within about $10^{-55}$ in sampled checks; the nearest option below a support loses at least $1.4\cdot10^{-4}$; the new option at a step with $P(n)=0$ loses at least $9.4\cdot10^{-3}$; and the smallest probability on a support is about $4.6\cdot10^{-8}$. The tolerance $10^{-30}$ is thus far above the rounding noise and far below every genuine quantity, and the final sum is not close to a rounding boundary at the tenth decimal.

Python uses the decimal module with prec = 50, whose default precision is $28$ digits. ^[4] It evaluates the same blocks one after another. Java uses BigDecimal with MathContext(50, HALF_EVEN), which is required because a division without a context throws an exception when the quotient does not terminate. ^[5] An ExecutorService runs the Java blocks in parallel. All three versions print the same digits. On a 16-thread machine the C++ version takes about two seconds, the Java version about four seconds and the single-threaded Python version about forty seconds.

Complexity and Verification

In practice every $n$ needs one interval solution. Its recurrence, its $3\times3$ coefficients and its certificate each take $O(s)$ operations for a support of width $s$, so the total work is $O\bigl(\sum_{n\le N}s(n)\bigr)$. With $s(n)\approx(12n)^{1/3}$ this is $O(N^{4/3})$. For $N=10^5$ the support has $79.7$ numbers on average, about $8\cdot10^6$ support entries in total. The memory is $O(s)$ per thread. Restarting a block costs one search over about $s/2$ widths, which is negligible.

The program runs these checks with --self-test: ^[7]

For every $2\le n\le200$, the support probabilities agree to within $10^{-30}$ with an independent dense Gauss–Jordan solution of the payoff equations on the support, with the last equation replaced by the normalisation.
For the same $n$, every pure counter-strategy $1\le i\le n$ earns $E_i\le\varepsilon$.
$P(3)=\tfrac19$, $P(4)=\tfrac15$ and $P(8)=0$; $P(10)$, $S(10)$ and $S(100)$ agree with the rounded values of the statement to within $5\cdot10^{-11}$.
$S(1000)$ computed with one thread and with all threads agree to within $10^{-30}$. Python compares one block with as many blocks as there are processors, evaluated in sequence.

Separate computations checked the method independently. For $3\le n\le13$, an exhaustive search over all supports in exact rational arithmetic found exactly one equilibrium each time; its support was always an interval of odd size, and its values equal those of the fast method. For $n\le400$, an exact rational re-run of the sliding search, which solves each support by dense elimination without the recurrence, certified every equilibrium exactly, with strict inequalities outside the support, and matched the $50$-digit values to within $10^{-45}$. For every $3\le n\le300$, a linear-programming solution of the full $n\times n$ game with the HiGHS solver in SciPy, which makes no assumption about the support, agreed with $P(n)$ to within $10^{-12}$. ^[6] At $N=10^5$ the C++, Python and Java versions print identical results. These comparisons are additional validation, separate from the packaged checks.

Footnotes and References

The references below supply the background facts used above. The parity decomposition, the difference recurrence, the boundary system and the observations about the support are derived explicitly in this article.

Project Euler 1012 — Rock Paper Scissors. The official statement defines RPS(n), states that each RPS(n) has exactly one Nash equilibrium, and gives the checkpoints P(3), P(4), P(10), S(10) and S(100). These examples test the implementation without disclosing the requested final value.
Wikipedia — Skew-symmetric matrix. The article states that xᵀAx = 0 for every skew-symmetric matrix A, that every skew-symmetric matrix of odd dimension is singular, and that a real skew-symmetric matrix has an orthogonal normal form made of 2×2 blocks with entries ±λ and zeros. The even rank used in Section 5 follows from that normal form.
Boost.Multiprecision — cpp_dec_float. The typedef cpp_dec_float_50 provides 50 decimal digits of precision in radix 10, with internal guard digits beyond the stated precision.
Python — decimal module. The context precision prec, 28 digits by default, sets the precision of arithmetic results, which are rounded with ROUND_HALF_EVEN in the default context.
Java — BigDecimal. With a MathContext, results are rounded to the given precision; without one, a division whose exact quotient has a non-terminating decimal expansion throws an ArithmeticException. setScale with a RoundingMode formats the final result.
SciPy — scipy.optimize.linprog. Its default method is HiGHS, which was used for the independent linear-programming comparison; that comparison is not part of the published solutions.
C++ — Euler1012.cpp. The linked immutable C++ revision contains the interval solver, the support search, the worker threads and the checks run by --self-test. It is the implementation record for this article. The separate exact and linear-programming comparisons described above are additional validation, not tests packaged in that source file.

Problem 1012 source code

C++

#include <algorithm>
#include <array>
#include <boost/multiprecision/cpp_dec_float.hpp>
#include <cstdlib>
#include <exception>
#include <iomanip>
#include <iostream>
#include <pthread.h>
#include <stdexcept>
#include <string>
#include <thread>
#include <vector>

namespace {

using Real = boost::multiprecision::cpp_dec_float_50;
using Basis = std::array<Real, 4>;
constexpr int TARGET = 100'000;
const Real EPS("1e-30");

void require(const bool condition, const std::string& description) {
    if (!condition) throw std::runtime_error("Check failed: " + description);
}

int payoff(const int i, const int j) {
    if (i == j) return 0;
    const int winner = (i - j) % 2 != 0 ? std::min(i, j) : std::max(i, j);
    return (winner == i ? 1 : -1) * (2 * winner - 1);
}

struct Strategy {
    int first;
    std::vector<Real> probability;

    int last() const { return first + static_cast<int>(probability.size()) - 1; }
};

Real expected_payoff(const Strategy& strategy, const int i) {
    Real result = 0;
    for (int j = strategy.first; j <= strategy.last(); ++j) {
        result += payoff(i, j) * strategy.probability[j - strategy.first];
    }
    return result;
}

Strategy interval_strategy(const int first, const int last) {
    const int count = last - first + 1;
    require(count >= 3 && count % 2 != 0, "odd support size");
    std::vector<Basis> basis(count);
    basis[0][0] = 1;
    basis[1][0] = basis[1][1] = 1;
    basis[2][0] = basis[2][1] = basis[2][2] = 1;
    std::array<Basis, 2> difference{};
    difference[0][1] = difference[1][2] = 1;

    // For d_i = p_(i+1)-p_i, (2i+5)d_(i+2) = (2i-1)d_i - 4.
    for (int offset = 0; offset < count - 3; ++offset) {
        const int i = first + offset;
        Basis& d = difference[offset % 2];
        for (Real& value : d) value = value * (2 * i - 1) / (2 * i + 5);
        d[3] -= Real(4) / (2 * i + 5);
        for (int q = 0; q < 4; ++q) basis[offset + 3][q] = basis[offset + 2][q] + d[q];
    }

    std::array<Basis, 3> matrix{};
    for (int offset = 0; offset < count; ++offset) {
        for (int q = 0; q < 4; ++q) {
            matrix[0][q] += basis[offset][q];
            matrix[1][q] += Real(payoff(first, first + offset)) / (2 * first - 1) * basis[offset][q];
            matrix[2][q] += Real(payoff(first + 1, first + offset)) / (2 * first + 1) * basis[offset][q];
        }
    }
    matrix[0][3] = 1 - matrix[0][3];
    matrix[1][3] = -matrix[1][3];
    matrix[2][3] = -matrix[2][3];

    for (int col = 0; col < 3; ++col) {
        int pivot = col;
        for (int row = col + 1; row < 3; ++row) {
            if (abs(matrix[row][col]) > abs(matrix[pivot][col])) pivot = row;
        }
        std::swap(matrix[col], matrix[pivot]);
        require(matrix[col][col] != 0, "nonsingular boundary equations");
        const Real divisor = matrix[col][col];
        for (int q = col; q < 4; ++q) matrix[col][q] /= divisor;
        for (int row = 0; row < 3; ++row) {
            if (row == col) continue;
            const Real multiplier = matrix[row][col];
            for (int q = col; q < 4; ++q) matrix[row][q] -= multiplier * matrix[col][q];
        }
    }

    Strategy result{first, std::vector<Real>(count)};
    for (int offset = 0; offset < count; ++offset) {
        result.probability[offset] = basis[offset][3];
        for (int q = 0; q < 3; ++q) result.probability[offset] += basis[offset][q] * matrix[q][3];
    }
    return result;
}

bool lower_options_unprofitable(const Strategy& strategy) {
    for (int i = std::max(1, strategy.first - 2); i < strategy.first; ++i) {
        if (expected_payoff(strategy, i) > EPS) return false;
    }
    return true;
}

bool probabilities_nonnegative(const Strategy& strategy) {
    return std::all_of(strategy.probability.begin(), strategy.probability.end(),
        [](const Real& value) { return value >= -EPS; });
}

void check_strategy(const Strategy& strategy, const int n) {
    std::array<Real, 2> mass{}, weight{}, prefix_mass{}, prefix_weight{};
    for (int i = strategy.first; i <= strategy.last(); ++i) {
        const Real& p = strategy.probability[i - strategy.first];
        require(p >= -EPS, "nonnegative equilibrium probability");
        mass[i % 2] += p;
        weight[i % 2] += (2 * i - 1) * p;
    }
    require(abs(mass[0] + mass[1] - 1) < EPS, "normalized equilibrium");
    for (int i = strategy.first; i <= strategy.last(); ++i) {
        const int parity = i % 2;
        const Real& p = strategy.probability[i - strategy.first];
        const Real value = (2 * i - 1) * (prefix_mass[parity] + mass[1 - parity] - prefix_mass[1 - parity] + p)
            + prefix_weight[parity] - prefix_weight[1 - parity] - weight[parity];
        require(abs(value) < EPS, "zero payoff on equilibrium support");
        prefix_mass[parity] += p;
        prefix_weight[parity] += (2 * i - 1) * p;
    }
    require(lower_options_unprofitable(strategy), "unprofitable lower options");
    for (int i = strategy.last() + 1; i <= n; ++i) {
        require(expected_payoff(strategy, i) <= EPS, "unprofitable upper options");
    }
}

Strategy initial_strategy(const int n) {
    if (n <= 2) return {1, {Real(1)}};
    for (int count = 3; count <= n; count += 2) {
        for (const int last : {n, n - 1}) {
            if (last < count) continue;
            Strategy candidate = interval_strategy(last - count + 1, last);
            if (!probabilities_nonnegative(candidate) || !lower_options_unprofitable(candidate)) continue;
            if (last < n && expected_payoff(candidate, n) > EPS) continue;
            check_strategy(candidate, n);
            return candidate;
        }
    }
    throw std::runtime_error("No equilibrium support found");
}

Real advance(Strategy& strategy, const int n) {
    if (expected_payoff(strategy, n) <= EPS) return 0;
    int first = std::min(n - 2, strategy.first + (n - strategy.first) % 2);
    for (int attempt = 0; attempt <= n; ++attempt) {
        require(first >= 1 && first <= n - 2, "valid support endpoints");
        Strategy candidate = interval_strategy(first, n);
        if (!probabilities_nonnegative(candidate)) {
            first += 2;
        } else if (!lower_options_unprofitable(candidate)) {
            first -= 2;
        } else {
            check_strategy(candidate, n);
            strategy = std::move(candidate);
            return strategy.probability.back();
        }
    }
    throw std::runtime_error("Equilibrium support search did not converge");
}

struct Task {
    int first = 0;
    int last = 0;
    Real sum = 0;
    std::exception_ptr error;
};

void* sum_worker(void* argument) {
    Task& task = *static_cast<Task*>(argument);
    try {
        Strategy strategy = initial_strategy(task.first - 1);
        for (int n = task.first; n <= task.last; ++n) task.sum += advance(strategy, n);
    } catch (...) {
        task.error = std::current_exception();
    }
    return nullptr;
}

Real solve(const int n, unsigned thread_count) {
    if (n < 3) return 0;
    thread_count = std::min(thread_count, static_cast<unsigned>(n - 2));
    require(thread_count > 0, "positive thread count");
    std::vector<Task> tasks(thread_count);
    std::vector<pthread_t> threads(thread_count);
    unsigned created = 0;
    for (unsigned t = 0; t < thread_count; ++t) {
        tasks[t].first = 3 + static_cast<int>(static_cast<long long>(n - 2) * t / thread_count);
        tasks[t].last = 2 + static_cast<int>(static_cast<long long>(n - 2) * (t + 1) / thread_count);
        if (thread_count == 1) {
            sum_worker(&tasks[t]);
        } else {
            if (pthread_create(&threads[t], nullptr, sum_worker, &tasks[t]) != 0) break;
            ++created;
        }
    }
    bool joined = true;
    for (unsigned t = 0; t < created; ++t) joined = pthread_join(threads[t], nullptr) == 0 && joined;
    require(thread_count == 1 || created == thread_count, "pthread_create");
    require(joined, "pthread_join");
    Real sum = 0;
    for (const Task& task : tasks) {
        if (task.error) std::rethrow_exception(task.error);
        sum += task.sum;
    }
    return sum;
}

std::vector<Real> dense_equilibrium(const Strategy& strategy) {
    const int count = static_cast<int>(strategy.probability.size());
    std::vector<std::vector<Real>> matrix(count, std::vector<Real>(count + 1));
    for (int row = 0; row < count - 1; ++row) {
        for (int col = 0; col < count; ++col) matrix[row][col] = payoff(strategy.first + row, strategy.first + col);
    }
    std::fill(matrix.back().begin(), matrix.back().end(), Real(1));
    for (int col = 0; col < count; ++col) {
        int pivot = col;
        for (int row = col + 1; row < count; ++row) {
            if (abs(matrix[row][col]) > abs(matrix[pivot][col])) pivot = row;
        }
        std::swap(matrix[col], matrix[pivot]);
        require(matrix[col][col] != 0, "dense equilibrium pivot");
        const Real divisor = matrix[col][col];
        for (int q = col; q <= count; ++q) matrix[col][q] /= divisor;
        for (int row = 0; row < count; ++row) {
            if (row == col) continue;
            const Real multiplier = matrix[row][col];
            for (int q = col; q <= count; ++q) matrix[row][q] -= multiplier * matrix[col][q];
        }
    }
    std::vector<Real> result(count);
    for (int row = 0; row < count; ++row) result[row] = matrix[row][count];
    return result;
}

void run_tests(const unsigned thread_count) {
    Strategy strategy{1, {Real(1)}};
    Real sum = 0;
    const Real rounded_tolerance("5e-11");
    for (int n = 2; n <= 200; ++n) {
        const Real p = advance(strategy, n);
        if (n >= 3) sum += p;
        const std::vector<Real> independent = dense_equilibrium(strategy);
        for (std::size_t i = 0; i < independent.size(); ++i) {
            require(abs(independent[i] - strategy.probability[i]) < EPS, "dense payoff-matrix comparison");
        }
        for (int i = 1; i <= n; ++i) require(expected_payoff(strategy, i) <= EPS, "all pure counter-strategies");
        if (n == 3) require(abs(p - Real(1) / 9) < EPS, "P(3) = 1/9");
        if (n == 4) require(abs(p - Real(1) / 5) < EPS, "P(4) = 1/5");
        if (n == 8) require(p == 0, "unused final option for n=8");
        if (n == 10) {
            require(abs(p - Real("0.0479638009")) < rounded_tolerance, "P(10)");
            require(abs(sum - Real("1.1546112276")) < rounded_tolerance, "S(10)");
        }
        if (n == 100) require(abs(sum - Real("4.8779925686")) < rounded_tolerance, "S(100)");
    }
    require(abs(solve(1000, 1) - solve(1000, thread_count)) < EPS, "thread consistency");
    std::cout << "All checks passed.\n";
}

}  // namespace

int main(int argc, char* argv[]) {
    try {
        const unsigned thread_count = std::max(1U, std::thread::hardware_concurrency());
        if (argc == 2 && std::string(argv[1]) == "--self-test") {
            run_tests(thread_count);
            return EXIT_SUCCESS;
        }
        if (argc != 1) throw std::invalid_argument("Usage: Euler1012 [--self-test]");
        std::cout << std::fixed << std::setprecision(10) << solve(TARGET, thread_count) << '\n';
    } catch (const std::exception& error) {
        std::cerr << error.what() << '\n';
        return EXIT_FAILURE;
    }
    return EXIT_SUCCESS;
}

Python

#!/usr/bin/env python3
"""Project Euler Problem 1012 - Rock Paper Scissors."""

import os
import sys
from decimal import Decimal, getcontext

getcontext().prec = 50

TARGET = 100_000
EPS = Decimal("1e-30")
ZERO = Decimal(0)
ONE = Decimal(1)


def require(condition, description):
    if not condition:
        raise AssertionError("Check failed: " + description)


def payoff(i, j):
    if i == j:
        return 0
    winner = min(i, j) if (i - j) % 2 != 0 else max(i, j)
    return (1 if winner == i else -1) * (2 * winner - 1)


class Strategy:
    def __init__(self, first, probability):
        self.first = first
        self.probability = probability

    def last(self):
        return self.first + len(self.probability) - 1


def expected_payoff(strategy, i):
    result = ZERO
    j = strategy.first
    for p in strategy.probability:
        result += payoff(i, j) * p
        j += 1
    return result


def interval_strategy(first, last):
    count = last - first + 1
    require(count >= 3 and count % 2 != 0, "odd support size")
    # Unknowns x0, x1, x2 with p_first = x0, d_first = x1, d_(first+1) = x2; column 3 is the constant term.
    basis = [[ONE, ZERO, ZERO, ZERO], [ONE, ONE, ZERO, ZERO], [ONE, ONE, ONE, ZERO]]
    difference = [[ZERO, ONE, ZERO, ZERO], [ZERO, ZERO, ONE, ZERO]]

    # For d_i = p_(i+1)-p_i, (2i+5)d_(i+2) = (2i-1)d_i - 4.
    for offset in range(count - 3):
        i = first + offset
        d = difference[offset % 2]
        for q in range(4):
            d[q] = d[q] * (2 * i - 1) / (2 * i + 5)
        d[3] -= Decimal(4) / (2 * i + 5)
        previous = basis[offset + 2]
        basis.append([previous[q] + d[q] for q in range(4)])

    matrix = [[ZERO] * 4 for _ in range(3)]
    low, high = 2 * first - 1, 2 * first + 1
    for offset in range(count):
        row = basis[offset]
        a = Decimal(payoff(first, first + offset)) / low
        b = Decimal(payoff(first + 1, first + offset)) / high
        for q in range(4):
            matrix[0][q] += row[q]
            matrix[1][q] += a * row[q]
            matrix[2][q] += b * row[q]
    matrix[0][3] = 1 - matrix[0][3]
    matrix[1][3] = -matrix[1][3]
    matrix[2][3] = -matrix[2][3]

    for col in range(3):
        pivot = max(range(col, 3), key=lambda row: abs(matrix[row][col]))
        matrix[col], matrix[pivot] = matrix[pivot], matrix[col]
        require(matrix[col][col] != 0, "nonsingular boundary equations")
        divisor = matrix[col][col]
        for q in range(col, 4):
            matrix[col][q] /= divisor
        for row in range(3):
            if row == col:
                continue
            multiplier = matrix[row][col]
            for q in range(col, 4):
                matrix[row][q] -= multiplier * matrix[col][q]

    x0, x1, x2 = matrix[0][3], matrix[1][3], matrix[2][3]
    probability = [row[3] + row[0] * x0 + row[1] * x1 + row[2] * x2 for row in basis]
    return Strategy(first, probability)


def lower_options_unprofitable(strategy):
    for i in range(max(1, strategy.first - 2), strategy.first):
        if expected_payoff(strategy, i) > EPS:
            return False
    return True


def probabilities_nonnegative(strategy):
    return all(value >= -EPS for value in strategy.probability)


def check_strategy(strategy, n):
    mass, weight = [ZERO, ZERO], [ZERO, ZERO]
    prefix_mass, prefix_weight = [ZERO, ZERO], [ZERO, ZERO]
    for i, p in enumerate(strategy.probability, strategy.first):
        require(p >= -EPS, "nonnegative equilibrium probability")
        mass[i % 2] += p
        weight[i % 2] += (2 * i - 1) * p
    require(abs(mass[0] + mass[1] - 1) < EPS, "normalized equilibrium")
    for i, p in enumerate(strategy.probability, strategy.first):
        parity = i % 2
        value = ((2 * i - 1) * (prefix_mass[parity] + mass[1 - parity] - prefix_mass[1 - parity] + p)
                 + prefix_weight[parity] - prefix_weight[1 - parity] - weight[parity])
        require(abs(value) < EPS, "zero payoff on equilibrium support")
        prefix_mass[parity] += p
        prefix_weight[parity] += (2 * i - 1) * p
    require(lower_options_unprofitable(strategy), "unprofitable lower options")
    for i in range(strategy.last() + 1, n + 1):
        require(expected_payoff(strategy, i) <= EPS, "unprofitable upper options")


def initial_strategy(n):
    if n <= 2:
        return Strategy(1, [ONE])
    for count in range(3, n + 1, 2):
        for last in (n, n - 1):
            if last < count:
                continue
            candidate = interval_strategy(last - count + 1, last)
            if not probabilities_nonnegative(candidate) or not lower_options_unprofitable(candidate):
                continue
            if last < n and expected_payoff(candidate, n) > EPS:
                continue
            check_strategy(candidate, n)
            return candidate
    raise AssertionError("No equilibrium support found")


def advance(strategy, n):
    if expected_payoff(strategy, n) <= EPS:
        return ZERO
    first = min(n - 2, strategy.first + (n - strategy.first) % 2)
    for _ in range(n + 1):
        require(1 <= first <= n - 2, "valid support endpoints")
        candidate = interval_strategy(first, n)
        if not probabilities_nonnegative(candidate):
            first += 2
        elif not lower_options_unprofitable(candidate):
            first -= 2
        else:
            check_strategy(candidate, n)
            strategy.first, strategy.probability = candidate.first, candidate.probability
            return strategy.probability[-1]
    raise AssertionError("Equilibrium support search did not converge")


def solve(n, blocks=1):
    # The C++ version gives each thread one block of n and restarts it with initial_strategy;
    # this port evaluates the same blocks one after another.
    if n < 3:
        return ZERO
    blocks = min(blocks, n - 2)
    total = ZERO
    for t in range(blocks):
        first = 3 + (n - 2) * t // blocks
        last = 2 + (n - 2) * (t + 1) // blocks
        strategy = initial_strategy(first - 1)
        part = ZERO
        for m in range(first, last + 1):
            part += advance(strategy, m)
        total += part
    return total


def dense_equilibrium(strategy):
    count = len(strategy.probability)
    matrix = [[Decimal(payoff(strategy.first + row, strategy.first + col)) for col in range(count)] + [ZERO]
              for row in range(count - 1)]
    matrix.append([ONE] * (count + 1))
    for col in range(count):
        pivot = max(range(col, count), key=lambda row: abs(matrix[row][col]))
        matrix[col], matrix[pivot] = matrix[pivot], matrix[col]
        require(matrix[col][col] != 0, "dense equilibrium pivot")
        divisor = matrix[col][col]
        for q in range(col, count + 1):
            matrix[col][q] /= divisor
        for row in range(count):
            if row == col:
                continue
            multiplier = matrix[row][col]
            for q in range(col, count + 1):
                matrix[row][q] -= multiplier * matrix[col][q]
    return [matrix[row][count] for row in range(count)]


def run_tests(blocks):
    strategy = Strategy(1, [ONE])
    total = ZERO
    rounded_tolerance = Decimal("5e-11")
    for n in range(2, 201):
        p = advance(strategy, n)
        if n >= 3:
            total += p
        independent = dense_equilibrium(strategy)
        for value, probability in zip(independent, strategy.probability):
            require(abs(value - probability) < EPS, "dense payoff-matrix comparison")
        for i in range(1, n + 1):
            require(expected_payoff(strategy, i) <= EPS, "all pure counter-strategies")
        if n == 3:
            require(abs(p - ONE / 9) < EPS, "P(3) = 1/9")
        if n == 4:
            require(abs(p - ONE / 5) < EPS, "P(4) = 1/5")
        if n == 8:
            require(p == 0, "unused final option for n=8")
        if n == 10:
            require(abs(p - Decimal("0.0479638009")) < rounded_tolerance, "P(10)")
            require(abs(total - Decimal("1.1546112276")) < rounded_tolerance, "S(10)")
        if n == 100:
            require(abs(total - Decimal("4.8779925686")) < rounded_tolerance, "S(100)")
    require(abs(solve(1000, 1) - solve(1000, blocks)) < EPS, "block consistency")
    print("All checks passed.")


def main():
    blocks = max(1, os.cpu_count() or 1)
    if len(sys.argv) == 2 and sys.argv[1] == "--self-test":
        run_tests(blocks)
        return
    if len(sys.argv) != 1:
        raise AssertionError("Usage: Euler1012.py [--self-test]")
    print(f"{solve(TARGET):.10f}")


if __name__ == "__main__":
    try:
        main()
    except AssertionError as error:
        print(error, file=sys.stderr)
        sys.exit(1)

Java

import java.math.BigDecimal;
import java.math.MathContext;
import java.math.RoundingMode;
import java.util.ArrayList;
import java.util.List;
import java.util.concurrent.ExecutionException;
import java.util.concurrent.ExecutorService;
import java.util.concurrent.Executors;
import java.util.concurrent.Future;

public class Euler1012 {
    private static final int TARGET = 100_000;
    // Fifty significant decimal digits, like boost::multiprecision::cpp_dec_float_50 in the C++ version.
    private static final MathContext MC = new MathContext(50, RoundingMode.HALF_EVEN);
    private static final BigDecimal EPS = new BigDecimal("1e-30");
    private static final BigDecimal MINUS_EPS = EPS.negate();
    private static final BigDecimal FOUR = BigDecimal.valueOf(4);

    private static final class Strategy {
        int first;
        BigDecimal[] probability;

        Strategy(int first, BigDecimal[] probability) {
            this.first = first;
            this.probability = probability;
        }

        int last() {
            return first + probability.length - 1;
        }
    }

    private static void require(boolean condition, String description) {
        if (!condition) throw new IllegalStateException("Check failed: " + description);
    }

    private static int payoff(int i, int j) {
        if (i == j) return 0;
        int winner = (i - j) % 2 != 0 ? Math.min(i, j) : Math.max(i, j);
        return (winner == i ? 1 : -1) * (2 * winner - 1);
    }

    private static BigDecimal expectedPayoff(Strategy strategy, int i) {
        BigDecimal result = BigDecimal.ZERO;
        for (int j = strategy.first; j <= strategy.last(); ++j) {
            BigDecimal term = BigDecimal.valueOf(payoff(i, j)).multiply(strategy.probability[j - strategy.first], MC);
            result = result.add(term, MC);
        }
        return result;
    }

    private static BigDecimal[] vector(int a, int b, int c, int d) {
        return new BigDecimal[] {BigDecimal.valueOf(a), BigDecimal.valueOf(b), BigDecimal.valueOf(c), BigDecimal.valueOf(d)};
    }

    private static Strategy intervalStrategy(int first, int last) {
        int count = last - first + 1;
        require(count >= 3 && count % 2 != 0, "odd support size");
        BigDecimal[][] basis = new BigDecimal[count][];
        basis[0] = vector(1, 0, 0, 0);
        basis[1] = vector(1, 1, 0, 0);
        basis[2] = vector(1, 1, 1, 0);
        BigDecimal[][] difference = {vector(0, 1, 0, 0), vector(0, 0, 1, 0)};

        // For d_i = p_(i+1)-p_i, (2i+5)d_(i+2) = (2i-1)d_i - 4.
        for (int offset = 0; offset < count - 3; ++offset) {
            int i = first + offset;
            BigDecimal[] d = difference[offset % 2];
            BigDecimal factor = BigDecimal.valueOf(2L * i - 1);
            BigDecimal divisor = BigDecimal.valueOf(2L * i + 5);
            for (int q = 0; q < 4; ++q) d[q] = d[q].multiply(factor, MC).divide(divisor, MC);
            d[3] = d[3].subtract(FOUR.divide(divisor, MC), MC);
            BigDecimal[] next = new BigDecimal[4];
            for (int q = 0; q < 4; ++q) next[q] = basis[offset + 2][q].add(d[q], MC);
            basis[offset + 3] = next;
        }

        BigDecimal[][] matrix = {vector(0, 0, 0, 0), vector(0, 0, 0, 0), vector(0, 0, 0, 0)};
        BigDecimal low = BigDecimal.valueOf(2L * first - 1);
        BigDecimal high = BigDecimal.valueOf(2L * first + 1);
        for (int offset = 0; offset < count; ++offset) {
            BigDecimal a = BigDecimal.valueOf(payoff(first, first + offset)).divide(low, MC);
            BigDecimal b = BigDecimal.valueOf(payoff(first + 1, first + offset)).divide(high, MC);
            for (int q = 0; q < 4; ++q) {
                matrix[0][q] = matrix[0][q].add(basis[offset][q], MC);
                matrix[1][q] = matrix[1][q].add(a.multiply(basis[offset][q], MC), MC);
                matrix[2][q] = matrix[2][q].add(b.multiply(basis[offset][q], MC), MC);
            }
        }
        matrix[0][3] = BigDecimal.ONE.subtract(matrix[0][3], MC);
        matrix[1][3] = matrix[1][3].negate();
        matrix[2][3] = matrix[2][3].negate();

        for (int col = 0; col < 3; ++col) {
            int pivot = col;
            for (int row = col + 1; row < 3; ++row) {
                if (matrix[row][col].abs().compareTo(matrix[pivot][col].abs()) > 0) pivot = row;
            }
            BigDecimal[] swap = matrix[col];
            matrix[col] = matrix[pivot];
            matrix[pivot] = swap;
            require(matrix[col][col].signum() != 0, "nonsingular boundary equations");
            BigDecimal divisor = matrix[col][col];
            for (int q = col; q < 4; ++q) matrix[col][q] = matrix[col][q].divide(divisor, MC);
            for (int row = 0; row < 3; ++row) {
                if (row == col) continue;
                BigDecimal multiplier = matrix[row][col];
                for (int q = col; q < 4; ++q) {
                    matrix[row][q] = matrix[row][q].subtract(multiplier.multiply(matrix[col][q], MC), MC);
                }
            }
        }

        BigDecimal[] probability = new BigDecimal[count];
        for (int offset = 0; offset < count; ++offset) {
            BigDecimal value = basis[offset][3];
            for (int q = 0; q < 3; ++q) value = value.add(basis[offset][q].multiply(matrix[q][3], MC), MC);
            probability[offset] = value;
        }
        return new Strategy(first, probability);
    }

    private static boolean lowerOptionsUnprofitable(Strategy strategy) {
        for (int i = Math.max(1, strategy.first - 2); i < strategy.first; ++i) {
            if (expectedPayoff(strategy, i).compareTo(EPS) > 0) return false;
        }
        return true;
    }

    private static boolean probabilitiesNonnegative(Strategy strategy) {
        for (BigDecimal value : strategy.probability) {
            if (value.compareTo(MINUS_EPS) < 0) return false;
        }
        return true;
    }

    private static void checkStrategy(Strategy strategy, int n) {
        BigDecimal[] mass = {BigDecimal.ZERO, BigDecimal.ZERO};
        BigDecimal[] weight = {BigDecimal.ZERO, BigDecimal.ZERO};
        BigDecimal[] prefixMass = {BigDecimal.ZERO, BigDecimal.ZERO};
        BigDecimal[] prefixWeight = {BigDecimal.ZERO, BigDecimal.ZERO};
        for (int i = strategy.first; i <= strategy.last(); ++i) {
            BigDecimal p = strategy.probability[i - strategy.first];
            require(p.compareTo(MINUS_EPS) >= 0, "nonnegative equilibrium probability");
            mass[i % 2] = mass[i % 2].add(p, MC);
            weight[i % 2] = weight[i % 2].add(BigDecimal.valueOf(2L * i - 1).multiply(p, MC), MC);
        }
        require(mass[0].add(mass[1], MC).subtract(BigDecimal.ONE, MC).abs().compareTo(EPS) < 0, "normalized equilibrium");
        for (int i = strategy.first; i <= strategy.last(); ++i) {
            int parity = i % 2;
            BigDecimal p = strategy.probability[i - strategy.first];
            BigDecimal w = BigDecimal.valueOf(2L * i - 1);
            BigDecimal mixed = prefixMass[parity].add(mass[1 - parity], MC).subtract(prefixMass[1 - parity], MC).add(p, MC);
            BigDecimal value = w.multiply(mixed, MC).add(prefixWeight[parity], MC)
                .subtract(prefixWeight[1 - parity], MC).subtract(weight[parity], MC);
            require(value.abs().compareTo(EPS) < 0, "zero payoff on equilibrium support");
            prefixMass[parity] = prefixMass[parity].add(p, MC);
            prefixWeight[parity] = prefixWeight[parity].add(w.multiply(p, MC), MC);
        }
        require(lowerOptionsUnprofitable(strategy), "unprofitable lower options");
        for (int i = strategy.last() + 1; i <= n; ++i) {
            require(expectedPayoff(strategy, i).compareTo(EPS) <= 0, "unprofitable upper options");
        }
    }

    private static Strategy initialStrategy(int n) {
        if (n <= 2) return new Strategy(1, new BigDecimal[] {BigDecimal.ONE});
        for (int count = 3; count <= n; count += 2) {
            for (int last : new int[] {n, n - 1}) {
                if (last < count) continue;
                Strategy candidate = intervalStrategy(last - count + 1, last);
                if (!probabilitiesNonnegative(candidate) || !lowerOptionsUnprofitable(candidate)) continue;
                if (last < n && expectedPayoff(candidate, n).compareTo(EPS) > 0) continue;
                checkStrategy(candidate, n);
                return candidate;
            }
        }
        throw new IllegalStateException("No equilibrium support found");
    }

    private static BigDecimal advance(Strategy strategy, int n) {
        if (expectedPayoff(strategy, n).compareTo(EPS) <= 0) return BigDecimal.ZERO;
        int first = Math.min(n - 2, strategy.first + (n - strategy.first) % 2);
        for (int attempt = 0; attempt <= n; ++attempt) {
            require(first >= 1 && first <= n - 2, "valid support endpoints");
            Strategy candidate = intervalStrategy(first, n);
            if (!probabilitiesNonnegative(candidate)) {
                first += 2;
            } else if (!lowerOptionsUnprofitable(candidate)) {
                first -= 2;
            } else {
                checkStrategy(candidate, n);
                strategy.first = candidate.first;
                strategy.probability = candidate.probability;
                return strategy.probability[strategy.probability.length - 1];
            }
        }
        throw new IllegalStateException("Equilibrium support search did not converge");
    }

    private static BigDecimal sumBlock(int first, int last) {
        Strategy strategy = initialStrategy(first - 1);
        BigDecimal sum = BigDecimal.ZERO;
        for (int n = first; n <= last; ++n) sum = sum.add(advance(strategy, n), MC);
        return sum;
    }

    private static BigDecimal solve(int n, int threadCount) throws InterruptedException, ExecutionException {
        if (n < 3) return BigDecimal.ZERO;
        int threads = Math.min(threadCount, n - 2);
        require(threads > 0, "positive thread count");
        ExecutorService pool = Executors.newFixedThreadPool(threads);
        try {
            List<Future<BigDecimal>> parts = new ArrayList<>();
            for (int t = 0; t < threads; ++t) {
                int first = 3 + (int) ((long) (n - 2) * t / threads);
                int last = 2 + (int) ((long) (n - 2) * (t + 1) / threads);
                parts.add(pool.submit(() -> sumBlock(first, last)));
            }
            BigDecimal sum = BigDecimal.ZERO;
            for (Future<BigDecimal> part : parts) sum = sum.add(part.get(), MC);
            return sum;
        } finally {
            pool.shutdown();
        }
    }

    private static BigDecimal[] denseEquilibrium(Strategy strategy) {
        int count = strategy.probability.length;
        BigDecimal[][] matrix = new BigDecimal[count][count + 1];
        for (int row = 0; row < count - 1; ++row) {
            for (int col = 0; col < count; ++col) {
                matrix[row][col] = BigDecimal.valueOf(payoff(strategy.first + row, strategy.first + col));
            }
            matrix[row][count] = BigDecimal.ZERO;
        }
        for (int col = 0; col <= count; ++col) matrix[count - 1][col] = BigDecimal.ONE;
        for (int col = 0; col < count; ++col) {
            int pivot = col;
            for (int row = col + 1; row < count; ++row) {
                if (matrix[row][col].abs().compareTo(matrix[pivot][col].abs()) > 0) pivot = row;
            }
            BigDecimal[] swap = matrix[col];
            matrix[col] = matrix[pivot];
            matrix[pivot] = swap;
            require(matrix[col][col].signum() != 0, "dense equilibrium pivot");
            BigDecimal divisor = matrix[col][col];
            for (int q = col; q <= count; ++q) matrix[col][q] = matrix[col][q].divide(divisor, MC);
            for (int row = 0; row < count; ++row) {
                if (row == col) continue;
                BigDecimal multiplier = matrix[row][col];
                for (int q = col; q <= count; ++q) {
                    matrix[row][q] = matrix[row][q].subtract(multiplier.multiply(matrix[col][q], MC), MC);
                }
            }
        }
        BigDecimal[] result = new BigDecimal[count];
        for (int row = 0; row < count; ++row) result[row] = matrix[row][count];
        return result;
    }

    private static boolean close(BigDecimal a, BigDecimal b, BigDecimal tolerance) {
        return a.subtract(b, MC).abs().compareTo(tolerance) < 0;
    }

    private static void runTests(int threadCount) throws InterruptedException, ExecutionException {
        Strategy strategy = new Strategy(1, new BigDecimal[] {BigDecimal.ONE});
        BigDecimal sum = BigDecimal.ZERO;
        BigDecimal roundedTolerance = new BigDecimal("5e-11");
        for (int n = 2; n <= 200; ++n) {
            BigDecimal p = advance(strategy, n);
            if (n >= 3) sum = sum.add(p, MC);
            BigDecimal[] independent = denseEquilibrium(strategy);
            for (int i = 0; i < independent.length; ++i) {
                require(close(independent[i], strategy.probability[i], EPS), "dense payoff-matrix comparison");
            }
            for (int i = 1; i <= n; ++i) {
                require(expectedPayoff(strategy, i).compareTo(EPS) <= 0, "all pure counter-strategies");
            }
            if (n == 3) require(close(p, BigDecimal.ONE.divide(BigDecimal.valueOf(9), MC), EPS), "P(3) = 1/9");
            if (n == 4) require(close(p, BigDecimal.ONE.divide(BigDecimal.valueOf(5), MC), EPS), "P(4) = 1/5");
            if (n == 8) require(p.signum() == 0, "unused final option for n=8");
            if (n == 10) {
                require(close(p, new BigDecimal("0.0479638009"), roundedTolerance), "P(10)");
                require(close(sum, new BigDecimal("1.1546112276"), roundedTolerance), "S(10)");
            }
            if (n == 100) require(close(sum, new BigDecimal("4.8779925686"), roundedTolerance), "S(100)");
        }
        require(close(solve(1000, 1), solve(1000, threadCount), EPS), "thread consistency");
        System.out.println("All checks passed.");
    }

    public static void main(String[] args) {
        try {
            int threadCount = Math.max(1, Runtime.getRuntime().availableProcessors());
            if (args.length == 1 && args[0].equals("--self-test")) {
                runTests(threadCount);
                return;
            }
            if (args.length != 0) throw new IllegalArgumentException("Usage: Euler1012 [--self-test]");
            BigDecimal answer = solve(TARGET, threadCount);
            System.out.println(answer.setScale(10, RoundingMode.HALF_EVEN).toPlainString());
        } catch (IllegalStateException | IllegalArgumentException | InterruptedException | ExecutionException error) {
            System.err.println(error.getMessage());
            System.exit(1);
        }
    }
}

\(n\)	support	\(P(n)\)
\(3\)	\(\{1,2,3\}\)	\(1/9\)
\(4\)	\(\{2,3,4\}\)	\(1/5\)
\(5\)	\(\{3,4,5\}\)	\(5/21\)
\(6\)	\(\{4,5,6\}\)	\(7/27\)
\(7\)	\(\{5,6,7\}\)	\(3/11\)
\(8\)	\(\{5,6,7\}\)	\(0\)
\(9\)	\(\{5,\dots,9\}\)	\(7/275\)
\(10\)	\(\{6,\dots,10\}\)	\(53/1105\)
\(11\)	\(\{7,\dots,11\}\)	\(31/475\)
\(12\)	\(\{8,\dots,12\}\)	\(47/595\)
\(13\)	\(\{9,\dots,13\}\)	\(197/2185\)