Problem 1012: Rock Paper Scissors

Project Euler Problem 1012 Solution

Prove that the unique mixed equilibrium of RPS(n) lives on an odd interval of options ending at n or n−1, derive a closed form for the probability of the largest number and a cubic formula for the points where the support widens, and sum each block of equal width with odd harmonic numbers; for very large N the blocks are expanded around their mean pole with a certified error bound. 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

The equilibrium has a rigid structure. Its support, the set of numbers that are played with positive probability, is always an interval of an odd number $2k+1$ of consecutive numbers ending at $n$ or at $n-1$. On such an interval the probability of the largest number has a closed form, the half-width $k$ changes at points given by a cubic formula, and the sum over all $n$ with the same $k$ reduces to differences of harmonic numbers of odd numbers. For $N=10^5$ only $53$ such blocks occur, with $1431$ terms in total. For much larger $N$ the blocks with $k\ge64$ are summed with an expansion around their mean pole and a certified error bound, which handles $N$ up to $10^{20}$ in seconds.

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 the winner receives $2w-1$. Checking the four cases ($j\lt i$ or $j\gt i$, equal or different parity) gives a single formula for all $i$ and $j$:

$$a_{ij}=(i-j)+\operatorname{sgn}(i-j)\,(-1)^{i+j}\,(i+j-1).$$

For example, if $j\lt i$ and $i-j$ is even, the larger number $i$ wins and the formula gives $(i-j)+(i+j-1)=2i-1$; if $j\lt i$ and $i-j$ is odd, the smaller number $j$ wins and it gives $(i-j)-(i+j-1)=-(2j-1)$. Swapping the players changes the sign, so $a_{ji}=-a_{ij}$: the matrix $A=(a_{ij})$ is skew-symmetric. Against a mixed strategy $p$, the pure counter-strategy $i$ earns

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

A mixed counter-strategy $r$ earns $\sum_i r_if_i\le\max_if_i$, so $p$ is an equilibrium exactly when $f_i\le0$ for every $i$. Moreover $\sum_ip_if_i=p^{\mathsf T}Ap=0$, because $x^{\mathsf T}Ax=0$ for every skew-symmetric matrix. ^[2] A sum of nonpositive terms $p_if_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\ f_i=0,\qquad p_i=0\ \Rightarrow\ f_i\le0.$$

For $n=3$ the payoffs are $f_1=p_2-5p_3$, $f_2=-p_1+3p_3$ and $f_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$. The statement guarantees that every $RPS(n)$ has exactly one equilibrium; we use its existence and its uniqueness below.

2. The support is an odd interval ending at $n$ or $n-1$

Let $X_i$ be the probability that the option $i$ wins, that is, the total probability of the numbers $j\lt i$ of the same parity as $i$ and of the numbers $j\gt i$ of the other parity. Compare the options $i$ and $i+2$, which have the same parity. Against every number $j\notin\{i,i+1,i+2\}$ they obtain the same result, and when both win, $i+2$ receives $4$ dollars more. The three remaining numbers contribute the other terms, and for every probability vector $p$

$$f_{i+2}-f_i=4\,\bigl(X_i-p_{i+1}\bigr)+(2i+3)\,\bigl(p_i+p_{i+2}\bigr)-4i\,p_{i+1},$$

where $X_i-p_{i+1}\ge0$ is the probability of the numbers other than $i+1$ that $i$ beats. Now suppose that $p_i\gt0$ and $p_{i+1}=0$ for some $i\le n-2$. Then $f_i=0$, and the identity gives $f_{i+2}\ge(2i+3)\,p_i\gt0$, which is impossible in an equilibrium. Hence every chosen number below $n-1$ is followed by another chosen number, and the support is an interval $\{a,\dots,\ell\}$ with $\ell=n-1$ or $\ell=n$. A single number $s\ge2$ is not an equilibrium, because $s-1$ beats it, so the interval has at least two elements.

Its size is odd. The submatrix $A_S$ of the payoffs among the chosen numbers satisfies $A_Sp_S=0$ with $p_S\ne0$, so it is singular. All entries of $A$ off the diagonal are odd, so modulo $2$ every principal submatrix of size $m$ equals $J-I\equiv J+I$, where $J$ is the all-ones matrix. For even $m$ we get $(J+I)^2=(m+2)J+I$${}\equiv I\pmod 2$, so the determinant of such a submatrix is odd and in particular nonzero. A support of even size is therefore impossible, and we write its size as $2k+1$. Every skew-symmetric matrix of odd order is singular, ^[2] and the principal submatrices of size $2k$ of $A_S$ are invertible, so $A_S$ has rank exactly $2k$. Its kernel is a line, and together with the normalisation the support determines the equilibrium completely.

3. Two equations on the support

Assume now that the support is $\{a,\dots,n\}$ with $a=n-2k$, and put $L=2k+1$. We use the first moment $M$ and the alternating partial sums $q_i$:

$$M=\sum_{j}j\,p_j,\qquad q_i=\sum_{j=a}^{i}(-1)^{i-j}p_j,\qquad Q=q_n .$$

Thus $q_{a-1}=0$, $p_i=q_i+q_{i-1}$, and $Q$ is the alternating sum of the whole strategy, with the sign $+$ at $n$. The numbers of the parity of $n$ carry the probability $\tfrac{1+Q}{2}$ and the others $\tfrac{1-Q}{2}$.

The first equation comes from the column sums of the support block. Pairing consecutive terms of $\sum_i\operatorname{sgn}(i-j)(-1)^{i+j}(i+j-1)$ gives, for $a\le j\le n$,

$$\sum_{i=a}^{n}a_{ij}=L\Bigl(\frac{a+n}{2}-j\Bigr)+\begin{cases}\phantom{-}k, & j\equiv n\pmod 2,\\ -(k+1), & j\not\equiv n\pmod 2.\end{cases}$$

Adding the $L$ support equations $f_i=0$ therefore gives $L\bigl(\tfrac{a+n}{2}-M\bigr)$${}+k\,\tfrac{1+Q}{2}$${}-(k+1)\,\tfrac{1-Q}{2}=0$, which is the moment equation

$$M-\frac Q2=\frac{a+n}{2}-\frac1{2L}.$$

The second equation comes from adding neighbouring rows. With $\varepsilon_j=(-1)^j$ and $\sum_jp_j=1$, the formula for $a_{ij}$ gives

$$f_i=i-M+\varepsilon_i\bigl[(i-1)\,T_i+U_i\bigr],\qquad T_i=\sum_{j\lt i}\varepsilon_jp_j-\sum_{j\gt i}\varepsilon_jp_j,\qquad U_i=\sum_{j\lt i}\varepsilon_j\,j\,p_j-\sum_{j\gt i}\varepsilon_j\,j\,p_j .$$

Going from $i$ to $i+1$ changes $T_i$ and $U_i$ only in the terms of $p_i$ and $p_{i+1}$, and $\varepsilon_{i+1}=-\varepsilon_i$, so

$$f_i+f_{i+1}=(2i+1)-2M-\varepsilon_iT_i-2i\,p_i+(2i+1)\,p_{i+1}.$$

Since $\varepsilon_iq_i=\sum_{j\le i}\varepsilon_jp_j$, we have $\varepsilon_iT_i=q_i-q_{i-1}-(-1)^{n-i}Q$. Writing $p_i=q_i+q_{i-1}$ and $p_{i+1}=q_{i+1}+q_i$ turns the last line into

$$f_i+f_{i+1}=(2i+1)-2M+(-1)^{n-i}Q+(2i+1)\,q_{i+1}-(2i-1)\,q_{i-1}.$$

This identity holds for every $i$, also outside the support, where $q_{n+1}=-Q$. For $a\le i\le n-1$ both payoffs vanish, which gives the neighbour equations

$$(2i+1)\,q_{i+1}-(2i-1)\,q_{i-1}=2M-(2i+1)-(-1)^{n-i}Q\qquad(a\le i\le n-1).$$

Conversely, the neighbour equations give $f_a=-f_{a+1}=f_{a+2}=\dots=f_n$, and because $L$ is odd, the sum of the support payoffs is then $f_a$. So the moment equation and the $2k$ neighbour equations together are equivalent to the support equations.

4. Solving along the two parity classes

A neighbour equation links $q_{i-1}$ with $q_{i+1}$, so the equations split into two chains: the indices $a-1,a+1,\dots,n-1$, which have the parity of $n-1$, and the indices $a,a+2,\dots,n$. On both chains the homogeneous equation $(2i+1)\,h_{i+1}=(2i-1)\,h_{i-1}$ multiplies by $\tfrac{2i-1}{2i+1}$ per step, and a linear function $\gamma-\tfrac m3$ is a particular solution, because $(2i+1)(i+1)$${}-(2i-1)(i-1)=6i$. Over a whole chain the step factors multiply to

$$R_k(n)=\prod_{j=1}^{k}\frac{2n-4j-1}{2n-4j+1},\qquad R_k(n+1)=\prod_{j=1}^{k}\frac{2n-4j+1}{2n-4j+3},$$

the first one from $a-1$ to $n-1$ and the second one from $a$ to $n$. Their product telescopes: $R_k(n)\,R_k(n+1)=\frac{2a-1}{2n-1}$.

On the chain of $n-1$ the equations have $(-1)^{n-i}=1$, the particular solution is $K-\tfrac m3$ with

$$K=M-\frac{1+Q}{2}=\frac{a+n-1}{2}-\frac1{2L}$$

by the moment equation, and the chain starts at $q_{a-1}=0$. Hence

$$q_{n-1}=K-\frac{n-1}{3}-\Bigl(K-\frac{a-1}{3}\Bigr)R_k(n).$$

On the chain of $n$ the particular solution is $K+Q-\tfrac m3$, and the chain ends at $q_n=Q$. Running it downwards to $q_a=p_a$ and using the telescoping product gives

$$(2a-1)\,p_a=(2a-1)\Bigl(K+Q-\frac a3\Bigr)+(2n-1)\,R_k(n)\Bigl(\frac n3-K\Bigr).$$

The last ingredient is the normalisation. When the neighbour equations for $a\le i\le n-1$ are added, the left sides telescope to $(2n-1)\,Q$${}+(2n-3)\,q_{n-1}$${}-(2a+1)\,q_a$${}-2\sum_{m=a+1}^{n-2}q_m$, the alternating terms $\pm Q$ on the right cancel, and $\sum_{i=a}^{n-1}(2i+1)=n^2-a^2$. The normalisation $1=\sum_jp_j=Q+2\sum_{m=a}^{n-1}q_m$ removes the remaining sum, and with the moment equation the result becomes

$$(2n-1)\,P-(2a-1)\,p_a=(2k-1)\,Q+\frac1L,\qquad P=p_n=Q+q_{n-1}.$$

Substituting $P=Q+q_{n-1}$ and the two chain solutions leaves a linear equation for $Q$, in which $Q$ has the coefficient $(2n-1)-(2a-1)$${}-(2k-1)=L$. With $a=n-2k$ its solution is

$$Q=\frac{(2n-1)\bigl(R_k(n)-1\bigr)+2k}{3}+\frac1L ,$$

and $P=Q+q_{n-1}$ gives the closed form for the probability of the largest number when the support is $\{n-2k,\dots,n\}$:

$$\boxed{\;P_k(n)=\frac{(2k^2+1)\,R_k(n)-2\,(k^2-1)}{3\,(2k+1)}\;}$$

For $k=1$ it reads $P_1(n)=\frac{2n-5}{3(2n-3)}$, which gives $P(3)=\tfrac19$ and $P(4)=\tfrac15$.

5. The options next to the support

The identity for $f_i+f_{i+1}$ also holds next to the support, so it yields the payoffs of the two neighbouring options without any further sums. At $i=n$ we use $q_{n+1}=-Q$ and $f_n=0$; at $i=a-1$ we use $q_{a-2}=q_{a-1}=0$ and $f_a=0$:

$$f_{n+1}=(2n+1)-2M-Q-(2n-1)\,P,\qquad f_{a-1}=(2a-1)\,(1+p_a)-2M-Q .$$

Here $f_{n+1}$ is the payoff of the new option of the larger game $RPS(n+1)$ against the same strategy. Inserting the formulas of Section 4, both payoffs turn out to be multiples of the closed form for the next half-width:

$$f_{a-1}=\frac{2k+3}{2k+1}\,(2n-1)\,R_k(n)\,P_{k+1}(n+1),\qquad f_{n+1}=-\frac{2k+3}{2k+1}\,(2n+1)\,P_{k+1}(n+2).$$

The second formula uses $R_{k+1}(n+2)=\frac{2n-1}{2n+1}R_k(n)$, and the first one $R_{k+1}(n+1)=\frac{2a-3}{(2n-1)\,R_k(n)}$. As $R_k(n)\gt0$ for $a\ge1$, the option $a-1$ loses exactly when $P_{k+1}(n+1)\le0$, and in $RPS(n+1)$ the new option $n+1$ loses exactly when $P_{k+1}(n+2)\ge0$.

6. Where the support widens

For fixed $k$, every factor of $R_k(n)$ increases with $n$, so $P_k(n)$ increases, and it is positive exactly when $R_k(n)\gt\frac{2(k^2-1)}{2k^2+1}$. For $k=1$ this holds for all $n\ge3$, and we set $t_1=3$. For $k\ge2$ let $\alpha_k$ be the root of $P_k$ and $t_k$ the first integer $n$ with $P_k(n)\gt0$. The root can be located through an integral. Put $x=n-k-1$ and let $U$ be the union of the $k$ unit intervals centred at $k-1,k-3,\dots,1-k$:

$$U=\bigcup_{j=1}^{k}\Bigl[k+1-2j-\tfrac12,\ k+1-2j+\tfrac12\Bigr],\qquad G(x)=-\log R_k(n)=\int_U\frac{dt}{x+t}.$$

Indeed, the interval centred at $k+1-2j$ contributes $\log\frac{x+k+\frac32-2j}{x+k+\frac12-2j}=\log\frac{2n-4j+1}{2n-4j-1}$. $G$ decreases, and $P_k(n)\gt0$ exactly when $G(x)\lt\lambda_k=\log\frac{2k^2+1}{2k^2-2}$. Let $x_0=\frac{k(4k^2-1)}{6}$. For the central interval of the same length,

$$\int_{-k/2}^{k/2}\frac{dt}{x_0+t}=\log\frac{x_0+k/2}{x_0-k/2}=\log\frac{k(2k^2+1)/3}{2k(k^2-1)/3}=\lambda_k .$$

Since $U$ is symmetric, pairing $t$ with $-t$ gives $G(x)=\int_{U_+}\frac{2x\,dt}{x^2-t^2}$ with $U_+=U\cap[0,\infty)$, a set of length $k/2$, and the integrand increases with $t$. The increasing map that sends $s\in[0,k/2]$ to the point of $U_+$ with exactly $s$ units of $U_+$ below it never moves a point to the left, and it moves a set of positive length to the right when $k\ge2$. Hence $G(x_0)\gt\lambda_k$. In the other direction, write $y=x_0+\tfrac12$, $b=k-\tfrac12$ and $v=\frac{4k^2-3}{12}$. Every $t\in U$ has $|t|\le b$, and $\int_Ut^2\,dt$${}=\sum_j\bigl((k+1-2j)^2+\tfrac1{12}\bigr)$${}=kv$, so

$$G(y)=\int_U\frac{y\,dt}{y^2-t^2}\le\frac1y\int_U\Bigl(1+\frac{t^2}{y^2-b^2}\Bigr)dt=\frac ky\Bigl(1+\frac{v}{y^2-b^2}\Bigr)\lt\frac k{x_0}\lt\lambda_k .$$

The middle inequality is equivalent to $2x_0v\lt y^2-b^2$${}=x_0^2+x_0-(k^2-k)$, which holds because $x_0\gt2v$ and $x_0\gt k^2-k$ for $k\ge2$. The last one reads $\frac{6}{4k^2-1}\lt\log(1+u)$ with $u=\frac{3}{2k^2-2}$, that is $\frac{2u}{2+u}\lt\log(1+u)$: the difference of the two sides vanishes at $u=0$ and has the derivative $\frac{u^2}{(1+u)(2+u)^2}\gt0$. Since $G$ decreases, its root lies strictly between $x_0$ and $x_0+\tfrac12$, that is

$$\frac{4k^3+5k+6}{6}\lt\alpha_k\lt\frac{4k^3+5k+6}{6}+\frac12 .$$

Because $4k^3+5k=4(k-1)k(k+1)+9k$, the number $\frac{4k^3+5k}{6}$ is an integer for even $k$ and a half-integer for odd $k$. In both cases the interval contains no integer, and the first integer above it is

$$\boxed{\;t_k=\Bigl\lfloor\frac{4k^3+5k}{6}\Bigr\rfloor+2\;}$$

This also gives $t_1=3$; the next values are $9,22,48,89,151,236,350$.

7. Which support belongs to which $n$

Now combine Sections 2, 5 and 6. Let $\{a,\dots,\ell\}$ be the support of the equilibrium of $RPS(n)$, of size $2k+1$.

If $\ell=n$, then $P(n)=P_k(n)\gt0$, so $n\ge t_k$. If $a\ge2$, the option $a-1$ must not gain, so $P_{k+1}(n+1)\le0$, that is $n+1\lt t_{k+1}$. If $a=1$, then $n=2k+1$, which is at most $t_{k+1}-2$ anyway. Hence $t_k\le n$${}\le t_{k+1}-2$.
If $\ell=n-1$, the same strategy satisfies all conditions of $RPS(n-1)$, so by uniqueness it is the equilibrium of $RPS(n-1)$, and the first case gives $n-1\le t_{k+1}-2$. The option $n$ must not gain, so $P_{k+1}(n+1)\ge0$, that is $n+1\ge t_{k+1}$. Hence $n=t_{k+1}-1$.

The ranges $t_k\le n$${}\le t_{k+1}-2$ and the single values $t_{k+1}-1$ cover every $n\ge3$ exactly once, and an equilibrium exists for every $n$. Therefore

$$P(n)=\begin{cases}P_k(n), & t_k\le n\le t_{k+1}-2,\\[2pt] 0, & n=t_{k+1}-1.\end{cases}$$

In words: the support keeps its width $2k+1$ and slides upwards with $n$ through the block $t_k\le n$${}\le t_{k+1}-2$; at $n=t_{k+1}-1$ the new option is not used, and at $n=t_{k+1}$ the support widens by two. The zeros occur at $n=8,21,47,88,150,235,\dots$. Up to $10^5$ there are $52$ of them, and the last of the $53$ blocks is cut off at $N$.

8. Summing a whole block

The closed form still contains a product of $k$ factors, but this product is a rational function of $n$ with simple poles, and its partial fraction decomposition has explicit coefficients. With the central binomial ratios $c_j=\binom{2j}{j}/4^j$, so that $c_0=1$ and $c_j=c_{j-1}\frac{2j-1}{2j}$,

$$R_k(n)=1-\sum_{j=1}^{k}\frac{w_{k,j}}{2n-4j+1},\qquad w_{k,j}=2\,(2k-2j+1)\,c_{j-1}\,c_{k-j}.$$

To see this, write $u_i=2n-4i+1$, so that $R_k(n)=\prod_i\frac{u_i-2}{u_i}$. At the pole $u_j=0$ every other factor equals $\frac{2(j-i)-1}{2(j-i)}$; their product over $i\lt j$ is $c_{j-1}$, over $i\gt j$ it is $(2k-2j+1)\,c_{k-j}$, and the factor $u_j-2$ contributes $-2$. The constant $1$ is the limit for $n\to\infty$. For example, $R_1(n)=1-\frac{2}{2n-3}$ and $R_2(n)=1-\frac{3}{2n-3}-\frac{1}{2n-7}$.

With the harmonic numbers of odd numbers $H^{\mathrm{odd}}_m=\sum_{r=1}^{m}\frac1{2r-1}$, the term $\frac1{2n-4j+1}$ sums over a block $A\le n\le B$ to a difference of two of them, and

$$\sum_{n=A}^{B}P_k(n)=\frac{B-A+1}{2k+1}-\frac{2k^2+1}{3\,(2k+1)}\sum_{j=1}^{k}w_{k,j}\Bigl(H^{\mathrm{odd}}_{B-2j+1}-H^{\mathrm{odd}}_{A-2j}\Bigr).$$

$S(N)$ is the sum of these block sums for $k=1,2,\dots$ with $A=t_k$ and $B=\min\bigl(N,\,t_{k+1}-2\bigr)$, as long as $t_k\le N$. Since $t_k\approx\tfrac23k^3$, there are about $(3N/2)^{1/3}$ blocks, $53$ for $N=10^5$. Finally, the recurrence $\psi(x+1)=\psi(x)+\frac1x$ of the digamma function ^[3] gives $H^{\mathrm{odd}}_m=\tfrac12\bigl(\psi(m+\tfrac12)-\psi(\tfrac12)\bigr)$, so a difference of two odd harmonic numbers costs one evaluation of $\psi$ instead of a long sum.

9. Large blocks: an expansion around the mean pole

The block formula needs $k$ terms for block $k$. That is cheap for $N=10^5$, but for $N=10^{20}$ there are $5\,313\,292$ blocks and about $1.4\cdot10^{13}$ terms. For $k\ge64$ the code therefore evaluates the weighted harmonic sums differently. The weights are positive and add up to $2k$, so $\pi_j=w_{k,j}/(2k)$ is a probability distribution of $X=j-1\in\{0,\dots,k-1\}$. It is the beta-binomial distribution with parameters $k-1$, $\tfrac12$ and $\tfrac32$, whose mean is $\frac{k-1}{4}$. ^[4] With $h=k-1$, $c_n=2n-k-2$ and the centred variable $Z=\frac{4X}{h}-1\in[-1,3]$, every denominator becomes $2n-4j+1=c_n-hZ$, and since $c_n\gt3h$ on these blocks,

$$\sum_{j=1}^{k}\frac{w_{k,j}}{2n-4j+1}=2k\,\mathbb E\Bigl[\frac{1}{c_n-hZ}\Bigr]=2k\sum_{i\ge0}\mu_i\,\frac{h^i}{c_n^{\,i+1}},\qquad\mu_i=\mathbb E\bigl[Z^i\bigr].$$

The moments are polynomials in $1/h$, for example $\mu_0=1$, $\mu_1=0$, $\mu_2=1+\frac2h$, $\mu_3=1+\frac3h+\frac2{h^2}$ and $\mu_4=3+\frac{12}{h}+\frac{10}{h^2}-\frac4{h^3}$; the code stores their integer coefficients up to $\mu_8$. Because $|Z|\le3$, the terms with $i\ge9$ contribute at most $\frac{B-A+1}{a}\cdot\frac{\rho^9}{1-\rho}$ over a block $A\le n\le B$, where $a=2A-k-2$ is the smallest $c_n$ and $\rho=\frac{3h}{a}$. For $k\ge64$ we have $\rho\lt6\cdot10^{-4}$.

The remaining sums of $h^ic^{-(i+1)}$ over the arithmetic progression $c=a,a+2,\dots,b-2$, with $b=2(B+1)-k-2$, follow from the Euler–Maclaurin formula with step $2$: ^[5]

$$\sum_{c=a,\,a+2,\,\dots,\,b-2}f(c)=\frac12\int_a^b f(c)\,dc+\frac{f(a)-f(b)}{2}+\frac{f'(b)-f'(a)}{6}-\frac{f'''(b)-f'''(a)}{90}+\frac{f^{(5)}(b)-f^{(5)}(a)}{945}+R .$$

For $f(c)=h^ic^{-(i+1)}$ the integral is $\frac{h^i}{2i}\bigl(a^{-i}-b^{-i}\bigr)$ for $i\ge1$ and $\frac12\log\frac ba$ for $i=0$. In the remainder of the formula, the derivative $f^{(2m)}$ is multiplied by a periodic Bernoulli function minus $B_{2m}$, which has a constant sign. Every even derivative of $f$ is positive, so consecutive remainders have opposite signs, and since their difference is the next term, the remainder after the $B_4$ term lies between $0$ and the $B_6$ term. The code keeps the $B_6$ term, so $|R|$ is at most its size. It adds this size, weighted by $|\mu_i|$, and the geometric tail to the error bound of the block. With $\Sigma_i$ the resulting value of $\sum_nh^ic_n^{-(i+1)}$, the block sum is

$$\sum_{n=A}^{B}P_k(n)\approx\frac{B-A+1}{2k+1}-\frac{2k\,(2k^2+1)}{3\,(2k+1)}\sum_{i=0}^{8}\mu_i\,\Sigma_i .$$

For $N=10^{20}$ the bounds of the $5\,313\,229$ series blocks add up to about $1.2\cdot10^{-27}$. The first series blocks dominate this total, so it is practically the same already for $N=10^6$, and the program stops with an error if it ever reaches $10^{-24}$. For $N=10^5$ every block has $k\le53$, so the requested sum uses only the exact formula of Section 8.

10. Worked example: $S(10)$

Two blocks contribute to $S(10)$. In the first one, $k=1$ and $3\le n\le7$, the support is $\{n-2,n-1,n\}$ and $P_1(n)=\frac{2n-5}{3(2n-3)}$, which gives $\tfrac19,\tfrac15,\tfrac5{21},\tfrac7{27},\tfrac3{11}$. At $n=8=t_2-1$ the equilibrium of $RPS(7)$ remains: by Section 5 the new option earns $f_8=-\tfrac53\cdot15\cdot P_2(9)=-\tfrac7{11}\lt0$. From $n=9$ on the support has five numbers. For $n=10$ it is $\{6,\dots,10\}$, with $R_2(10)=1-\tfrac3{17}-\tfrac1{13}=\tfrac{165}{221}$ and

$$P(10)=\frac{9\cdot\frac{165}{221}-6}{15}=\frac{53}{1105}\approx0.0479638009 .$$

The whole equilibrium of $RPS(10)$ is

$$\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),$$

with $Q=-\tfrac{79}{1105}$ and $M=\tfrac{1738}{221}$, so that $M-\tfrac Q2=\tfrac{79}{10}=\tfrac{6+10}{2}-\tfrac1{10}$, as the moment equation requires. The option $5$ below the support earns $f_5=-\tfrac{4123}{1105}$. The block formula with $w_{1,1}=2$, $w_{2,1}=3$ and $w_{2,2}=1$ gives

$$\sum_{n=3}^{7}P(n)=\frac53-\frac23\Bigl(\frac13+\frac15+\frac17+\frac19+\frac1{11}\Bigr)=\frac{11239}{10395},$$

$$\sum_{n=9}^{10}P(n)=\frac25-\frac35\Bigl(3\Bigl(\frac1{15}+\frac1{17}\Bigr)+\frac1{11}+\frac1{13}\Bigr)=\frac{4462}{60775},$$

and $S(10)=\tfrac{13262413}{11486475}\approx1.1546112276$, as in the statement. The table lists the first blocks.

$k$	block $t_k\le n$${}\le t_{k+1}-2$	$P(n)=0$ at	$w_{k,j}$
$1$	$3$–$7$	$8$	$2$
$2$	$9$–$20$	$21$	$3$, $1$
$3$	$22$–$46$	$47$	$\tfrac{15}4$, $\tfrac32$, $\tfrac34$
$4$	$48$–$87$	$88$	$\tfrac{35}8$, $\tfrac{15}8$, $\tfrac98$, $\tfrac58$
$5$	$89$–$149$	$150$	$\tfrac{315}{64}$, $\tfrac{35}{16}$, $\tfrac{45}{32}$, $\tfrac{15}{16}$, $\tfrac{35}{64}$
$6$	$151$–$234$	$235$	$\dots$
$53$	$99\,297$–$105\,021$	$105\,022$	$\dots$

For $N=10^5$ the last block is used only for $99\,297\le n\le100\,000$.

How the Code Works

$N$ and the block boundaries are unsigned 128-bit integers, and the program accepts $N\le10^{20}$: parse_limit reads the limit digit by digit and rejects anything else. threshold(k) evaluates $t_k$, and probability(n, k) evaluates $P_k(n)$ directly from the product. block_count(n) finds the number of blocks with $t_k\le n$ by a binary search. evaluate(n, threads) computes the ratios $c_j$ for $j\lt64$ by $c_j=c_{j-1}\frac{2j-1}{2j}$ and distributes chunks of up to $256$ consecutive blocks round-robin over POSIX threads, one per hardware thread. Blocks with $k\lt64$ go to block_sum, the others to series_block_sum. Each block returns a value and an error bound, which is zero for the exact blocks. evaluate adds both chunk by chunk in the order of $k$ and rejects a total bound of $10^{-24}$ or more, and solve returns the value.

block_sum(k, n, central) implements Section 8. It first computes the sum for $j=1$, $\sum_{m=A}^{B}\frac1{2m-3}$${}=\frac12\bigl(\psi(B-\frac12)-\psi(A-\frac32)\bigr)$, with the digamma function of Boost.Math. ^[6] Passing from $j$ to $j+1$ moves the window of odd denominators $2m-4j+1$ down by four, so the code adds the two new terms $\frac1{2A-4j-3}$ and $\frac1{2A-4j-1}$ at the lower end and removes $\frac1{2B-4j-1}$ and $\frac1{2B-4j+1}$ at the upper end. A block therefore costs one digamma evaluation and $O(k)$ further operations. central_moments(k) evaluates $\mu_0,\dots,\mu_8$ by Horner's rule in $1/(k-1)$, and series_block_sum(k, n) implements Section 9 with the Euler–Maclaurin corrections through $B_6$. The arithmetic type is Boost.Multiprecision's cpp_dec_float_50, a radix-10 type with $50$ decimal digits and extra internal guard digits. ^[7] Without arguments the program prints $S(10^5)$ with ten decimals. --limit N prints $S(N)$ instead; --table prints $S(10^e)$ for $e=1,\dots,20$ together with $S(N)/N^{1/3}$, the median time of three runs and the error bound; and --self-test runs the checks listed below.

Precision and the three implementations

The probabilities are rational numbers with quickly growing denominators, so all three versions use $50$-digit decimal arithmetic instead of exact fractions. For $N=10^5$ the results agree to $45$ decimal places with the earlier implementation of this page, which solved one support system for every $n$, so the final sum is far from a rounding boundary at the tenth decimal.

Python uses the decimal module with prec = 50, whose default precision is $28$ digits. ^[8] 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. ^[9] Neither standard library has a digamma function, so both versions implement it. They move the argument above $40$ with $\psi(x)=\psi(x+1)-\frac1x$ and then use the asymptotic expansion $\psi(x)\sim\ln x-\frac1{2x}-\sum_{m\ge1}\frac{B_{2m}}{2m\,x^{2m}}$ through $B_{44}$, whose remainder for real $x\gt0$ is smaller than the first omitted term. ^[10] Java has no logarithm for BigDecimal either; it writes $x=y\cdot2^e$ with $y$ close to $1$ and computes $\ln x=2\operatorname{artanh}\frac{y-1}{y+1}+e\ln2$ from the series of $\operatorname{artanh}$. Java keeps $N$ and the block boundaries in BigInteger, since long ends below $10^{19}$, and runs the chunks on as many threads as there are processors; Python integers are unbounded, and Python evaluates the blocks one after another. All three versions print the same ten decimals, and at $N=10^{20}$ they agree to within $4\cdot10^{-34}$. On a 16-thread machine $N=10^5$ takes a few milliseconds in C++ including start-up, about $0.03$ seconds in Python and about $0.05$ seconds in Java, mostly start-up; $N=10^{16}$ takes $0.4$, $9$ and $1.9$ seconds, and $N=10^{20}$ takes about $8$ seconds in C++, $30$ seconds in Java and $3.5$ minutes in Python.

Complexity and Verification

There are $K\approx(3N/2)^{1/3}$ blocks. The $63$ exact blocks cost at most $\sum_{k\lt64}k=2016$ terms and $63$ digamma evaluations, and every series block costs a fixed number of operations: nine moments, nine Euler–Maclaurin sums and one logarithm. The work is therefore $O(K)=O(N^{1/3})$, and the memory is $O(K)$ at most, one result per chunk. For $N=10^5$ the $53$ exact blocks need $1431$ terms. Solving a support system for every $n$ would instead take $O(N^{4/3})$ operations, because the support has about $(12n)^{1/3}$ elements.

The C++ program runs these checks with --self-test: ^[11]

parse_limit accepts $10^{20}$ and $2^{64}$ and rejects an empty string, signs, an exponent, a trailing space, $10^{20}+1$ and a 41-digit number; block_count returns the right number of blocks for $N=0,2,3,2^{64},10^{20}$; and the 128-bit thresholds of the last block below $10^{20}$ and of the next one agree with unbounded integer arithmetic.
For every block up to the first one beyond $10^5$, the signs $P_k(t_k-1)\lt0\lt P_k(t_k)$ are verified exactly, with the integer numerator $(2k^2+1)\prod(2n-4i-1)$${}-2(k^2-1)\prod(2n-4i+1)$ in cpp_int.
For every $3\le n\le200$, an independent dense Gauss–Jordan solution of the payoff equations on the predicted support, with the last equation replaced by the normalisation, has positive entries summing to $1$; every pure counter-strategy $1\le i\le n$ earns at most $\varepsilon=10^{-40}$; and the probability of $n$ agrees with the closed form to within $\varepsilon$.
The block sums agree with the accumulated dense values of $S(n)$ for $n\le10$ and for $n=21,47,88,100,200$.
$P(3)=\tfrac19$ and $P(4)=\tfrac15$; $P(10)$, $S(10)$ and $S(100)$ agree with the rounded values of the statement to within $5\cdot10^{-11}$.
For $2\le k\le12$, the moments $\mu_0,\dots,\mu_8$ from the coefficient tables agree with the moments computed directly from the weights $w_{k,j}$.
For $k=64$, $65$, $1000$ and $246\,621$, the block that contains $10^{16}$, the series value at the start, the middle and the end of the block differs from the exact block formula by less than its error bound plus $10^{-35}$.
$S(10^5)$, $S(10^6)$ and $S(10^8)$ computed with one thread and with all threads agree to within $\varepsilon$; the last two sums use series blocks, and the third one is split into several chunks per thread. The Java version runs the same list; the Python version, which has no threads, compares the block sums with the direct sum of the closed form for $N=10^4$ instead.

Separate computations checked the derivation itself. In exact rational arithmetic, for every $3\le n\le400$, the predicted support solved by dense elimination is an equilibrium, every option outside it loses strictly, and the closed form, the moment equation, the identity for $f_i+f_{i+1}$ at every $i$, the formulas for $q_{n-1}$, $p_a$, $Q$, $f_{a-1}$ and $f_{n+1}$, and the block sums hold exactly. The partial fractions were checked exactly for $k\le40$, the threshold signs for $k\le60$, and the root bounds of Section 6, including every intermediate inequality, with $60$-digit arithmetic for $2\le k\le400$ and for $k=500,1000,2000,5000$. For Section 9, the normalised weights equal the beta-binomial probabilities exactly for $k\lt60$, the moment polynomials equal the moments of the weights exactly for $k\le150$, and for $k=64,65,70,100,200,1000,3000$ the series value differs from an $80$-digit evaluation of the exact block formula by at most $1.4\%$ of its error bound. For $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 $3\cdot10^{-12}$. ^[12] These comparisons are additional validation, separate from the packaged checks.

Footnotes and References

The references below supply the background facts used above. The payoff formula, the structure of the support, the closed form, the payoffs next to the support, the proof of the transition formula, the block sums and the error bounds of the series 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 and that every skew-symmetric matrix of odd dimension is singular. Sections 1 and 2 use both facts.
NIST DLMF §5.5 — Gamma function, functional relations. Equation 5.5.2 is the recurrence ψ(z+1) = ψ(z) + 1/z, which turns a sum of 1/(r − ½) over consecutive r into a difference of two digamma values.
Wikipedia — Beta-binomial distribution. The article gives the probability mass function of the distribution with parameters n, α and β and its mean nα/(α+β). With n = k − 1, α = 1/2 and β = 3/2 its probabilities are the normalised weights of Section 9, and the mean is (k − 1)/4.
NIST DLMF §2.10 — Sums and sequences. Equation 2.10.1 is the Euler–Maclaurin formula with its integral remainder, in which f⁽²ᵐ⁾ is multiplied by a periodic Bernoulli function; the examples of the section bound a remainder by the next term in the same way as Section 9.
Boost.Math — Digamma. digamma(x) returns ψ(x), the logarithmic derivative of the gamma function; the C++ program calls it with the type cpp_dec_float_50.
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.
NIST DLMF §5.11 — Gamma function, asymptotic expansions. Equation 5.11.2 is the expansion ψ(z) ~ ln z − 1/(2z) − Σ B₂ₖ/(2k z²ᵏ); for real positive z its remainder is bounded in magnitude by the first neglected term and has the same sign.
C++ — Euler1012.cpp. The linked immutable C++ revision contains the exact and the series block sums, the worker threads, the options --limit and --table, 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.
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.

Problem 1012 source code

C++

#include <algorithm>
#include <array>
#include <boost/math/special_functions/digamma.hpp>
#include <boost/multiprecision/cpp_dec_float.hpp>
#include <boost/multiprecision/cpp_int.hpp>
#include <chrono>
#include <cstdint>
#include <cstdlib>
#include <exception>
#include <iomanip>
#include <iostream>
#include <pthread.h>
#include <stdexcept>
#include <string>
#include <thread>
#include <vector>
#ifdef _WIN32
#ifndef NOMINMAX
#define NOMINMAX
#endif
#define WIN32_LEAN_AND_MEAN
#include <windows.h>
#endif

namespace {

using Real = boost::multiprecision::cpp_dec_float_50;
using u64 = std::uint64_t;
using Wide = boost::multiprecision::uint128_t;
using Integer = boost::multiprecision::cpp_int;
constexpr u64 TARGET = 100'000;
const Wide MAX_LIMIT = Wide(10'000'000'000ULL) * 10'000'000'000ULL;
constexpr unsigned SERIES_START = 64;
constexpr unsigned SERIES_ORDER = 8;
const Real EPS("1e-40");

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

Wide parse_limit(const std::string& text) {
    require(!text.empty(), "integer limit");
    Wide n = 0;
    for (const char digit : text) {
        require(digit >= '0' && digit <= '9', "integer limit");
        n = 10 * n + (digit - '0');
        require(n <= MAX_LIMIT, "supported limit through 10^20");
    }
    return n;
}

Wide threshold(const u64 k) {
    return (Wide(4) * k * k * k + 5 * k) / 6 + 2;
}

unsigned block_count(const Wide& n) {
    unsigned low = 0, high = 1;
    while (threshold(high) <= n) high *= 2;
    while (high - low > 1) {
        const unsigned middle = low + (high - low) / 2;
        if (threshold(middle) <= n) low = middle;
        else high = middle;
    }
    return low;
}

Real probability(const Wide& n, const u64 k) {
    Real product = 1;
    for (u64 i = 1; i <= k; ++i) product *= Real(2 * n - 4 * i - 1) / Real(2 * n - 4 * i + 1);
    return ((2 * k * k + 1) * product - 2 * (k * k - 1)) / (3 * (2 * k + 1));
}

Real block_sum(const u64 k, const Wide& n, const std::vector<Real>& central) {
    const Wide first = threshold(k);
    const Wide last = std::min(n, Wide(threshold(k + 1) - 2));
    if (last < first) return 0;

    Real harmonic = (boost::math::digamma(Real(last) - Real("0.5"))
        - boost::math::digamma(Real(first) - Real("1.5"))) / 2;
    Real weighted = 0;
    for (u64 i = 1; i <= k; ++i) {
        const Real residue = 2 * central[i - 1] * (2 * k - 2 * i + 1) * central[k - i];
        weighted += residue * harmonic;
        if (i < k) {
            harmonic += Real(1) / Real(2 * first - 4 * i - 3) + Real(1) / Real(2 * first - 4 * i - 1)
                - Real(1) / Real(2 * last - 4 * i - 1) - Real(1) / Real(2 * last - 4 * i + 1);
        }
    }
    return (Real(last - first + 1) - Real(2 * k * k + 1) * weighted / 3) / (2 * k + 1);
}

struct Result {
    Real value = 0;
    Real error_bound = 0;
};

std::array<Real, SERIES_ORDER + 1> central_moments(const u64 k) {
    constexpr int coefficients[SERIES_ORDER + 1][SERIES_ORDER] = {
        {1}, {0}, {1, 2}, {1, 3, 2}, {3, 12, 10, -4},
        {6, 30, 40, 0, -16}, {15, 90, 150, 0, -104, 32},
        {36, 252, 525, 105, -560, -84, 272},
        {91, 728, 1792, 560, -2618, -616, 2248, -544}
    };
    const Real inverse = Real(1) / (k - 1);
    std::array<Real, SERIES_ORDER + 1> moments{};
    moments[0] = 1;
    for (unsigned j = 2; j <= SERIES_ORDER; ++j) {
        for (int degree = static_cast<int>(j) - 1; degree >= 0; --degree) {
            moments[j] = moments[j] * inverse + coefficients[j][degree];
        }
    }
    return moments;
}

Result series_block_sum(const u64 k, const Wide& n) {
    const Wide first = threshold(k);
    const Wide last = std::min(n, Wide(threshold(k + 1) - 2));
    if (last < first) return {};
    const Real h = k - 1;
    const Real a = 2 * Real(first) - k - 2;
    const Real b = 2 * Real(last + 1) - k - 2;
    const Real ia = 1 / a, ib = 1 / b;
    const Real ua = h * ia, ub = h * ib;
    const Real ia2 = ia * ia, ib2 = ib * ib;
    const Real ia4 = ia2 * ia2, ib4 = ib2 * ib2;
    const Real ia6 = ia4 * ia2, ib6 = ib4 * ib2;
    const auto moments = central_moments(k);
    Real pa = 1, pb = 1, weighted = 0, bound = 0;
    for (unsigned j = 0; j <= SERIES_ORDER; ++j) {
        const Real integral = j == 0 ? Real(log(b / a) / 2) : Real((pa - pb) / (2 * j));
        const Real second = Real(j + 1) / 6;
        const Real fourth = Real((j + 1) * (j + 2) * (j + 3)) / 90;
        const Real sixth = Real((j + 1) * (j + 2) * (j + 3) * (j + 4) * (j + 5)) / 945;
        const Real remainder = sixth * (pa * ia6 - pb * ib6);
        const Real harmonic = integral + (pa * ia - pb * ib) / 2
            + second * (pa * ia2 - pb * ib2) - fourth * (pa * ia4 - pb * ib4) + remainder;
        weighted += moments[j] * harmonic;
        bound += abs(moments[j]) * remainder;
        pa *= ua;
        pb *= ub;
    }
    const Real rho = 3 * h / a;
    bound += Real(last - first + 1) / a * pow(rho, SERIES_ORDER + 1) / (1 - rho);
    const Real coefficient = Real(2 * k * k + 1) * (2 * k) / (3 * (2 * k + 1));
    return {Real(last - first + 1) / (2 * k + 1) - coefficient * weighted, coefficient * bound};
}

struct Task {
    Wide n = 0;
    unsigned index = 0;
    unsigned stride = 1;
    unsigned blocks = 0;
    unsigned chunk_size = 1;
    const std::vector<Real>* central = nullptr;
    std::vector<Result>* sums = nullptr;
    std::exception_ptr error;
};

void* sum_worker(void* argument) {
    Task& task = *static_cast<Task*>(argument);
    try {
        for (unsigned chunk = task.index; chunk < task.sums->size(); chunk += task.stride) {
            const unsigned first = chunk * task.chunk_size + 1;
            const unsigned last = std::min(task.blocks, first + task.chunk_size - 1);
            Result sum;
            for (unsigned k = first; k <= last; ++k) {
                const Result current = k < SERIES_START
                    ? Result{block_sum(k, task.n, *task.central), 0}
                    : series_block_sum(k, task.n);
                sum.value += current.value;
                sum.error_bound += current.error_bound;
            }
            (*task.sums)[chunk] = sum;
        }
    } catch (...) {
        task.error = std::current_exception();
    }
    return nullptr;
}

Result evaluate(const Wide& n, unsigned thread_count) {
    require(n <= MAX_LIMIT, "supported limit through 10^20");
    if (n < 3) return {};
    const unsigned blocks = block_count(n);
    thread_count = std::min(thread_count, blocks);
    require(thread_count > 0, "positive thread count");
    const unsigned chunk_size = std::min(256U, std::max(1U, blocks / (4 * thread_count)));
    std::vector<Real> central(std::min(blocks, SERIES_START));
    std::vector<Result> sums((blocks + chunk_size - 1) / chunk_size);
    central[0] = 1;
    for (unsigned j = 1; j < central.size(); ++j) central[j] = central[j - 1] * (2 * j - 1) / (2 * j);
    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] = {n, t, thread_count, blocks, chunk_size, &central, &sums, {}};
        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");
    for (const Task& task : tasks) if (task.error) std::rethrow_exception(task.error);
    Result result;
    for (const Result& sum : sums) {
        result.value += sum.value;
        result.error_bound += sum.error_bound;
    }
    require(result.error_bound < Real("1e-24"), "series truncation bound");
    return result;
}

Real solve(const Wide& n, const unsigned thread_count) {
    return evaluate(n, thread_count).value;
}

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

std::vector<Real> dense_equilibrium(const int first, const int last) {
    const int count = last - first + 1;
    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(first + row, 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;
}

Integer threshold_numerator(const u64 k, const Wide& n) {
    Integer numerator = 1, denominator = 1;
    for (u64 i = 1; i <= k; ++i) {
        numerator *= 2 * Integer(n) - 4 * i - 1;
        denominator *= 2 * Integer(n) - 4 * i + 1;
    }
    return (2 * k * k + 1) * numerator - 2 * (k * k - 1) * denominator;
}

void run_tests(const unsigned thread_count) {
    require(parse_limit("100000000000000000000") == MAX_LIMIT, "parse 10^20");
    require(parse_limit("18446744073709551616") == (Wide(1) << 64), "parse above 64 bits");
    for (const std::string text : {"", "-1", "+1", "1e20", "1 ", "100000000000000000001",
            "99999999999999999999999999999999999999999"}) {
        bool rejected = false;
        try { parse_limit(text); }
        catch (const std::runtime_error&) { rejected = true; }
        require(rejected, "reject invalid or unsupported limit");
    }
    for (const Wide& n : {Wide(0), Wide(2), Wide(3), Wide(1) << 64, MAX_LIMIT}) {
        const unsigned count = block_count(n);
        require((count == 0 || threshold(count) <= n) && threshold(count + 1) > n,
            "block search boundaries");
    }
    const u64 largest_block = block_count(MAX_LIMIT);
    for (const u64 block : {largest_block, largest_block + 1}) {
        require(Integer(threshold(block)) == (4 * Integer(block) * block * block + 5 * block) / 6 + 2,
            "wide cubic threshold versus unbounded integer");
    }
    int max_block = 1;
    while (threshold(max_block) <= TARGET) ++max_block;
    for (int k = 1; k <= max_block; ++k) {
        const int first = static_cast<int>(threshold(k));
        require(threshold_numerator(k, first - 1) < 0, "exact lower threshold sign");
        require(threshold_numerator(k, first) > 0, "exact upper threshold sign");
    }
    Real sum = 0;
    int k = 1;
    for (int n = 3; n <= 200; ++n) {
        if (threshold(k + 1) <= static_cast<u64>(n)) ++k;
        const int last = static_cast<u64>(n + 1) == threshold(k + 1) ? n - 1 : n;
        const int first = last - 2 * k;
        const auto p = dense_equilibrium(first, last);
        Real mass = 0;
        for (const Real& value : p) {
            require(value > 0, "positive active probabilities");
            mass += value;
        }
        require(abs(mass - 1) < EPS, "normalized dense equilibrium");
        for (int i = 1; i <= n; ++i) {
            Real value = 0;
            for (int j = first; j <= last; ++j) value += payoff(i, j) * p[j - first];
            require(value <= EPS, "all pure counter-strategies");
        }
        const Real final_probability = last == n ? p.back() : Real(0);
        require(abs(final_probability - (last == n ? probability(n, k) : Real(0))) < EPS,
            "closed probability versus dense payoff matrix");
        sum += final_probability;
        if (n <= 10 || n == 21 || n == 47 || n == 88 || n == 100 || n == 200) {
            require(abs(sum - solve(n, 1)) < EPS, "harmonic block sum versus dense equilibria");
        }
    }
    require(abs(probability(3, 1) - Real(1) / 9) < EPS, "P(3) = 1/9");
    require(abs(probability(4, 1) - Real(1) / 5) < EPS, "P(4) = 1/5");
    require(abs(probability(10, 2) - Real("0.0479638009")) < Real("5e-11"), "P(10)");
    require(abs(solve(10, 1) - Real("1.1546112276")) < Real("5e-11"), "S(10)");
    require(abs(solve(100, 1) - Real("4.8779925686")) < Real("5e-11"), "S(100)");
    require(abs(solve(TARGET, 1) - solve(TARGET, thread_count)) < EPS, "thread consistency");
    for (u64 block = 2; block <= 12; ++block) {
        std::vector<Real> central(block);
        central[0] = 1;
        for (u64 j = 1; j < block; ++j) central[j] = central[j - 1] * (2 * j - 1) / (2 * j);
        std::array<Real, SERIES_ORDER + 1> direct{};
        for (u64 j = 0; j < block; ++j) {
            const Real mass = central[j] * (2 * block - 2 * j - 1) * central[block - j - 1] / block;
            const Real z = 4 * Real(j) / (block - 1) - 1;
            Real power = 1;
            for (unsigned order = 0; order <= SERIES_ORDER; ++order) {
                direct[order] += mass * power;
                power *= z;
            }
        }
        const auto moments = central_moments(block);
        for (unsigned order = 0; order <= SERIES_ORDER; ++order) {
            require(abs(direct[order] - moments[order]) < EPS, "central moments versus partial-fraction residues");
        }
    }
    const u64 reference_block = block_count(Wide(10'000'000'000'000'000ULL));
    for (const u64 block : std::array<u64, 4>{64, 65, 1000, reference_block}) {
        std::vector<Real> central(block);
        central[0] = 1;
        for (u64 j = 1; j < block; ++j) central[j] = central[j - 1] * (2 * j - 1) / (2 * j);
        const Wide first = threshold(block), last = threshold(block + 1) - 2;
        for (const Wide& end : {first, Wide(first + (last - first) / 2), last}) {
            const Result approximation = series_block_sum(block, end);
            require(abs(approximation.value - block_sum(block, end, central))
                < approximation.error_bound + Real("1e-35"), "bounded series versus exact harmonic block");
        }
    }
    require(abs(solve(1'000'000, 1) - solve(1'000'000, thread_count)) < EPS,
        "series thread consistency");
    require(abs(solve(100'000'000, 1) - solve(100'000'000, thread_count)) < EPS,
        "chunked thread consistency");
    std::cout << "All checks passed.\n";
}

unsigned logical_processor_count() {
#ifdef _WIN32
    return static_cast<unsigned>(GetActiveProcessorCount(ALL_PROCESSOR_GROUPS));
#else
    return std::thread::hardware_concurrency();
#endif
}

void print_table(const unsigned thread_count) {
    std::cout << "Hardware threads: " << thread_count << "; three runs per input.\n";
    std::cout << "| N | S(N) | S(N)/N^(1/3) | Median seconds | Series error bound |\n"
              << "|---:|---:|---:|---:|---:|\n";
    solve(TARGET, thread_count);
    Wide n = 1;
    for (unsigned exponent = 1; exponent <= 20; ++exponent) {
        n *= 10;
        std::array<double, 3> seconds{};
        Result result;
        for (unsigned repeat = 0; repeat < seconds.size(); ++repeat) {
            const auto start = std::chrono::steady_clock::now();
            const Result current = evaluate(n, thread_count);
            seconds[repeat] = std::chrono::duration<double>(std::chrono::steady_clock::now() - start).count();
            if (repeat != 0) require(abs(current.value - result.value) < EPS, "table repeat consistency");
            result = current;
        }
        std::sort(seconds.begin(), seconds.end());
        std::cout << "| 10^" << exponent << " | " << std::fixed << std::setprecision(24)
                  << result.value << " | " << std::setprecision(23) << result.value / pow(Real(n), Real(1) / 3)
                  << " | " << std::setprecision(6) << seconds[1] << " | "
                  << std::scientific << std::setprecision(2) << result.error_bound << " |" << std::endl;
    }
}

}

int main(int argc, char* argv[]) {
    try {
        const unsigned thread_count = std::max(1U, logical_processor_count());
        if (argc == 2 && std::string(argv[1]) == "--self-test") {
            run_tests(thread_count);
            return EXIT_SUCCESS;
        }
        if (argc == 2 && std::string(argv[1]) == "--table") {
            print_table(thread_count);
            return EXIT_SUCCESS;
        }
        Wide n = TARGET;
        if (argc == 3 && std::string(argv[1]) == "--limit") {
            n = parse_limit(argv[2]);
        } else if (argc != 1) {
            throw std::invalid_argument("Usage: Euler1012 [--self-test | --table | --limit N]");
        }
        std::cout << std::fixed << std::setprecision(10) << solve(n, 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 sys
import time
from decimal import Decimal, getcontext
from fractions import Fraction
from math import comb

getcontext().prec = 50

TARGET = 100_000
MAX_LIMIT = 10 ** 20
SERIES_START = 64
SERIES_ORDER = 8
EPS = Decimal("1e-40")
ZERO = Decimal(0)
ONE = Decimal(1)
HALF = Decimal("0.5")


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


def bernoulli_numbers(count):
    numbers = [Fraction(1)]
    for m in range(1, count + 1):
        numbers.append(-sum(comb(m + 1, j) * numbers[j] for j in range(m)) / (m + 1))
    return numbers


DIGAMMA_SHIFT = 40
DIGAMMA_TERMS = [Decimal(b.numerator) / b.denominator / (2 * m)
                 for m, b in enumerate(bernoulli_numbers(44)[2::2], 1)]


def digamma(x):
    correction = ZERO
    while x < DIGAMMA_SHIFT:
        correction += ONE / x
        x += 1
    inverse_square = ONE / (x * x)
    power = inverse_square
    series = ZERO
    for term in DIGAMMA_TERMS:
        series += term * power
        power *= inverse_square
    return x.ln() - ONE / (2 * x) - series - correction


def parse_limit(text):
    require(text != "", "integer limit")
    n = 0
    for digit in text:
        require("0" <= digit <= "9", "integer limit")
        n = 10 * n + ord(digit) - ord("0")
        require(n <= MAX_LIMIT, "supported limit through 10^20")
    return n


def threshold(k):
    return (4 * k ** 3 + 5 * k) // 6 + 2


def block_count(n):
    low, high = 0, 1
    while threshold(high) <= n:
        high *= 2
    while high - low > 1:
        middle = low + (high - low) // 2
        if threshold(middle) <= n:
            low = middle
        else:
            high = middle
    return low


def probability(n, k):
    product = ONE
    for i in range(1, k + 1):
        product *= Decimal(2 * n - 4 * i - 1) / (2 * n - 4 * i + 1)
    return ((2 * k * k + 1) * product - 2 * (k * k - 1)) / (3 * (2 * k + 1))


def block_sum(k, n, central):
    first = threshold(k)
    last = min(n, threshold(k + 1) - 2)
    if last < first:
        return ZERO

    harmonic = (digamma(Decimal(last) - HALF) - digamma(Decimal(first) - Decimal("1.5"))) / 2
    weighted = ZERO
    for i in range(1, k + 1):
        residue = 2 * central[i - 1] * (2 * k - 2 * i + 1) * central[k - i]
        weighted += residue * harmonic
        if i < k:
            harmonic += (ONE / (2 * first - 4 * i - 3) + ONE / (2 * first - 4 * i - 1)
                         - ONE / (2 * last - 4 * i - 1) - ONE / (2 * last - 4 * i + 1))
    return (Decimal(last - first + 1) - (2 * k * k + 1) * weighted / 3) / (2 * k + 1)


MOMENT_COEFFICIENTS = [
    [1], [0], [1, 2], [1, 3, 2], [3, 12, 10, -4],
    [6, 30, 40, 0, -16], [15, 90, 150, 0, -104, 32],
    [36, 252, 525, 105, -560, -84, 272],
    [91, 728, 1792, 560, -2618, -616, 2248, -544],
]


def central_moments(k):
    inverse = ONE / (k - 1)
    moments = [ONE] + [ZERO] * SERIES_ORDER
    for j in range(2, SERIES_ORDER + 1):
        for coefficient in reversed(MOMENT_COEFFICIENTS[j]):
            moments[j] = moments[j] * inverse + coefficient
    return moments


def series_block_sum(k, n):
    first = threshold(k)
    last = min(n, threshold(k + 1) - 2)
    if last < first:
        return ZERO, ZERO
    h = Decimal(k - 1)
    a = Decimal(2 * first - k - 2)
    b = Decimal(2 * (last + 1) - k - 2)
    ia, ib = ONE / a, ONE / b
    ua, ub = h * ia, h * ib
    ia2, ib2 = ia * ia, ib * ib
    ia4, ib4 = ia2 * ia2, ib2 * ib2
    ia6, ib6 = ia4 * ia2, ib4 * ib2
    moments = central_moments(k)
    pa = pb = ONE
    weighted = bound = ZERO
    for j in range(SERIES_ORDER + 1):
        integral = (b / a).ln() / 2 if j == 0 else (pa - pb) / (2 * j)
        second = Decimal(j + 1) / 6
        fourth = Decimal((j + 1) * (j + 2) * (j + 3)) / 90
        sixth = Decimal((j + 1) * (j + 2) * (j + 3) * (j + 4) * (j + 5)) / 945
        remainder = sixth * (pa * ia6 - pb * ib6)
        harmonic = (integral + (pa * ia - pb * ib) / 2
                    + second * (pa * ia2 - pb * ib2) - fourth * (pa * ia4 - pb * ib4) + remainder)
        weighted += moments[j] * harmonic
        bound += abs(moments[j]) * remainder
        pa *= ua
        pb *= ub
    rho = 3 * h / a
    bound += Decimal(last - first + 1) / a * rho ** (SERIES_ORDER + 1) / (1 - rho)
    coefficient = Decimal((2 * k * k + 1) * 2 * k) / (3 * (2 * k + 1))
    return Decimal(last - first + 1) / (2 * k + 1) - coefficient * weighted, coefficient * bound


def evaluate(n):
    require(n <= MAX_LIMIT, "supported limit through 10^20")
    if n < 3:
        return ZERO, ZERO
    blocks = block_count(n)
    central = [ONE]
    for j in range(1, min(blocks, SERIES_START)):
        central.append(central[-1] * (2 * j - 1) / (2 * j))
    value = bound = ZERO
    for k in range(1, blocks + 1):
        if k < SERIES_START:
            value += block_sum(k, n, central)
        else:
            part, error = series_block_sum(k, n)
            value += part
            bound += error
    require(bound < Decimal("1e-24"), "series truncation bound")
    return value, bound


def solve(n):
    return evaluate(n)[0]


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)


def dense_equilibrium(first, last):
    count = last - first + 1
    matrix = [[Decimal(payoff(first + row, 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 threshold_numerator(k, n):
    numerator = denominator = 1
    for i in range(1, k + 1):
        numerator *= 2 * n - 4 * i - 1
        denominator *= 2 * n - 4 * i + 1
    return (2 * k * k + 1) * numerator - 2 * (k * k - 1) * denominator


def direct_sum(n):
    result = ZERO
    k = 1
    for m in range(3, n + 1):
        if threshold(k + 1) <= m:
            k += 1
        if m + 1 != threshold(k + 1):
            result += probability(m, k)
    return result


def central_ratios(count):
    central = [ONE]
    for j in range(1, count):
        central.append(central[-1] * (2 * j - 1) / (2 * j))
    return central


def run_tests():
    require(parse_limit("100000000000000000000") == MAX_LIMIT, "parse 10^20")
    require(parse_limit("18446744073709551616") == 2 ** 64, "parse above 64 bits")
    for text in ("", "-1", "+1", "1e20", "1 ", "100000000000000000001", "99999999999999999999999999999999999999999"):
        rejected = False
        try:
            parse_limit(text)
        except AssertionError:
            rejected = True
        require(rejected, "reject invalid or unsupported limit")
    for n in (0, 2, 3, 2 ** 64, MAX_LIMIT):
        count = block_count(n)
        require((count == 0 or threshold(count) <= n) and threshold(count + 1) > n, "block search boundaries")
    largest_block = block_count(MAX_LIMIT)
    for block in (largest_block, largest_block + 1):
        require(threshold(block) == (4 * block * block * block + 5 * block) // 6 + 2,
                "wide cubic threshold versus unbounded integer")
    max_block = 1
    while threshold(max_block) <= TARGET:
        max_block += 1
    for k in range(1, max_block + 1):
        first = threshold(k)
        require(threshold_numerator(k, first - 1) < 0, "exact lower threshold sign")
        require(threshold_numerator(k, first) > 0, "exact upper threshold sign")
    total = ZERO
    k = 1
    for n in range(3, 201):
        if threshold(k + 1) <= n:
            k += 1
        last = n - 1 if n + 1 == threshold(k + 1) else n
        first = last - 2 * k
        p = dense_equilibrium(first, last)
        require(all(value > 0 for value in p), "positive active probabilities")
        require(abs(sum(p) - 1) < EPS, "normalized dense equilibrium")
        for i in range(1, n + 1):
            value = sum(payoff(i, j) * p[j - first] for j in range(first, last + 1))
            require(value <= EPS, "all pure counter-strategies")
        final_probability = p[-1] if last == n else ZERO
        require(abs(final_probability - (probability(n, k) if last == n else ZERO)) < EPS,
                "closed probability versus dense payoff matrix")
        total += final_probability
        if n <= 10 or n in (21, 47, 88, 100, 200):
            require(abs(total - solve(n)) < EPS, "harmonic block sum versus dense equilibria")
    rounded_tolerance = Decimal("5e-11")
    require(abs(probability(3, 1) - ONE / 9) < EPS, "P(3) = 1/9")
    require(abs(probability(4, 1) - ONE / 5) < EPS, "P(4) = 1/5")
    require(abs(probability(10, 2) - Decimal("0.0479638009")) < rounded_tolerance, "P(10)")
    require(abs(solve(10) - Decimal("1.1546112276")) < rounded_tolerance, "S(10)")
    require(abs(solve(100) - Decimal("4.8779925686")) < rounded_tolerance, "S(100)")
    require(abs(solve(10_000) - direct_sum(10_000)) < EPS, "harmonic block sum versus direct closed-form sum")
    for block in range(2, 13):
        central = central_ratios(block)
        direct = [ZERO] * (SERIES_ORDER + 1)
        for j in range(block):
            mass = central[j] * (2 * block - 2 * j - 1) * central[block - j - 1] / block
            z = 4 * Decimal(j) / (block - 1) - 1
            power = ONE
            for order in range(SERIES_ORDER + 1):
                direct[order] += mass * power
                power *= z
        moments = central_moments(block)
        for order in range(SERIES_ORDER + 1):
            require(abs(direct[order] - moments[order]) < EPS, "central moments versus partial-fraction residues")
    reference_block = block_count(10 ** 16)
    for block in (64, 65, 1000, reference_block):
        central = central_ratios(block)
        first, last = threshold(block), threshold(block + 1) - 2
        for end in (first, first + (last - first) // 2, last):
            value, error = series_block_sum(block, end)
            require(abs(value - block_sum(block, end, central)) < error + Decimal("1e-35"),
                    "bounded series versus exact harmonic block")
    print("All checks passed.")


def print_table():
    print("Python evaluates the blocks one after another; one run per input.")
    print("| N | S(N) | S(N)/N^(1/3) | Seconds | Series error bound |")
    print("|---:|---:|---:|---:|---:|")
    n = 1
    for exponent in range(1, 21):
        n *= 10
        start = time.perf_counter()
        value, bound = evaluate(n)
        seconds = time.perf_counter() - start
        ratio = value / Decimal(n) ** (ONE / 3)
        print(f"| 10^{exponent} | {value:.24f} | {ratio:.23f} | {seconds:.6f} | {float(bound):.2e} |", flush=True)


def main():
    if len(sys.argv) == 2 and sys.argv[1] == "--self-test":
        run_tests()
        return
    if len(sys.argv) == 2 and sys.argv[1] == "--table":
        print_table()
        return
    n = TARGET
    if len(sys.argv) == 3 and sys.argv[1] == "--limit":
        n = parse_limit(sys.argv[2])
    elif len(sys.argv) != 1:
        raise AssertionError("Usage: Euler1012.py [--self-test | --table | --limit N]")
    print(f"{solve(n):.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.BigInteger;
import java.math.MathContext;
import java.math.RoundingMode;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;
import java.util.Locale;
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 long TARGET = 100_000L;
    private static final BigInteger MAX_LIMIT = BigInteger.TEN.pow(20);
    private static final int SERIES_START = 64;
    private static final int SERIES_ORDER = 8;
    private static final MathContext MC = new MathContext(50, RoundingMode.HALF_EVEN);
    private static final MathContext WIDE = new MathContext(70, RoundingMode.HALF_EVEN);
    private static final BigDecimal EPS = new BigDecimal("1e-40");
    private static final BigDecimal TWO = BigDecimal.valueOf(2);
    private static final BigDecimal THREE = BigDecimal.valueOf(3);
    private static final BigDecimal HALF = new BigDecimal("0.5");
    private static final BigDecimal ONE_AND_HALF = new BigDecimal("1.5");
    private static final BigDecimal LN2 = twiceAtanh(BigDecimal.ONE.divide(THREE, WIDE));
    private static final String[] BERNOULLI = {
        "1/6", "-1/30", "1/42", "-1/30", "5/66", "-691/2730", "7/6", "-3617/510", "43867/798", "-174611/330",
        "854513/138", "-236364091/2730", "8553103/6", "-23749461029/870", "8615841276005/14322",
        "-7709321041217/510", "2577687858367/6", "-26315271553053477373/1919190", "2929993913841559/6",
        "-261082718496449122051/13530", "1520097643918070802691/1806", "-27833269579301024235023/690"
    };
    private static final int DIGAMMA_SHIFT = 40;
    private static final BigDecimal[] DIGAMMA_TERMS = digammaTerms();
    private static final int[][] MOMENT_COEFFICIENTS = {
        {1}, {0}, {1, 2}, {1, 3, 2}, {3, 12, 10, -4},
        {6, 30, 40, 0, -16}, {15, 90, 150, 0, -104, 32},
        {36, 252, 525, 105, -560, -84, 272},
        {91, 728, 1792, 560, -2618, -616, 2248, -544}
    };

    private static final class Result {
        final BigDecimal value;
        final BigDecimal errorBound;

        Result(BigDecimal value, BigDecimal errorBound) {
            this.value = value;
            this.errorBound = errorBound;
        }
    }

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

    private static BigDecimal twiceAtanh(BigDecimal z) {
        BigDecimal square = z.multiply(z, WIDE);
        BigDecimal power = z;
        BigDecimal sum = BigDecimal.ZERO;
        BigDecimal negligible = new BigDecimal("1e-75");
        for (int m = 1; power.abs().compareTo(negligible) > 0; m += 2) {
            sum = sum.add(power.divide(BigDecimal.valueOf(m), WIDE), WIDE);
            power = power.multiply(square, WIDE);
        }
        return sum.multiply(TWO, WIDE);
    }

    private static BigDecimal ln(BigDecimal x) {
        require(x.signum() > 0, "positive logarithm argument");
        int exponent = (int) Math.round(Math.log(x.doubleValue()) / Math.log(2));
        BigDecimal y = exponent >= 0 ? x.divide(TWO.pow(exponent), WIDE) : x.multiply(TWO.pow(-exponent), WIDE);
        BigDecimal z = y.subtract(BigDecimal.ONE).divide(y.add(BigDecimal.ONE), WIDE);
        return twiceAtanh(z).add(LN2.multiply(BigDecimal.valueOf(exponent), WIDE), WIDE).round(MC);
    }

    private static BigDecimal[] digammaTerms() {
        BigDecimal[] terms = new BigDecimal[BERNOULLI.length];
        for (int m = 1; m <= BERNOULLI.length; ++m) {
            String[] fraction = BERNOULLI[m - 1].split("/");
            BigDecimal denominator = new BigDecimal(fraction[1]).multiply(BigDecimal.valueOf(2L * m));
            terms[m - 1] = new BigDecimal(fraction[0]).divide(denominator, MC);
        }
        return terms;
    }

    private static BigDecimal digamma(BigDecimal x) {
        BigDecimal correction = BigDecimal.ZERO;
        BigDecimal shift = BigDecimal.valueOf(DIGAMMA_SHIFT);
        while (x.compareTo(shift) < 0) {
            correction = correction.add(BigDecimal.ONE.divide(x, MC), MC);
            x = x.add(BigDecimal.ONE);
        }
        BigDecimal inverseSquare = BigDecimal.ONE.divide(x.multiply(x, MC), MC);
        BigDecimal power = inverseSquare;
        BigDecimal series = BigDecimal.ZERO;
        for (BigDecimal term : DIGAMMA_TERMS) {
            series = series.add(term.multiply(power, MC), MC);
            power = power.multiply(inverseSquare, MC);
        }
        return ln(x).subtract(BigDecimal.ONE.divide(x.multiply(TWO), MC), MC).subtract(series, MC).subtract(correction, MC);
    }

    private static BigInteger parseLimit(String text) {
        require(!text.isEmpty(), "integer limit");
        BigInteger n = BigInteger.ZERO;
        for (char digit : text.toCharArray()) {
            require(digit >= '0' && digit <= '9', "integer limit");
            n = n.multiply(BigInteger.TEN).add(BigInteger.valueOf(digit - '0'));
            require(n.compareTo(MAX_LIMIT) <= 0, "supported limit through 10^20");
        }
        return n;
    }

    private static BigInteger threshold(long k) {
        BigInteger block = BigInteger.valueOf(k);
        return block.pow(3).shiftLeft(2).add(block.multiply(BigInteger.valueOf(5))).divide(BigInteger.valueOf(6)).add(BigInteger.TWO);
    }

    private static int blockCount(BigInteger n) {
        int low = 0;
        int high = 1;
        while (threshold(high).compareTo(n) <= 0) high *= 2;
        while (high - low > 1) {
            int middle = low + (high - low) / 2;
            if (threshold(middle).compareTo(n) <= 0) low = middle;
            else high = middle;
        }
        return low;
    }

    private static BigDecimal probability(long n, long k) {
        BigDecimal product = BigDecimal.ONE;
        for (long i = 1; i <= k; ++i) {
            BigDecimal ratio = BigDecimal.valueOf(2 * n - 4 * i - 1).divide(BigDecimal.valueOf(2 * n - 4 * i + 1), MC);
            product = product.multiply(ratio, MC);
        }
        BigDecimal numerator = product.multiply(BigDecimal.valueOf(2 * k * k + 1), MC)
            .subtract(BigDecimal.valueOf(2 * (k * k - 1)), MC);
        return numerator.divide(BigDecimal.valueOf(3 * (2 * k + 1)), MC);
    }

    private static BigDecimal reciprocal(BigInteger value) {
        return BigDecimal.ONE.divide(new BigDecimal(value), MC);
    }

    private static BigDecimal blockSum(long k, BigInteger n, BigDecimal[] central) {
        BigInteger first = threshold(k);
        BigInteger last = n.min(threshold(k + 1).subtract(BigInteger.TWO));
        if (last.compareTo(first) < 0) return BigDecimal.ZERO;

        BigDecimal harmonic = digamma(new BigDecimal(last).subtract(HALF))
            .subtract(digamma(new BigDecimal(first).subtract(ONE_AND_HALF)), MC).divide(TWO, MC);
        BigInteger twiceFirst = first.shiftLeft(1);
        BigInteger twiceLast = last.shiftLeft(1);
        BigDecimal weighted = BigDecimal.ZERO;
        for (long i = 1; i <= k; ++i) {
            BigDecimal residue = central[(int) (i - 1)].multiply(central[(int) (k - i)], MC)
                .multiply(BigDecimal.valueOf(2 * (2 * k - 2 * i + 1)), MC);
            weighted = weighted.add(residue.multiply(harmonic, MC), MC);
            if (i < k) {
                harmonic = harmonic.add(reciprocal(twiceFirst.subtract(BigInteger.valueOf(4 * i + 3))), MC)
                    .add(reciprocal(twiceFirst.subtract(BigInteger.valueOf(4 * i + 1))), MC)
                    .subtract(reciprocal(twiceLast.subtract(BigInteger.valueOf(4 * i + 1))), MC)
                    .subtract(reciprocal(twiceLast.subtract(BigInteger.valueOf(4 * i - 1))), MC);
            }
        }
        BigDecimal correction = weighted.multiply(BigDecimal.valueOf(2 * k * k + 1), MC).divide(THREE, MC);
        return new BigDecimal(last.subtract(first).add(BigInteger.ONE)).subtract(correction, MC)
            .divide(BigDecimal.valueOf(2 * k + 1), MC);
    }

    private static BigDecimal[] centralMoments(long k) {
        BigDecimal inverse = BigDecimal.ONE.divide(BigDecimal.valueOf(k - 1), MC);
        BigDecimal[] moments = new BigDecimal[SERIES_ORDER + 1];
        Arrays.fill(moments, BigDecimal.ZERO);
        moments[0] = BigDecimal.ONE;
        for (int j = 2; j <= SERIES_ORDER; ++j) {
            for (int degree = j - 1; degree >= 0; --degree) {
                int coefficient = degree < MOMENT_COEFFICIENTS[j].length ? MOMENT_COEFFICIENTS[j][degree] : 0;
                moments[j] = moments[j].multiply(inverse, MC).add(BigDecimal.valueOf(coefficient), MC);
            }
        }
        return moments;
    }

    private static Result seriesBlockSum(long k, BigInteger n) {
        BigInteger first = threshold(k);
        BigInteger last = n.min(threshold(k + 1).subtract(BigInteger.TWO));
        if (last.compareTo(first) < 0) return new Result(BigDecimal.ZERO, BigDecimal.ZERO);
        BigDecimal count = new BigDecimal(last.subtract(first).add(BigInteger.ONE));
        BigDecimal h = BigDecimal.valueOf(k - 1);
        BigDecimal a = new BigDecimal(first.shiftLeft(1).subtract(BigInteger.valueOf(k + 2)));
        BigDecimal b = new BigDecimal(last.add(BigInteger.ONE).shiftLeft(1).subtract(BigInteger.valueOf(k + 2)));
        BigDecimal ia = BigDecimal.ONE.divide(a, MC);
        BigDecimal ib = BigDecimal.ONE.divide(b, MC);
        BigDecimal ua = h.multiply(ia, MC);
        BigDecimal ub = h.multiply(ib, MC);
        BigDecimal ia2 = ia.multiply(ia, MC);
        BigDecimal ib2 = ib.multiply(ib, MC);
        BigDecimal ia4 = ia2.multiply(ia2, MC);
        BigDecimal ib4 = ib2.multiply(ib2, MC);
        BigDecimal ia6 = ia4.multiply(ia2, MC);
        BigDecimal ib6 = ib4.multiply(ib2, MC);
        BigDecimal[] moments = centralMoments(k);
        BigDecimal pa = BigDecimal.ONE;
        BigDecimal pb = BigDecimal.ONE;
        BigDecimal weighted = BigDecimal.ZERO;
        BigDecimal bound = BigDecimal.ZERO;
        for (int j = 0; j <= SERIES_ORDER; ++j) {
            BigDecimal integral = j == 0 ? ln(b.divide(a, WIDE)).divide(TWO, MC)
                : pa.subtract(pb, MC).divide(BigDecimal.valueOf(2L * j), MC);
            BigDecimal second = BigDecimal.valueOf(j + 1).divide(BigDecimal.valueOf(6), MC);
            BigDecimal fourth = BigDecimal.valueOf((long) (j + 1) * (j + 2) * (j + 3)).divide(BigDecimal.valueOf(90), MC);
            BigDecimal sixth = BigDecimal.valueOf((long) (j + 1) * (j + 2) * (j + 3) * (j + 4) * (j + 5))
                .divide(BigDecimal.valueOf(945), MC);
            BigDecimal remainder = sixth.multiply(pa.multiply(ia6, MC).subtract(pb.multiply(ib6, MC), MC), MC);
            BigDecimal harmonic = integral
                .add(pa.multiply(ia, MC).subtract(pb.multiply(ib, MC), MC).divide(TWO, MC), MC)
                .add(second.multiply(pa.multiply(ia2, MC).subtract(pb.multiply(ib2, MC), MC), MC), MC)
                .subtract(fourth.multiply(pa.multiply(ia4, MC).subtract(pb.multiply(ib4, MC), MC), MC), MC)
                .add(remainder, MC);
            weighted = weighted.add(moments[j].multiply(harmonic, MC), MC);
            bound = bound.add(moments[j].abs().multiply(remainder, MC), MC);
            pa = pa.multiply(ua, MC);
            pb = pb.multiply(ub, MC);
        }
        BigDecimal rho = THREE.multiply(h).divide(a, MC);
        BigDecimal tail = count.divide(a, MC).multiply(rho.pow(SERIES_ORDER + 1, MC), MC)
            .divide(BigDecimal.ONE.subtract(rho, MC), MC);
        bound = bound.add(tail, MC);
        BigDecimal coefficient = BigDecimal.valueOf(2 * k * k + 1).multiply(BigDecimal.valueOf(2 * k))
            .divide(BigDecimal.valueOf(3 * (2 * k + 1)), MC);
        BigDecimal value = count.divide(BigDecimal.valueOf(2 * k + 1), MC)
            .subtract(coefficient.multiply(weighted, MC), MC);
        return new Result(value, coefficient.multiply(bound, MC));
    }

    private static BigDecimal[] centralRatios(int count) {
        BigDecimal[] central = new BigDecimal[count];
        central[0] = BigDecimal.ONE;
        for (int j = 1; j < count; ++j) {
            central[j] = central[j - 1].multiply(BigDecimal.valueOf(2L * j - 1), MC).divide(BigDecimal.valueOf(2L * j), MC);
        }
        return central;
    }

    private static Result evaluate(BigInteger n, int threadCount) throws InterruptedException, ExecutionException {
        require(n.compareTo(MAX_LIMIT) <= 0, "supported limit through 10^20");
        if (n.compareTo(BigInteger.valueOf(3)) < 0) return new Result(BigDecimal.ZERO, BigDecimal.ZERO);
        int blocks = blockCount(n);
        int threads = Math.min(threadCount, blocks);
        require(threads > 0, "positive thread count");
        int chunkSize = Math.min(256, Math.max(1, blocks / (4 * threads)));
        BigDecimal[] central = centralRatios(Math.min(blocks, SERIES_START));
        Result[] sums = new Result[(blocks + chunkSize - 1) / chunkSize];
        ExecutorService pool = Executors.newFixedThreadPool(threads);
        try {
            List<Future<?>> workers = new ArrayList<>();
            for (int t = 0; t < threads; ++t) {
                int index = t;
                workers.add(pool.submit(() -> {
                    for (int chunk = index; chunk < sums.length; chunk += threads) {
                        int first = chunk * chunkSize + 1;
                        int last = Math.min(blocks, first + chunkSize - 1);
                        BigDecimal value = BigDecimal.ZERO;
                        BigDecimal errorBound = BigDecimal.ZERO;
                        for (int k = first; k <= last; ++k) {
                            Result current = k < SERIES_START
                                ? new Result(blockSum(k, n, central), BigDecimal.ZERO)
                                : seriesBlockSum(k, n);
                            value = value.add(current.value, MC);
                            errorBound = errorBound.add(current.errorBound, MC);
                        }
                        sums[chunk] = new Result(value, errorBound);
                    }
                }));
            }
            for (Future<?> worker : workers) worker.get();
        } finally {
            pool.shutdown();
        }
        BigDecimal value = BigDecimal.ZERO;
        BigDecimal errorBound = BigDecimal.ZERO;
        for (Result sum : sums) {
            value = value.add(sum.value, MC);
            errorBound = errorBound.add(sum.errorBound, MC);
        }
        require(errorBound.compareTo(new BigDecimal("1e-24")) < 0, "series truncation bound");
        return new Result(value, errorBound);
    }

    private static BigDecimal solve(BigInteger n, int threadCount) throws InterruptedException, ExecutionException {
        return evaluate(n, threadCount).value;
    }

    private static BigDecimal solve(long n, int threadCount) throws InterruptedException, ExecutionException {
        return solve(BigInteger.valueOf(n), threadCount);
    }

    private static BigDecimal cubeRoot(BigDecimal x) {
        BigDecimal y = BigDecimal.valueOf(Math.cbrt(x.doubleValue()));
        for (int iteration = 0; iteration < 6; ++iteration) {
            y = y.multiply(TWO).add(x.divide(y.multiply(y, WIDE), WIDE)).divide(THREE, WIDE);
        }
        return y.round(MC);
    }

    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[] denseEquilibrium(int first, int last) {
        int count = last - first + 1;
        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(first + row, 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 BigInteger thresholdNumerator(long k, BigInteger n) {
        BigInteger numerator = BigInteger.ONE;
        BigInteger denominator = BigInteger.ONE;
        for (long i = 1; i <= k; ++i) {
            numerator = numerator.multiply(n.shiftLeft(1).subtract(BigInteger.valueOf(4 * i + 1)));
            denominator = denominator.multiply(n.shiftLeft(1).subtract(BigInteger.valueOf(4 * i - 1)));
        }
        return numerator.multiply(BigInteger.valueOf(2 * k * k + 1))
            .subtract(denominator.multiply(BigInteger.valueOf(2 * k * k - 2)));
    }

    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 {
        require(parseLimit("100000000000000000000").equals(MAX_LIMIT), "parse 10^20");
        require(parseLimit("18446744073709551616").equals(BigInteger.ONE.shiftLeft(64)), "parse above 64 bits");
        for (String text : new String[] {"", "-1", "+1", "1e20", "1 ", "100000000000000000001",
                "99999999999999999999999999999999999999999"}) {
            boolean rejected = false;
            try {
                parseLimit(text);
            } catch (IllegalStateException error) {
                rejected = true;
            }
            require(rejected, "reject invalid or unsupported limit");
        }
        for (BigInteger n : new BigInteger[] {BigInteger.ZERO, BigInteger.TWO, BigInteger.valueOf(3), BigInteger.ONE.shiftLeft(64), MAX_LIMIT}) {
            int count = blockCount(n);
            require((count == 0 || threshold(count).compareTo(n) <= 0) && threshold(count + 1).compareTo(n) > 0,
                "block search boundaries");
        }
        long largestBlock = blockCount(MAX_LIMIT);
        for (long block : new long[] {largestBlock, largestBlock + 1}) {
            BigInteger exact = BigInteger.valueOf(4).multiply(BigInteger.valueOf(block)).multiply(BigInteger.valueOf(block))
                .multiply(BigInteger.valueOf(block)).add(BigInteger.valueOf(5 * block)).divide(BigInteger.valueOf(6)).add(BigInteger.TWO);
            require(threshold(block).equals(exact), "wide cubic threshold versus unbounded integer");
        }
        int maxBlock = 1;
        while (threshold(maxBlock).compareTo(BigInteger.valueOf(TARGET)) <= 0) ++maxBlock;
        for (int k = 1; k <= maxBlock; ++k) {
            BigInteger first = threshold(k);
            require(thresholdNumerator(k, first.subtract(BigInteger.ONE)).signum() < 0, "exact lower threshold sign");
            require(thresholdNumerator(k, first).signum() > 0, "exact upper threshold sign");
        }
        BigDecimal sum = BigDecimal.ZERO;
        int k = 1;
        for (int n = 3; n <= 200; ++n) {
            if (threshold(k + 1).compareTo(BigInteger.valueOf(n)) <= 0) ++k;
            int last = threshold(k + 1).equals(BigInteger.valueOf(n + 1)) ? n - 1 : n;
            int first = last - 2 * k;
            BigDecimal[] p = denseEquilibrium(first, last);
            BigDecimal mass = BigDecimal.ZERO;
            for (BigDecimal value : p) {
                require(value.signum() > 0, "positive active probabilities");
                mass = mass.add(value, MC);
            }
            require(close(mass, BigDecimal.ONE, EPS), "normalized dense equilibrium");
            for (int i = 1; i <= n; ++i) {
                BigDecimal value = BigDecimal.ZERO;
                for (int j = first; j <= last; ++j) value = value.add(p[j - first].multiply(BigDecimal.valueOf(payoff(i, j)), MC), MC);
                require(value.compareTo(EPS) <= 0, "all pure counter-strategies");
            }
            BigDecimal finalProbability = last == n ? p[p.length - 1] : BigDecimal.ZERO;
            require(close(finalProbability, last == n ? probability(n, k) : BigDecimal.ZERO, EPS),
                "closed probability versus dense payoff matrix");
            sum = sum.add(finalProbability, MC);
            if (n <= 10 || n == 21 || n == 47 || n == 88 || n == 100 || n == 200) {
                require(close(sum, solve(n, 1), EPS), "harmonic block sum versus dense equilibria");
            }
        }
        BigDecimal roundedTolerance = new BigDecimal("5e-11");
        require(close(probability(3, 1), BigDecimal.ONE.divide(BigDecimal.valueOf(9), MC), EPS), "P(3) = 1/9");
        require(close(probability(4, 1), BigDecimal.ONE.divide(BigDecimal.valueOf(5), MC), EPS), "P(4) = 1/5");
        require(close(probability(10, 2), new BigDecimal("0.0479638009"), roundedTolerance), "P(10)");
        require(close(solve(10, 1), new BigDecimal("1.1546112276"), roundedTolerance), "S(10)");
        require(close(solve(100, 1), new BigDecimal("4.8779925686"), roundedTolerance), "S(100)");
        require(close(solve(TARGET, 1), solve(TARGET, threadCount), EPS), "thread consistency");
        for (int block = 2; block <= 12; ++block) {
            BigDecimal[] central = centralRatios(block);
            BigDecimal[] direct = new BigDecimal[SERIES_ORDER + 1];
            Arrays.fill(direct, BigDecimal.ZERO);
            for (int j = 0; j < block; ++j) {
                BigDecimal mass = central[j].multiply(BigDecimal.valueOf(2L * block - 2L * j - 1), MC)
                    .multiply(central[block - j - 1], MC).divide(BigDecimal.valueOf(block), MC);
                BigDecimal z = BigDecimal.valueOf(4L * j).divide(BigDecimal.valueOf(block - 1), MC).subtract(BigDecimal.ONE, MC);
                BigDecimal power = BigDecimal.ONE;
                for (int order = 0; order <= SERIES_ORDER; ++order) {
                    direct[order] = direct[order].add(mass.multiply(power, MC), MC);
                    power = power.multiply(z, MC);
                }
            }
            BigDecimal[] moments = centralMoments(block);
            for (int order = 0; order <= SERIES_ORDER; ++order) {
                require(close(direct[order], moments[order], EPS), "central moments versus partial-fraction residues");
            }
        }
        long referenceBlock = blockCount(BigInteger.TEN.pow(16));
        for (long block : new long[] {64, 65, 1000, referenceBlock}) {
            BigDecimal[] central = centralRatios((int) block);
            BigInteger first = threshold(block);
            BigInteger last = threshold(block + 1).subtract(BigInteger.TWO);
            for (BigInteger end : new BigInteger[] {first, first.add(last.subtract(first).shiftRight(1)), last}) {
                Result approximation = seriesBlockSum(block, end);
                require(close(approximation.value, blockSum(block, end, central),
                    approximation.errorBound.add(new BigDecimal("1e-35"), MC)), "bounded series versus exact harmonic block");
            }
        }
        require(close(solve(1_000_000, 1), solve(1_000_000, threadCount), EPS), "series thread consistency");
        require(close(solve(100_000_000, 1), solve(100_000_000, threadCount), EPS), "chunked thread consistency");
        System.out.println("All checks passed.");
    }

    private static void printTable(int threadCount) throws InterruptedException, ExecutionException {
        System.out.println("Hardware threads: " + threadCount + "; three runs per input.");
        System.out.println("| N | S(N) | S(N)/N^(1/3) | Median seconds | Series error bound |");
        System.out.println("|---:|---:|---:|---:|---:|");
        solve(TARGET, threadCount);
        BigInteger n = BigInteger.ONE;
        for (int exponent = 1; exponent <= 20; ++exponent) {
            n = n.multiply(BigInteger.TEN);
            double[] seconds = new double[3];
            Result result = null;
            for (int repeat = 0; repeat < seconds.length; ++repeat) {
                long start = System.nanoTime();
                Result current = evaluate(n, threadCount);
                seconds[repeat] = (System.nanoTime() - start) / 1e9;
                if (result != null) require(close(current.value, result.value, EPS), "table repeat consistency");
                result = current;
            }
            Arrays.sort(seconds);
            BigDecimal ratio = result.value.divide(cubeRoot(new BigDecimal(n)), MC);
            System.out.printf(Locale.ROOT, "| 10^%d | %s | %s | %.6f | %.2e |%n", exponent,
                result.value.setScale(24, RoundingMode.HALF_EVEN).toPlainString(),
                ratio.setScale(23, RoundingMode.HALF_EVEN).toPlainString(), seconds[1], result.errorBound.doubleValue());
        }
    }

    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 == 1 && args[0].equals("--table")) {
                printTable(threadCount);
                return;
            }
            BigInteger n = BigInteger.valueOf(TARGET);
            if (args.length == 2 && args[0].equals("--limit")) {
                n = parseLimit(args[1]);
            } else if (args.length != 0) {
                throw new IllegalArgumentException("Usage: Euler1012 [--self-test | --table | --limit N]");
            }
            BigDecimal answer = solve(n, 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);
        }
    }
}

\(k\)	block \(t_k\le n\)\({}\le t_{k+1}-2\)	\(P(n)=0\) at	\(w_{k,j}\)
\(1\)	\(3\)–\(7\)	\(8\)	\(2\)
\(2\)	\(9\)–\(20\)	\(21\)	\(3\), \(1\)
\(3\)	\(22\)–\(46\)	\(47\)	\(\tfrac{15}4\), \(\tfrac32\), \(\tfrac34\)
\(4\)	\(48\)–\(87\)	\(88\)	\(\tfrac{35}8\), \(\tfrac{15}8\), \(\tfrac98\), \(\tfrac58\)
\(5\)	\(89\)–\(149\)	\(150\)	\(\tfrac{315}{64}\), \(\tfrac{35}{16}\), \(\tfrac{45}{32}\), \(\tfrac{15}{16}\), \(\tfrac{35}{64}\)
\(6\)	\(151\)–\(234\)	\(235\)	\(\dots\)
\(53\)	\(99\,297\)–\(105\,021\)	\(105\,022\)	\(\dots\)