Problem 1015: Partisan Kayles

Project Euler Problem 1015 Solution

A row of blue and red pins is fair exactly when it has as many blue as red pins and the game of knocking over mixed blue-red pairs, which is Dawson's Kayles on its alternating blocks, has nim-sum zero. Sixteen Walsh signs and a square-root generating function count these rows for 1000 pins. Implementations are available in C++, Python and Java.

Detailed mathematical approach

Problem Summary

Blue and red bowling pins stand in a row, and two players, Blue and Red, take turns. On a turn a player knocks over one pin of their own colour, or two neighbouring pins of which at least one has their colour. A knocked-over pin leaves a gap, so the pins on its two sides are no longer neighbours. The player who cannot move loses. A starting row is called fair if the player who moves second wins when both players play perfectly: $\mathtt{BRRB}$ is fair, while $\mathtt{BRBR}$, $\mathtt{BBRR}$ and $\mathtt{BRBB}$ are not ($\mathtt B$ is a blue pin and $\mathtt R$ a red pin). With $f(n)$ the number of fair rows of $n$ pins, we are given $f(6)=8$ and $f(20)=49844$, and we must find the last $9$ digits of $f(1000)$. ^[1] The game is a two-coloured cousin of Kayles, a game invented by Henry Dudeney in which players take turns knocking over one pin or two neighbouring pins from a row of bowling pins. ^[2]

The whole solution in one sentence: a row is fair exactly when it has as many blue pins as red pins and a small hidden game inside it, the game of knocking over mixed pairs, is lost by the player who has to move; both conditions can then be counted for $1000$ pins with a generating function. The sections below build this up slowly, starting from tiny rows, so that every step can be checked by hand.

Mathematical Approach

We first find out who wins any row, without playing the game to the end. Then we count the rows in which the second player wins.

1. Small rows

Perfect play means that each player always picks a move that keeps a win if such a move exists. Every row then has exactly one of four results: Blue wins no matter who starts, Red wins no matter who starts, the first player wins, or the second player wins. Only the last kind is fair. Three tiny examples:

$\mathtt B$: if Blue moves first, Blue knocks the pin over and Red has no move left, so Red loses. If Red moves first, Red has nothing to knock over and loses at once. Blue wins either way.
$\mathtt{BR}$: whoever moves first knocks over both pins, which is allowed because one of them has the mover's colour, and the other player is left with nothing. The first player wins, so the row is not fair.
$\mathtt{BBRR}$: the first player knocks over the middle pair $\mathtt{BR}$, which either player may do. One blue pin and one red pin are left, far apart. Now each player can only knock over their own last pin, and the second player has to do it first, so the second player is the one left without a move. The first player wins, so this row is not fair either.

The fair row $\mathtt{BRRB}$ can be checked by hand as well. Suppose Blue moves first. If Blue knocks over one blue pin, Red knocks over the mixed pair that holds the other blue pin, and Blue has no move left. If Blue knocks over one end pair, $\mathtt{BR}$ or $\mathtt{RB}$, Red knocks over the other end pair, and again Blue has no move left. Blue may not touch the middle pair $\mathtt{RR}$, which has no blue pin. So Blue loses. In the same way Red loses when Red moves first; Red's one extra move, knocking over $\mathtt{RR}$, leaves two blue pins and no red pin, so after Blue's reply Red has nothing to knock over. The second player always wins, so the row is fair.

Checking longer rows by hand quickly becomes hopeless, because every move creates a new position with its own possible moves. So we look for a pattern instead.

2. Pins are spare moves

Think of every blue pin as a sweet that belongs to Blue and every red pin as a sweet that belongs to Red. A player who knocks over a single pin of their own colour eats one of their own sweets. There are exactly three kinds of moves:

knock over one of your own pins: you lose one sweet;
knock over two neighbouring pins of your own colour: you lose two sweets;
knock over a mixed pair, one blue pin next to one red pin: each player loses one sweet.

Let $d$ be the number of blue pins minus the number of red pins, the lead of Blue. The first two kinds of move bring the lead closer to the opponent's side: Blue's own moves lower $d$ by $1$ or $2$, and Red's own moves raise it by $1$ or $2$. A mixed pair does not change $d$ at all; it is a fair trade. In the language of combinatorial game theory an integer counts free moves, the move advantage of one player, ^[6] and $d$ will play exactly that role here.

Mixed pairs have one more special property: both players may knock over any mixed pair, because it always contains a pin of each colour. Moves that are available to both players in the same way form an impartial game, and impartial games are well understood. ^[5]

3. Mixed pairs live in alternating blocks

Cut the row between every two neighbouring pins that have the same colour, and, once pins have been knocked over, also at every gap. The pieces are alternating blocks, in which the colours take turns. For example, $\mathtt{BBRRBR}$ is cut into $\mathtt B$, $\mathtt{BR}$ and $\mathtt{RBR}$, blocks of lengths $1$, $2$ and $3$. A mixed pair is exactly a pair of neighbouring pins inside one block.

If a mixed pair is knocked out of a block of length $L$, two shorter blocks of lengths $i$ and $L-2-i$ remain, separated by the new gap, and each of them still alternates. Colours no longer matter inside a block: the mixed-pair game on one block is a row of $L$ counters from which a move removes two neighbouring counters. This game is known as Dawson's Kayles, the octal game $0.07$. ^[3]

4. A number for every block

For impartial games there is a wonderful shortcut, the Sprague–Grundy theorem: every position gets a number, its Grundy value, and a position with value $0$ is lost by the player who has to move in it. ^[5] The value of a position is the smallest number that is not among the values of the positions reachable in one move; this smallest missing number is called the mex. Writing $\oplus$ for the way values of separate pieces are combined, which Section 5 explains, the value of a block of $L$ pins is

$$g(L)=\operatorname{mex}\{\,g(i)\oplus g(L-2-i)\,\},$$

where $i$ runs over $0,1,\dots,L-2$, and $g(0)=g(1)=0$ because no mixed pair fits into a block of $0$ or $1$ pin, and the smallest number missing from an empty list is $0$. A block of $2$ pins has a single move, which leaves an empty block of value $0$, so $g(2)=1$. A block of $3$ pins always leaves one single pin, of value $0$, so $g(3)=1$ as well. Next, from a block of $4$ pins one can knock over pins one and two, leaving a block of $2$ with value $1$; pins two and three, leaving two single pins with value $0\oplus0=0$; or pins three and four, again with value $1$. The reachable values are $0$ and $1$, the smallest missing number is $2$, and so $g(4)=2$. Continuing in this way gives the values for $L=0,1,\dots,20$:

$0, 0, 1, 1, 2, 0, 3, 1, 1, 0, 3,$
$3, 2, 2, 4, 0, 5, 2, 2, 3, 3.$

Up to $L=1000$ the largest value is $9$, reached for the first time at $L=86$. The last irregular value is at $L=52$; from $L=53$ on the values repeat with period $34$; the On-Line Encyclopedia of Integer Sequences lists the same values, shifted by one place, under Dawson's Chess, together with this period. ^[4]

5. Adding blocks with the nim-sum

Write each value in binary, with only the digits $0$ and $1$, where each place is worth twice the place to its right: $2=10_2$, $3=11_2$, $5=101_2$ and $6=110_2$. Think of every binary digit as a light switch. Combining two values flips the switches of the first value wherever the second value has a $1$; this is the nim-sum $\oplus$, also called exclusive or. For example $1\oplus1=0$; $2\oplus3=1$, because flipping both switches of $10_2$ gives $01_2$; and $5\oplus3=6$, because flipping the last two switches of $101_2$ gives $110_2$. A game made of several separate pieces, where each move happens in one piece, has as its value the nim-sum of the values of the pieces. ^[5] The mixed-pair game of a whole row, with blocks of lengths $L_1,\dots,L_k$, therefore has the value

$$x=g(L_1)\oplus g(L_2)\oplus\dots\oplus g(L_k).$$

Two facts follow from the definition of the mex. If $x\ne0$, some mixed-pair move makes the value $0$. If $x=0$, every mixed-pair move makes it different from $0$.

6. The winner rule

Now the whole game can be decided with two numbers, the lead $d$ and the value $x$:

if $d\gt0$, Blue wins, no matter who starts;
if $d\lt0$, Red wins, no matter who starts;
if $d=0$ and $x\ne0$, the first player wins;
if $d=0$ and $x=0$, the second player wins, so the row is fair.

Put simply: whoever has more sweets wins, because the other player runs out of moves first. When both have the same number of sweets, nobody wants to eat one of their own, since that hands the lead to the other player, so only the mixed-pair game matters.

The precise proof goes by induction on the number of pins, assuming the rule for all smaller positions; the rule is about any position, also one with gaps. Suppose $d\gt0$. If Blue moves first and the row has a mixed pair, Blue knocks it over: $d$ stays positive and Red is to move, so Blue wins by the induction hypothesis. If there is no mixed pair, Blue knocks over one blue pin, which exists because $d\gt0$. Then $d-1\ge0$, and there is still no mixed pair, because knocking pins over never creates new neighbours, so $x=0$; with Red to move, both $d-1\gt0$ and the case $d-1=0$, $x=0$ are wins for Blue. If Red moves first, every move of Red leaves $d$ positive or raises it, and Blue to move wins again; a Red player without moves simply loses. The case $d\lt0$ is the mirror image. Finally suppose $d=0$. If $x\ne0$, the player to move knocks over a mixed pair that makes $x=0$ and leaves $d=0$, and the opponent then loses. If $x=0$, a mixed-pair move leaves $d=0$ and $x\ne0$, so the opponent wins; any other move hands the opponent a positive lead, which also wins for the opponent. So the player to move loses exactly when $x=0$.

The examples of the statement agree. $\mathtt{BRRB}$ has $d=0$ and the blocks $\mathtt{BR}$, $\mathtt{RB}$, so $x=1\oplus1=0$: fair. $\mathtt{BRBR}$ is a single block of length $4$ with $x=2$, and $\mathtt{BBRR}$ has the blocks $\mathtt B$, $\mathtt{BR}$, $\mathtt R$ with $x=0\oplus1\oplus0=1$: in both the first player wins. $\mathtt{BRBB}$ has $d=2$, so Blue always wins. In game-theory language, every row has the same winner as the number $d$ plus the nimber $*x$. ^[6]

7. What has to be counted

A row of $n$ pins is fair exactly when it has $n/2$ blue pins and $x=0$. Rows of odd length are never fair. For $n=6$ there are eight fair rows, which is the given $f(6)=8$: $\mathtt{BBRRBR}$, $\mathtt{BBRRRB}$, $\mathtt{BRBBRR}$, $\mathtt{BRRRBB}$, $\mathtt{RBBBRR}$, $\mathtt{RBRRBB}$, $\mathtt{RRBBBR}$ and $\mathtt{RRBBRB}$. For instance $\mathtt{BRBBRR}$ has the blocks $\mathtt{BRB}$, $\mathtt{BR}$, $\mathtt R$ with $x=1\oplus1\oplus0=0$.

Write $n=2m$, and call a row balanced when it has as many blue pins as red pins, that is, $d=0$. Without the condition $x=0$ the count would be easy: we only have to choose which $m$ of the $2m$ places are blue, so there are $\binom{2m}{m}$ balanced rows of length $2m$, the central binomial coefficient. ^[7] For $n=6$ there are $\binom63=20$ balanced rows, and only $8$ of them are fair. The difficulty is to count only the balanced rows with $x=0$.

8. The sign trick

Every block value up to $1000$ pins is below $16$, so it has at most four binary digits. For two numbers $h$ and $v$ between $0$ and $15$, let $h\cdot v$ be the number of binary places in which both have the digit $1$; here the dot is not ordinary multiplication. For example, $3=0011_2$ and $5=0101_2$ both have a $1$ only in the last place, so for $h=3$ and $v=5$ we get $h\cdot v=1$. Then

$$\frac1{16}\sum_{h=0}^{15}(-1)^{h\cdot v}$$

equals $1$ when $v=0$ and $0$ otherwise. The reason is simple. If $v=0$, every term is $+1$. If $v$ has a digit $1$ in some place, flipping that digit of $h$ changes the sign of the term, so the sixteen terms split into pairs of $+1$ and $-1$ that cancel, like one step forward and one step back. These signs are the characters behind the Walsh–Hadamard transform. ^[8]

The sign of a nim-sum is the product of the signs of its parts. Only the places where $h$ has a $1$ matter. If exactly one part has a $1$ in such a place, the nim-sum has a $1$ there too, and both sides get one minus sign. If both parts have a $1$ there, the product gets two minus signs, which cancel, and the nim-sum has a $0$ there and gets no minus sign. So for a row with blocks of lengths $L_1,\dots,L_k$,

$$(-1)^{h\cdot x}=w_h(L_1)\,w_h(L_2)\cdots w_h(L_k),$$

with $w_h(L)=(-1)^{h\cdot g(L)}$. Let $c_h(m)$ be the sum of these products over all balanced rows of length $2m$. Averaging over $h$ keeps exactly the rows with $x=0$:

$$f(2m)=\frac1{16}\sum_{h=0}^{15}c_h(m).$$

For each of the $16$ signs $h$ we now need $c_h(m)$, a sum of products over blocks, and such sums are what generating functions are good at. A generating function keeps a whole list of numbers in one expression: the number in front of $t^n$ is the entry for $n$ pins. When the list goes on forever, as in $1+t+t^2+t^3+\dots$, the expression is called a power series.

9. Building a row block by block

A row is the colour of its first pin followed by a list of block lengths; two neighbouring blocks always meet with two pins of the same colour. A block of odd length starts and ends with the same colour, so the next block starts with that colour again, and the block changes the lead by $+1$ if that colour is blue and by $-1$ if it is red. A block of even length ends with the other colour, so the next block starts with the other colour, and it does not change the lead.

Let $t$ count pins and $y$ count the lead, so that a row of $n$ pins and lead $d$ contributes $t^n y^d$ times its product of signs. For example, $\mathtt{BRB}$ has $3$ pins and lead $1$, so it contributes $t^3y$ times its sign, and $\mathtt{RBR}$ contributes $t^3y^{-1}$ times its sign, where $y^{-1}=1/y$. Collect the signed blocks of each kind in two power series:

$$E=\sum_{j\ge1}w_h(2j)\,t^{2j},$$

$$O=\sum_{j\ge0}w_h(2j+1)\,t^{2j+1}.$$

Remember the colour with which the next block starts, blue or red. One block then acts like the matrix

$$P=\begin{pmatrix}yO&E\\E&y^{-1}O\end{pmatrix},$$

a small table whose row is the colour with which a block starts and whose column is the colour with which the next block starts: the first row describes a block that starts with blue and the second row a block that starts with red. Multiplying such tables chains blocks together, so $P^2$ describes two blocks one after the other, $P^3$ three blocks, and so on. Adding up $I$, $P$, $P^2$, $P^3$ and so on, where the identity matrix $I$ stands for no block at all, over both colours of the first pin, gives the signed rows of every length. Just as $1+q+q^2+\dots=1/(1-q)$, this sum is $(I-P)^{-1}$. The empty row appears once for each starting colour, and one copy is removed:

$$F=\begin{pmatrix}1&1\end{pmatrix}(I-P)^{-1}\begin{pmatrix}1\\1\end{pmatrix}-1.$$

Inverting the $2\times2$ matrix and simplifying gives

$$F=\frac{(1+E)^2-O^2}{1+O^2-E^2-(y+y^{-1})O}.$$

10. Keeping the colours balanced

The balanced rows are the terms with $y^0$. Picture the lead as a walk that steps up for each blue odd block and down for each red one: we only keep the walks that come back to where they started. For power series $a$ with constant term $1$ and $b$ with constant term $0$, where the constant term is the part without any $t$, the rule $1/(1-q)=1+q+q^2+\dots$ with $q=b(y+y^{-1})/a$ gives

$$\frac1{a-b(y+y^{-1})}=\sum_{k\ge0}\frac{b^k(y+y^{-1})^k}{a^{k+1}}.$$

The power $(y+y^{-1})^k$ has a term with $y^0$ only when $k=2j$ is even, and then its coefficient is $\binom{2j}{j}$, the number of ways to take $j$ steps up and $j$ steps down. Since $\sum_j\binom{2j}{j}z^j=1/\sqrt{1-4z}$, ^[7] the part with $y^0$ is

$$\frac1{a}\cdot\frac1{\sqrt{1-4b^2/a^2}}=\frac1{\sqrt{a^2-4b^2}}.$$

Here $a=1+O^2-E^2$ and $b=O$, and $a^2-4b^2$ factors as $\bigl((1-E)^2-O^2\bigr)$$\bigl((1+E)^2-O^2\bigr)$. Cancelling the common factor with the numerator of $F$ gives the series of the balanced rows, which all have even length:

$$\sum_{m\ge0}c_h(m)\,t^{2m}=\sqrt{\frac{(1+E)^2-O^2}{(1-E)^2-O^2}}.$$

This is the formula the program evaluates, once for each of the $16$ values of $h$.

11. A check with all signs equal to plus one

For $h=0$ every sign is $+1$, and the formula must simply count all balanced rows. Indeed, with $u=t^2$ we get $E=u/(1-u)$ and $O^2=u/(1-u)^2$, so $(1+E)^2-O^2=1/(1-u)$ and $(1-E)^2-O^2$${}=(1-4u)/(1-u)$. The fraction under the root is $1/(1-4u)$, and its square root $1/\sqrt{1-4u}$ is the generating function of $\binom{2m}{m}$, exactly the number of balanced rows found in Section 7. ^[7]

12. Only the last nine digits

Only the last nine digits of $f(1000)$ are needed, that is, $f(1000)$ modulo $10^9$: the remainder left after dividing by $10^9$. At the end, however, we divide by $16$. So every $c_h(m)$ is computed modulo $M=16\cdot10^9$. The true sum of the sixteen values is $16f(1000)$, and $16$ divides $M$, so the sum modulo $M$ is divisible by $16$ as well, and dividing it by $16$ gives $f(1000)$ modulo $10^9$.

With $u=t^2$, write $E=\sum_{j\ge1}e_ju^j$ and $O=t\sum_{j\ge0}o_ju^j$, where $e_j=w_h(2j)$ and $o_j=w_h(2j+1)$, so that $O^2=u\,(\sum_{j\ge0}o_ju^j)^2$. There is no block of length $0$, so $E$ has no constant term. The numerator $(1+E)^2-O^2$ and the denominator $(1-E)^2-O^2$ are power series in $u$ with small integer coefficients, and the denominator starts with $1$. Their quotient $s=\sum_m s_mu^m$ is found term by term, modulo $2M$, like a long division: because the denominator starts with $1$, each $s_m$ follows from the numerator and the earlier $s_0,\dots,s_{m-1}$. Its square root $r$, with $r_0=1$, satisfies $r^2=s$. The coefficient of $u^m$ in $r^2$ is $r_0r_m+r_1r_{m-1}+\dots+r_mr_0$; its first and last terms together give $2r_m$, and setting it equal to $s_m$ gives

$$2r_m=s_m-\sum_{i=1}^{m-1}r_ir_{m-i}.$$

Halving needs care when we only know the last digits: the program computes the right-hand side modulo $2M$ and then halves it, which gives $r_m$ modulo $M$. This is enough. The terms $r_ir_{m-i}$ with $i\ne m-i$ occur in equal pairs, and adding $M$ to $r_i$ changes $2r_ir_{m-i}$ by $2Mr_{m-i}$, a multiple of $2M$; so $2r_ir_{m-i}$ modulo $2M$ depends only on $r_i$ and $r_{m-i}$ modulo $M$. When $m$ is even there is also a middle term $r_{m/2}^2$; adding $M$ to $r_{m/2}$ changes it by $2Mr_{m/2}+M^2$, and $M^2$ is a multiple of $2M$ because $M$ is even. So the middle term modulo $2M$ also depends only on $r_{m/2}$ modulo $M$. The coefficients $r_m$ are whole numbers, because they are signed counts of rows, so the right-hand side is always even.

How the Code Works

The C++ program ^[9] starts with grundy_values, which computes $g(L)$ for $L\le1000$ from the mex rule of Section 4, marking the reachable values of each $L$ in a table. character_sign returns $(-1)^{h\cdot v}$ by counting the binary ones of the common digits. balanced_character handles one value of $h$: it fills the signs $e_j=w_h(2j)$ and $o_j=w_h(2j+1)$, forms the numerator and denominator of Section 12 as integer coefficient lists, divides them modulo $2M$ and takes the square root with the halving rule. A check stops the program if the number to be halved is ever odd.

fair_counts finds that the largest value $9$ needs $16$ signs, sets $M=16\cdot10^9$, and computes the sixteen series. It adds them, checks that every sum is divisible by $16$ and divides. Without arguments the main program prints the last nine digits of $f(1000)$, with leading zeros.

With --self-test, run_tests checks $f(0)=1$, $f(6)=8$ and $f(20)=49844$. brute_outcomes then solves the original game directly, with all $3^{12}=531441$ positions of up to $12$ places, each place empty, blue or red, and predicted_outcome confirms that the winner rule of Section 6 gives the same result for every one of them. Counting the fair starting rows among these positions reproduces $f(n)$ for every $n\le12$, and for $h=0$ the series of Section 10 agrees with $\binom{2m}{m}$ modulo $10^9$ for every $2m\le1000$. The Python and Java versions follow the C++ program function by function; Python uses exact integers, and Java uses overflow-checked long arithmetic and BigInteger where the C++ program uses 128-bit integers.

Complexity and Verification

The Grundy table needs about $n^2/4$ steps for $n=1000$ pins. Each of the $16$ series has $500$ coefficients, and each coefficient of the quotient and of the square root needs at most $500$ multiplications, so the whole computation takes a few million arithmetic operations and finishes in a fraction of a second. The self-test, which also solves all $531441$ small positions by brute force, takes well under a second.

All three programs print the same nine digits and pass the self-test. Two independent checks, which are not part of the published solutions, confirm the result. A separately written solver of the original game checks the winner rule on every position of up to $20$ places, $3^{20}=3486784401$ positions, and counts the fair rows by brute force for every $n\le20$, including $f(20)=49844$. A direct dynamic programme walks along the row pin by pin, remembering the colour of the last pin, the lead, the nim-sum of the finished blocks and the length of the current block, without any generating function or sign trick; it gives the exact $f(n)$ for every even $n\le200$ and the same last nine digits of $f(n)$ for every even $n\le1000$.

Footnotes and References

The references below supply the background facts used above. The winner rule with the lead $d$ and the value $x$, the generating function of the signed rows and its balanced part, and the halving rule of Section 12 are derived explicitly in this article.

Project Euler 1015 — Partisan Kayles. The official statement: Blue and Red knock over one pin of their own colour or two adjacent pins of which at least one has their colour, a knocked-over pin leaves a gap, and the player who cannot move loses. A starting row is fair if, with perfect play by both players, the second player wins; BRRB is fair, BRBR, BBRR and BRBB are not, f(6) = 8, f(20) = 49844, and the task asks for the last 9 digits of f(1000).
Wikipedia — Kayles. Kayles, invented by Henry Dudeney, is played on a row of bowling pins; on a turn a player knocks out one pin or two adjacent pins, and the player without a legal move loses. In octal notation it is the game 0.77.
Wikipedia — Octal game. In the game 0.07, called Dawson's Kayles, a move removes exactly two tokens from a heap and leaves the rest as zero, one or two heaps. For a row, this is knocking over two neighbouring pins, as in the mixed-pair game of Section 3.
OEIS A002187 — Sprague–Grundy values for Dawson's Chess. The entry notes that Dawson's Kayles, the octal game .07, has the values a(n − 1), and that the sequence has period 34 with the only exceptions at n = 0, 14, 16, 17, 31, 34 and 51. In the numbering of this article, g(L) = a(L − 1) for L ≥ 1.
Wikipedia — Sprague–Grundy theorem. Every impartial game under the normal play convention is equivalent to a nimber; the value of a position is the mex of the values of its options, and the value of a sum of games is the nim-sum of the values of its parts.
Wikipedia — Combinatorial game theory. Numbers represent the number of free moves, or the move advantage of one player, with positive numbers favouring Left; sums of number games behave like the integers. Star, written ∗ = {0 | 0}, is a first-player win, and impartial games correspond to nimbers.
Wikipedia — Central binomial coefficient. The central binomial coefficients C(2n, n) are 1, 2, 6, 20, 70, 252, … and their ordinary generating function is 1/√(1 − 4x).
Wikipedia — Hadamard transform. The entries of the Walsh–Hadamard matrix are (−1)^(i·j), up to a normalising factor, where i·j is the bitwise dot product of the binary representations of i and j; the characters of (ℤ/2ℤ)ⁿ have the same form.
C++ — Euler1015.cpp. The linked immutable C++ revision contains grundy_values, balanced_character, fair_counts and the checks run by --self-test. It is the implementation record for this article.

Problem 1015 source code

C++

#include <algorithm>
#include <cstdint>
#include <cstdlib>
#include <exception>
#include <iomanip>
#include <iostream>
#include <pthread.h>
#include <stdexcept>
#include <string>
#include <thread>
#include <vector>

namespace {

using u64 = std::uint64_t;
using i128 = __int128_t;
using u128 = __uint128_t;
constexpr unsigned LIMIT = 1000;
constexpr u64 MOD = 1'000'000'000;

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

std::vector<unsigned> grundy_values(const unsigned limit) {
    std::vector<unsigned> grundy(limit + 1), seen(limit + 2);
    for (unsigned n = 2; n <= limit; ++n) {
        for (unsigned left = 0; left <= (n - 2) / 2; ++left) {
            const unsigned value = grundy[left] ^ grundy[n - 2 - left];
            if (value < seen.size()) seen[value] = n;
        }
        while (seen[grundy[n]] == n) ++grundy[n];
    }
    return grundy;
}

int character_sign(unsigned value) {
    int sign = 1;
    for (; value; value &= value - 1) sign = -sign;
    return sign;
}

u64 residue(i128 value, const u64 modulus) {
    value %= static_cast<i128>(modulus);
    if (value < 0) value += modulus;
    return static_cast<u64>(value);
}

std::vector<u64> balanced_character(const std::vector<unsigned>& grundy,
                                  const unsigned character, const u64 modulus) {
    const unsigned degree = static_cast<unsigned>(grundy.size() - 1) / 2;
    std::vector<int> odd(degree), even(degree + 1), numerator(degree + 1), denominator(degree + 1);
    for (unsigned n = 0; n < degree; ++n) odd[n] = character_sign(character & grundy[2 * n + 1]);
    for (unsigned n = 1; n <= degree; ++n) even[n] = character_sign(character & grundy[2 * n]);
    numerator[0] = denominator[0] = 1;
    for (unsigned n = 1; n <= degree; ++n) {
        int common = 0;
        for (unsigned i = 1; i < n; ++i) common += even[i] * even[n - i];
        for (unsigned i = 0; i < n; ++i) common -= odd[i] * odd[n - 1 - i];
        numerator[n] = common + 2 * even[n];
        denominator[n] = common - 2 * even[n];
    }

    // The balanced-row series is sqrt(((1+E)^2-A^2)/((1-E)^2-A^2)).
    const u64 double_modulus = 2 * modulus;
    std::vector<u64> square(degree + 1), root(degree + 1);
    square[0] = root[0] = 1;
    for (unsigned n = 1; n <= degree; ++n) {
        i128 coefficient = numerator[n];
        for (unsigned i = 1; i <= n; ++i) {
            coefficient -= static_cast<i128>(denominator[i]) * square[n - i];
        }
        square[n] = residue(coefficient, double_modulus);
        u128 products = 0;
        for (unsigned i = 1; 2 * i < n; ++i) products += u128{2} * root[i] * root[n - i];
        if (n % 2 == 0) products += static_cast<u128>(root[n / 2]) * root[n / 2];
        // For even M, this symmetric sum modulo 2M needs earlier roots only modulo M.
        const u64 twice_root = (square[n] + double_modulus
            - static_cast<u64>(products % double_modulus)) % double_modulus;
        require(twice_root % 2 == 0, "integral square-root coefficient");
        root[n] = twice_root / 2;
    }
    return root;
}

struct CharacterWork {
    const std::vector<unsigned>* grundy = nullptr;
    std::vector<std::vector<u64>>* coefficients = nullptr;
    u64 modulus = 0;
    unsigned first = 0, stride = 1;
    std::exception_ptr error;
};

void* calculate_characters(void* argument) {
    auto& work = *static_cast<CharacterWork*>(argument);
    try {
        for (unsigned h = work.first; h < work.coefficients->size(); h += work.stride) {
            (*work.coefficients)[h] = balanced_character(*work.grundy, h, work.modulus);
        }
    } catch (...) {
        work.error = std::current_exception();
    }
    return nullptr;
}

std::vector<u64> fair_counts(const unsigned limit, unsigned thread_count = 0) {
    if (limit > LIMIT) throw std::invalid_argument("Supported length is at most 1000");
    const auto grundy = grundy_values(limit);
    unsigned characters = 1;
    while (characters <= *std::max_element(grundy.begin(), grundy.end())) characters *= 2;
    const u64 modulus = MOD * characters;
    if (thread_count == 0) {
        thread_count = limit < 128 ? 1 : std::min(8U, std::max(1U, std::thread::hardware_concurrency()));
    }
    thread_count = std::min(thread_count, characters);
    std::vector<std::vector<u64>> coefficients(characters);
    std::vector<CharacterWork> work(thread_count);
    std::vector<pthread_t> threads(thread_count);
    std::vector<bool> started(thread_count);
    for (unsigned i = 0; i < thread_count; ++i) {
        work[i] = {&grundy, &coefficients, modulus, i, thread_count, nullptr};
    }
    for (unsigned i = 1; i < thread_count; ++i) {
        started[i] = pthread_create(&threads[i], nullptr, calculate_characters, &work[i]) == 0;
        if (!started[i]) calculate_characters(&work[i]);
    }
    calculate_characters(&work[0]);
    for (unsigned i = 1; i < thread_count; ++i) {
        if (started[i] && pthread_join(threads[i], nullptr) != 0) std::terminate();
    }
    for (const auto& job : work) {
        if (job.error) std::rethrow_exception(job.error);
    }
    std::vector<u64> counts(limit / 2 + 1);
    for (unsigned character = 0; character < characters; ++character) {
        for (unsigned n = 0; n < counts.size(); ++n) {
            counts[n] = (counts[n] + coefficients[character][n]) % modulus;
        }
    }
    for (u64& count : counts) {
        require(count % characters == 0, "integral Walsh average");
        count /= characters;
    }
    return counts;
}

std::vector<unsigned char> brute_outcomes(const unsigned length, std::vector<unsigned>& powers) {
    powers.assign(length + 1, 1);
    for (unsigned i = 1; i <= length; ++i) powers[i] = 3 * powers[i - 1];
    std::vector<unsigned char> outcomes(powers[length]);
    for (unsigned state = 1; state < outcomes.size(); ++state) {
        unsigned remaining = state;
        unsigned char result = 0;
        for (unsigned i = 0; i < length && result != 3; ++i) {
            const unsigned pin = remaining % 3;
            remaining /= 3;
            if (pin == 0) continue;
            const unsigned single = state - pin * powers[i];
            if (pin == 1 && !(outcomes[single] & 2)) result |= 1;
            if (pin == 2 && !(outcomes[single] & 1)) result |= 2;
            const unsigned next = remaining % 3;
            if (next == 0) continue;
            const unsigned pair = single - next * powers[i + 1];
            if ((pin == 1 || next == 1) && !(outcomes[pair] & 2)) result |= 1;
            if ((pin == 2 || next == 2) && !(outcomes[pair] & 1)) result |= 2;
        }
        outcomes[state] = result;
    }
    return outcomes;
}

unsigned predicted_outcome(unsigned state, const unsigned length, const std::vector<unsigned>& grundy) {
    int balance = 0;
    unsigned previous = 0, run = 0, nim_sum = 0;
    for (unsigned i = 0; i < length; ++i) {
        const unsigned pin = state % 3;
        state /= 3;
        if (pin == 1) ++balance;
        if (pin == 2) --balance;
        if (pin != 0 && previous != 0 && pin != previous) ++run;
        else {
            nim_sum ^= grundy[run];
            run = pin == 0 ? 0 : 1;
        }
        previous = pin;
    }
    nim_sum ^= grundy[run];
    if (balance > 0) return 1;
    if (balance < 0) return 2;
    return nim_sum == 0 ? 0 : 3;
}

void run_tests() {
    constexpr unsigned SMALL = 12;
    const auto counts = fair_counts(LIMIT);
    require(counts[0] == 1, "empty arrangement");
    require(counts[3] == 8, "given f(6)");
    require(counts[10] == 49'844, "given f(20)");
    require(fair_counts(LIMIT, 1) == counts, "serial and pthread agreement");

    std::vector<unsigned> powers;
    const auto outcomes = brute_outcomes(SMALL, powers);
    const auto grundy = grundy_values(LIMIT);
    for (unsigned state = 0; state < outcomes.size(); ++state) {
        require(outcomes[state] == predicted_outcome(state, SMALL, grundy),
                "original-game minimax agrees with balance and alternating blocks");
    }
    for (unsigned n = 0; n <= SMALL; ++n) {
        u64 expected = 0;
        for (unsigned word = 0; word < (1U << n); ++word) {
            unsigned state = 0;
            for (unsigned i = 0; i < n; ++i) state += (1 + ((word >> i) & 1U)) * powers[i];
            expected += outcomes[state] == 0;
        }
        require(expected == (n % 2 == 0 ? counts[n / 2] : 0), "exhaustive initial arrangements");
        require(fair_counts(n).back() == counts[n / 2], "short-prefix consistency");
    }

    std::vector<u64> binomial(LIMIT + 1);
    binomial[0] = 1;
    const auto unsigned_counts = balanced_character(grundy, 0, MOD);
    for (unsigned n = 1; n <= LIMIT; ++n) {
        for (unsigned k = n; k > 0; --k) binomial[k] = (binomial[k] + binomial[k - 1]) % MOD;
        if (n % 2 == 0) require(unsigned_counts[n / 2] == binomial[n / 2], "central binomial coefficients");
    }
    std::cout << "All checks passed; " << outcomes.size()
              << " original-game positions, given examples, and central binomial coefficients through " << LIMIT << ".\n";
}

}

int main(int argc, char* argv[]) {
    try {
        if (argc == 2 && std::string(argv[1]) == "--self-test") run_tests();
        else {
            if (argc != 1) throw std::invalid_argument("Usage: Euler1015 [--self-test]");
            std::cout << std::setfill('0') << std::setw(9) << fair_counts(LIMIT).back() << '\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 1015 - Partisan Kayles."""

import os
import sys

LIMIT = 1000
MOD = 1_000_000_000


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


def grundy_values(limit):
    grundy = [0] * (limit + 1)
    seen = [0] * (limit + 2)
    for n in range(2, limit + 1):
        for left in range((n - 2) // 2 + 1):
            value = grundy[left] ^ grundy[n - 2 - left]
            if value < len(seen):
                seen[value] = n
        while seen[grundy[n]] == n:
            grundy[n] += 1
    return grundy


def character_sign(value):
    sign = 1
    while value:
        sign = -sign
        value &= value - 1
    return sign


def residue(value, modulus):
    # Python's % already returns a value in [0, modulus).
    return value % modulus


def balanced_character(grundy, character, modulus):
    degree = (len(grundy) - 1) // 2
    odd = [0] * degree
    even = [0] * (degree + 1)
    numerator = [0] * (degree + 1)
    denominator = [0] * (degree + 1)
    for n in range(degree):
        odd[n] = character_sign(character & grundy[2 * n + 1])
    for n in range(1, degree + 1):
        even[n] = character_sign(character & grundy[2 * n])
    numerator[0] = denominator[0] = 1
    for n in range(1, degree + 1):
        common = 0
        for i in range(1, n):
            common += even[i] * even[n - i]
        for i in range(n):
            common -= odd[i] * odd[n - 1 - i]
        numerator[n] = common + 2 * even[n]
        denominator[n] = common - 2 * even[n]

    # The balanced-row series is sqrt(((1+E)^2-A^2)/((1-E)^2-A^2)).
    double_modulus = 2 * modulus
    square = [0] * (degree + 1)
    root = [0] * (degree + 1)
    square[0] = root[0] = 1
    for n in range(1, degree + 1):
        coefficient = numerator[n]
        for i in range(1, n + 1):
            coefficient -= denominator[i] * square[n - i]
        square[n] = residue(coefficient, double_modulus)
        products = 0
        i = 1
        while 2 * i < n:
            products += 2 * root[i] * root[n - i]
            i += 1
        if n % 2 == 0:
            products += root[n // 2] * root[n // 2]
        # For even M, this symmetric sum modulo 2M needs earlier roots only modulo M.
        twice_root = (square[n] + double_modulus - products % double_modulus) % double_modulus
        require(twice_root % 2 == 0, "integral square-root coefficient")
        root[n] = twice_root // 2
    return root


def calculate_characters(grundy, coefficients, modulus, first, stride):
    for h in range(first, len(coefficients), stride):
        coefficients[h] = balanced_character(grundy, h, modulus)


def fair_counts(limit, thread_count=0):
    if limit > LIMIT:
        raise ValueError("Supported length is at most 1000")
    grundy = grundy_values(limit)
    characters = 1
    while characters <= max(grundy):
        characters *= 2
    modulus = MOD * characters
    if thread_count == 0:
        thread_count = 1 if limit < 128 else min(8, max(1, os.cpu_count() or 0))
    thread_count = min(thread_count, characters)
    coefficients = [None] * characters
    # The C++ runs one pthread per stride class; here the same classes run one after another.
    for i in range(thread_count):
        calculate_characters(grundy, coefficients, modulus, i, thread_count)
    counts = [0] * (limit // 2 + 1)
    for character in range(characters):
        for n in range(len(counts)):
            counts[n] = (counts[n] + coefficients[character][n]) % modulus
    for n in range(len(counts)):
        require(counts[n] % characters == 0, "integral Walsh average")
        counts[n] //= characters
    return counts


def brute_outcomes(length, powers):
    powers[:] = [1] * (length + 1)
    for i in range(1, length + 1):
        powers[i] = 3 * powers[i - 1]
    outcomes = bytearray(powers[length])
    for state in range(1, len(outcomes)):
        remaining = state
        result = 0
        for i in range(length):
            if result == 3:
                break
            pin = remaining % 3
            remaining //= 3
            if pin == 0:
                continue
            single = state - pin * powers[i]
            if pin == 1 and not (outcomes[single] & 2):
                result |= 1
            if pin == 2 and not (outcomes[single] & 1):
                result |= 2
            next = remaining % 3
            if next == 0:
                continue
            pair = single - next * powers[i + 1]
            if (pin == 1 or next == 1) and not (outcomes[pair] & 2):
                result |= 1
            if (pin == 2 or next == 2) and not (outcomes[pair] & 1):
                result |= 2
        outcomes[state] = result
    return outcomes


def predicted_outcome(state, length, grundy):
    balance = 0
    previous = run = nim_sum = 0
    for _ in range(length):
        pin = state % 3
        state //= 3
        if pin == 1:
            balance += 1
        if pin == 2:
            balance -= 1
        if pin != 0 and previous != 0 and pin != previous:
            run += 1
        else:
            nim_sum ^= grundy[run]
            run = 0 if pin == 0 else 1
        previous = pin
    nim_sum ^= grundy[run]
    if balance > 0:
        return 1
    if balance < 0:
        return 2
    return 0 if nim_sum == 0 else 3


def run_tests():
    SMALL = 12
    counts = fair_counts(LIMIT)
    require(counts[0] == 1, "empty arrangement")
    require(counts[3] == 8, "given f(6)")
    require(counts[10] == 49_844, "given f(20)")
    require(fair_counts(LIMIT, 1) == counts, "serial and strided agreement")

    powers = []
    outcomes = brute_outcomes(SMALL, powers)
    grundy = grundy_values(LIMIT)
    for state in range(len(outcomes)):
        require(outcomes[state] == predicted_outcome(state, SMALL, grundy),
                "original-game minimax agrees with balance and alternating blocks")
    for n in range(SMALL + 1):
        expected = 0
        for word in range(1 << n):
            state = 0
            for i in range(n):
                state += (1 + ((word >> i) & 1)) * powers[i]
            expected += outcomes[state] == 0
        require(expected == (counts[n // 2] if n % 2 == 0 else 0), "exhaustive initial arrangements")
        require(fair_counts(n)[-1] == counts[n // 2], "short-prefix consistency")

    binomial = [0] * (LIMIT + 1)
    binomial[0] = 1
    unsigned_counts = balanced_character(grundy, 0, MOD)
    for n in range(1, LIMIT + 1):
        for k in range(n, 0, -1):
            binomial[k] = (binomial[k] + binomial[k - 1]) % MOD
        if n % 2 == 0:
            require(unsigned_counts[n // 2] == binomial[n // 2], "central binomial coefficients")
    print(f"All checks passed; {len(outcomes)} original-game positions, given examples, "
          f"and central binomial coefficients through {LIMIT}.")


def main(argv):
    try:
        if len(argv) == 2 and argv[1] == "--self-test":
            run_tests()
        else:
            if len(argv) != 1:
                raise ValueError("Usage: Euler1015.py [--self-test]")
            print(f"{fair_counts(LIMIT)[-1]:09d}")
    except (ValueError, RuntimeError) as error:
        print(error, file=sys.stderr)
        return 1
    return 0


if __name__ == "__main__":
    sys.exit(main(sys.argv))

Java

import java.math.BigInteger;
import java.util.ArrayList;
import java.util.Arrays;
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 Euler1015 {
    private static final int LIMIT = 1000;
    private static final long MOD = 1_000_000_000L;
    private static final BigInteger TWO = BigInteger.valueOf(2);

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

    static int[] grundyValues(int limit) {
        int[] grundy = new int[limit + 1], seen = new int[limit + 2];
        for (int n = 2; n <= limit; ++n) {
            for (int left = 0; left <= (n - 2) / 2; ++left) {
                int value = grundy[left] ^ grundy[n - 2 - left];
                if (value < seen.length) seen[value] = n;
            }
            while (seen[grundy[n]] == n) ++grundy[n];
        }
        return grundy;
    }

    private static int characterSign(int value) {
        int sign = 1;
        for (; value != 0; value &= value - 1) sign = -sign;
        return sign;
    }

    private static long residue(long value, long modulus) {
        value %= modulus;
        if (value < 0) value += modulus;
        return value;
    }

    static long[] balancedCharacter(int[] grundy, int character, long modulus) {
        int degree = (grundy.length - 1) / 2;
        int[] odd = new int[degree], even = new int[degree + 1];
        int[] numerator = new int[degree + 1], denominator = new int[degree + 1];
        for (int n = 0; n < degree; ++n) odd[n] = characterSign(character & grundy[2 * n + 1]);
        for (int n = 1; n <= degree; ++n) even[n] = characterSign(character & grundy[2 * n]);
        numerator[0] = denominator[0] = 1;
        for (int n = 1; n <= degree; ++n) {
            int common = 0;
            for (int i = 1; i < n; ++i) common += even[i] * even[n - i];
            for (int i = 0; i < n; ++i) common -= odd[i] * odd[n - 1 - i];
            numerator[n] = common + 2 * even[n];
            denominator[n] = common - 2 * even[n];
        }

        // The balanced-row series is sqrt(((1+E)^2-A^2)/((1-E)^2-A^2)).
        long doubleModulus = Math.multiplyExact(2L, modulus);
        BigInteger bigDoubleModulus = BigInteger.valueOf(doubleModulus);
        long[] square = new long[degree + 1], root = new long[degree + 1];
        square[0] = root[0] = 1;
        for (int n = 1; n <= degree; ++n) {
            // |denominator| <= 2n + 1 and square < 2M, so this sum is exact in a long;
            // the exact operations throw instead of wrapping if that ever stops being true.
            long coefficient = numerator[n];
            for (int i = 1; i <= n; ++i) {
                coefficient = Math.subtractExact(coefficient, Math.multiplyExact((long) denominator[i], square[n - i]));
            }
            square[n] = residue(coefficient, doubleModulus);
            // root < M, so one product can exceed 2^63: accumulate exactly, as the C++ does in 128 bits.
            BigInteger products = BigInteger.ZERO;
            for (int i = 1; 2 * i < n; ++i) {
                products = products.add(TWO.multiply(BigInteger.valueOf(root[i])).multiply(BigInteger.valueOf(root[n - i])));
            }
            if (n % 2 == 0) products = products.add(BigInteger.valueOf(root[n / 2]).multiply(BigInteger.valueOf(root[n / 2])));
            // For even M, this symmetric sum modulo 2M needs earlier roots only modulo M.
            long twiceRoot = Math.subtractExact(Math.addExact(square[n], doubleModulus),
                products.mod(bigDoubleModulus).longValueExact()) % doubleModulus;
            require(twiceRoot % 2 == 0, "integral square-root coefficient");
            root[n] = twiceRoot / 2;
        }
        return root;
    }

    private static final class CharacterWork {
        final int[] grundy;
        final long[][] coefficients;
        final long modulus;
        final int first, stride;
        RuntimeException error;

        CharacterWork(int[] grundy, long[][] coefficients, long modulus, int first, int stride) {
            this.grundy = grundy;
            this.coefficients = coefficients;
            this.modulus = modulus;
            this.first = first;
            this.stride = stride;
        }
    }

    private static void calculateCharacters(CharacterWork work) {
        try {
            for (int h = work.first; h < work.coefficients.length; h += work.stride) {
                work.coefficients[h] = balancedCharacter(work.grundy, h, work.modulus);
            }
        } catch (RuntimeException error) {
            work.error = error;
        }
    }

    static long[] fairCounts(int limit) {
        return fairCounts(limit, 0);
    }

    static long[] fairCounts(int limit, int threadCount) {
        if (limit > LIMIT) throw new IllegalArgumentException("Supported length is at most 1000");
        int[] grundy = grundyValues(limit);
        int characters = 1;
        while (characters <= Arrays.stream(grundy).max().getAsInt()) characters *= 2;
        long modulus = Math.multiplyExact(MOD, (long) characters);
        if (threadCount == 0) {
            threadCount = limit < 128 ? 1 : Math.min(8, Math.max(1, Runtime.getRuntime().availableProcessors()));
        }
        threadCount = Math.min(threadCount, characters);
        long[][] coefficients = new long[characters][];
        CharacterWork[] work = new CharacterWork[threadCount];
        for (int i = 0; i < threadCount; ++i) {
            work[i] = new CharacterWork(grundy, coefficients, modulus, i, threadCount);
        }
        // As in the C++, work[0] runs on the calling thread and the others on worker threads.
        ExecutorService pool = threadCount > 1 ? Executors.newFixedThreadPool(threadCount - 1) : null;
        try {
            List<Future<?>> started = new ArrayList<>();
            for (int i = 1; i < threadCount; ++i) {
                CharacterWork job = work[i];
                started.add(pool.submit(() -> calculateCharacters(job)));
            }
            calculateCharacters(work[0]);
            for (Future<?> thread : started) thread.get();
        } catch (InterruptedException error) {
            Thread.currentThread().interrupt();
            throw new IllegalStateException("Interrupted while waiting for character workers", error);
        } catch (ExecutionException error) {
            throw new IllegalStateException("Character worker failed: " + error.getCause(), error.getCause());
        } finally {
            if (pool != null) pool.shutdown();
        }
        for (CharacterWork job : work) {
            if (job.error != null) throw job.error;
        }
        long[] counts = new long[limit / 2 + 1];
        for (int character = 0; character < characters; ++character) {
            for (int n = 0; n < counts.length; ++n) {
                counts[n] = Math.addExact(counts[n], coefficients[character][n]) % modulus;
            }
        }
        for (int n = 0; n < counts.length; ++n) {
            require(counts[n] % characters == 0, "integral Walsh average");
            counts[n] /= characters;
        }
        return counts;
    }

    // Fills powers (length + 1 entries) with 3^i, like the C++ output parameter.
    static byte[] bruteOutcomes(int length, int[] powers) {
        Arrays.fill(powers, 1);
        for (int i = 1; i <= length; ++i) powers[i] = 3 * powers[i - 1];
        byte[] outcomes = new byte[powers[length]];
        for (int state = 1; state < outcomes.length; ++state) {
            int remaining = state;
            int result = 0;
            for (int i = 0; i < length && result != 3; ++i) {
                int pin = remaining % 3;
                remaining /= 3;
                if (pin == 0) continue;
                int single = state - pin * powers[i];
                if (pin == 1 && (outcomes[single] & 2) == 0) result |= 1;
                if (pin == 2 && (outcomes[single] & 1) == 0) result |= 2;
                int next = remaining % 3;
                if (next == 0) continue;
                int pair = single - next * powers[i + 1];
                if ((pin == 1 || next == 1) && (outcomes[pair] & 2) == 0) result |= 1;
                if ((pin == 2 || next == 2) && (outcomes[pair] & 1) == 0) result |= 2;
            }
            outcomes[state] = (byte) result;
        }
        return outcomes;
    }

    static int predictedOutcome(int state, int length, int[] grundy) {
        int balance = 0;
        int previous = 0, run = 0, nimSum = 0;
        for (int i = 0; i < length; ++i) {
            int pin = state % 3;
            state /= 3;
            if (pin == 1) ++balance;
            if (pin == 2) --balance;
            if (pin != 0 && previous != 0 && pin != previous) ++run;
            else {
                nimSum ^= grundy[run];
                run = pin == 0 ? 0 : 1;
            }
            previous = pin;
        }
        nimSum ^= grundy[run];
        if (balance > 0) return 1;
        if (balance < 0) return 2;
        return nimSum == 0 ? 0 : 3;
    }

    private static void runTests() {
        final int SMALL = 12;
        long[] counts = fairCounts(LIMIT);
        require(counts[0] == 1, "empty arrangement");
        require(counts[3] == 8, "given f(6)");
        require(counts[10] == 49_844, "given f(20)");
        require(Arrays.equals(fairCounts(LIMIT, 1), counts), "serial and thread-pool agreement");

        int[] powers = new int[SMALL + 1];
        byte[] outcomes = bruteOutcomes(SMALL, powers);
        int[] grundy = grundyValues(LIMIT);
        for (int state = 0; state < outcomes.length; ++state) {
            require(outcomes[state] == predictedOutcome(state, SMALL, grundy),
                "original-game minimax agrees with balance and alternating blocks");
        }
        for (int n = 0; n <= SMALL; ++n) {
            long expected = 0;
            for (int word = 0; word < (1 << n); ++word) {
                int state = 0;
                for (int i = 0; i < n; ++i) state += (1 + ((word >> i) & 1)) * powers[i];
                expected += outcomes[state] == 0 ? 1 : 0;
            }
            require(expected == (n % 2 == 0 ? counts[n / 2] : 0), "exhaustive initial arrangements");
            long[] prefix = fairCounts(n);
            require(prefix[prefix.length - 1] == counts[n / 2], "short-prefix consistency");
        }

        long[] binomial = new long[LIMIT + 1];
        binomial[0] = 1;
        long[] unsignedCounts = balancedCharacter(grundy, 0, MOD);
        for (int n = 1; n <= LIMIT; ++n) {
            for (int k = n; k > 0; --k) binomial[k] = (binomial[k] + binomial[k - 1]) % MOD;
            if (n % 2 == 0) require(unsignedCounts[n / 2] == binomial[n / 2], "central binomial coefficients");
        }
        System.out.println("All checks passed; " + outcomes.length
            + " original-game positions, given examples, and central binomial coefficients through " + LIMIT + ".");
    }

    public static void main(String[] args) {
        try {
            if (args.length == 1 && args[0].equals("--self-test")) runTests();
            else {
                if (args.length != 0) throw new IllegalArgumentException("Usage: java Euler1015 [--self-test]");
                long[] counts = fairCounts(LIMIT);
                System.out.println(String.format("%09d", counts[counts.length - 1]));
            }
        } catch (RuntimeException error) {
            System.err.println(error.getMessage());
            System.exit(1);
        }
    }
}