Problem 1018: Hex

View on Project Euler

Project Euler Problem 1018 Solution

Every finished Hex board has exactly one winner, so the games that end on a given move can be counted through the numbers C_k of k-cell sets that contain a winning path: White wins on move 2k+1 in M!/(M-2k-1)! (p_(k+1) - p_k) games with p_k = C_k / binom(M, k), and the games Black wins on move 2k+2 are their reflected twins, M - 2k - 1 times as many. Hex duality p_k + p_(M-k) = 1 halves the counting, and a cell-by-cell frontier sweep with connectivity labels counts all sets of at most 60 cells of the 11 x 11 board through 29,964 states. Implementations are available in C++, Python and Java.

Detailed mathematical approach

Problem Summary

Hex is played on a board of \(N\times N\) hexagonal cells in the shape of a rhombus. Two players take turns placing stones on empty cells: the first player places white stones, the second black ones, and a stone is never moved. White wins as soon as the white stones contain a path of adjacent cells joining the two white sides of the board, Black wins as soon as the black stones join the two black sides, and the game cannot end in a draw. Let \(H(N)\) be the number of different games on an \(N\times N\) board, that is, the number of sequences of moves that can be played until one of the players wins. For example \(H(2)=24\), \(H(3)=257760\) and \(H(4)=10663950528000\), and we must find \(H(11)\) modulo \(1\,000\,000\,007\). [1] The game was invented by Piet Hein in 1942 and later rediscovered by John Nash, and \(11\times11\) is its standard board size. [2]

The \(11\times11\) board has \(121\) cells, and the number of orders in which they can be filled, \(121!\), already has \(201\) digits, so the games cannot be listed one by one. The count splits into two parts. The first part is a formula. If \(C_k\) is the number of sets of \(k\) cells that contain a winning path for White, the number of games that end on a given move depends only on the numbers \(C_k\), and every game that White wins gives as many games won by Black one move later as there are empty cells left. Since a full board always has exactly one winner, the symmetry of the board halves the work: only sets of at most half the cells have to be counted. The second part counts these sets for all sizes at once, with a sweep over the board, cell by cell, that only remembers how the last \(N\) cells it has passed are connected to each other and to the top of the board.

Mathematical Approach

We fix coordinates, recall why a full board always has exactly one winner, count the games that end on a given move, and reduce everything to the numbers \(C_k\). The last two sections compute them.

1. Rows, columns and neighbours

Number the rows of the board \(r=0,1,\dots,N-1\) from top to bottom and the columns \(c=0,1,\dots,N-1\) from left to right, so that the cells are the pairs \((r,c)\), and let \(M=N^2\) be their number. Each row of the rhombus is shifted half a cell to the right of the row above it, and the neighbours of the cell \((r,c)\), that is, the cells that share an edge with it, are the cells

$$\begin{aligned}&(r,c-1),\ (r,c+1),\ (r-1,c),\\&(r-1,c+1),\ (r+1,c-1),\ (r+1,c)\end{aligned}$$

that lie on the board: two in the same row, two in the row above and two in the row below. For \(N\ge2\) the corner cells \((0,0)\) and \((N-1,N-1)\) have only two neighbours, and the corners \((0,N-1)\) and \((N-1,0)\) have three. We give White the top row \(r=0\) and the bottom row \(r=N-1\), and Black the left column \(c=0\) and the right column \(c=N-1\). A top–bottom path is a sequence of cells, each adjacent to the next, that starts in the top row and ends in the bottom row; a left–right path runs from the left column to the right column. White wins when the white stones contain a top–bottom path, and Black wins when the black stones contain a left–right path.

The reflection in the long diagonal of the rhombus, \((r,c)\mapsto(c,r)\), maps the board onto itself and keeps adjacency, because it maps the six neighbours of \((r,c)\) listed above onto the six neighbours of \((c,r)\). It exchanges rows and columns, so it turns top–bottom paths into left–right paths and back. Two consequences follow. For every \(k\), as many sets of \(k\) cells contain a top–bottom path as contain a left–right path. And \(H(N)\) does not depend on which pair of sides is white, since reflecting every move turns the games of one convention into the games of the other.

2. A full board has exactly one winner

The whole computation rests on the Hex theorem: if every cell of the board is white or black, then exactly one of two things happens, the white cells contain a top–bottom path or the black cells contain a left–right path. [2] [3]

  • At least one. Gale's proof follows the border between the two colours. Add a white strip above the top row and another below the bottom row, and a black strip on the left and another on the right of the board, so that neighbouring strips meet at the four corners of the board; the cells \((0,N-1)\) and \((N-1,0)\) each have one outer edge that faces a white and a black strip, and it is given to either of the two. Call each cell and each strip a region, and mark every edge of a hexagon of the board that separates a white region from a black one. At each corner point of a hexagon either none or exactly two of the hexagon edges that meet there are marked, since going once around the point the colour changes an even number of times; the exception is the four corners of the board, where a white strip meets a black one without a hexagon edge between them, so that exactly one is marked. The marked edges therefore form closed loops and lines whose ends are corners of the board. Follow the line that starts at the top left corner, between the top strip and the left strip; it never meets itself, so it ends at another corner of the board. Along it, the regions on one side are all white and those on the other side all black, and consecutive regions on the same side are equal or adjacent. So the white regions along the line, together with the top strip, form a connected white set, and the black regions, together with the left strip, a connected black set. If the line ends at the bottom left corner, the white set also reaches the bottom strip; if it ends at the top right corner, the black set reaches the right strip. It cannot end at the bottom right corner: the white regions stay on the same side of the line all the way, on its left as it leaves the top left corner, and at the bottom right corner the white strip would lie on its right. In either case the white cells join the top to the bottom or the black cells join the left to the right.
  • At most one. If the white cells contain a top–bottom path, a shortest one visits no cell twice, and the polygonal line through the centres of its cells, extended from its first cell to the midpoint of one of that cell's edges on the top side and from its last cell to the midpoint of one of that cell's edges on the bottom side, is a simple curve from the top side to the bottom side. By the Jordan curve theorem it separates the left side of the board from the right side, so every left–right path, drawn in the same way as a curve that ends on edges of the left and the right side, has to cross it. On a board of hexagons, cells that touch at all share a whole edge, so the segment between the centres of two neighbours stays inside these two cells, and two paths of cells can only cross in a common cell, which cannot be both white and black; on a chessboard with diagonal moves, by contrast, two paths can cross at a corner point without sharing a square.

The program checks the theorem on all \(2^{16}\) colourings of the \(4\times4\) board, and Section 6 turns it into identities between the numbers \(C_k\). The second half also holds in the middle of a game: there at most one player has a winning path, since the empty cells could be coloured arbitrarily and a path of each colour would survive on the full board. A game therefore ends exactly when the player who has just moved completes a path, and that player wins.

3. Games won by White

For \(0\le k\le M\) let \(C_k\) be the number of sets of \(k\) cells that contain a top–bottom path, and let

$$p_k=\frac{C_k}{\binom{M}{k}}$$

be the probability that \(k\) cells chosen at random contain one. By Section 1, \(C_k\) is also the number of sets of \(k\) cells that contain a left–right path.

A game is a sequence of distinct cells \(x_1,x_2,\dots,x_t\): White places stones on \(x_1,x_3,x_5,\dots\) and Black on \(x_2,x_4,\dots\), no player has a winning path after any of the moves \(1,\dots,t-1\), and one of them has one after move \(t\). Write \(G_t\) for the number of games of length \(t\), so that \(H(N)\) is the sum of all \(G_t\).

A game that ends on move \(t=2k+1\), where \(2k+1\le M\), is won by White with the \((k+1)\)-th white stone. It consists of the white stones \(w_1,\dots,w_{k+1}\) and the black stones \(b_1,\dots,b_k\), in the order in which they were placed, and such a pair of lists is a game exactly when the set \(\{w_1,\dots,w_{k+1}\}\) contains a top–bottom path and the set \(\{w_1,\dots,w_k\}\) does not; then no earlier set of white stones contains one either. Black never had a winning path: the black stones only accumulate, and on the board after the last move, where White has a path, Black has none by Section 2.

Now count. The ordered lists of \(k+1\) cells whose set contains a top–bottom path number \((k+1)!\,C_{k+1}\). Among them, those whose first \(k\) cells already contain one number \(k!\,C_k\,(M-k)\): an ordered list of \(k\) cells with a path, followed by any of the \(M-k\) other cells. The black stones form an ordered list of \(k\) of the remaining \(M-k-1\) cells, which can be chosen in \((M-k-1)!/(M-2k-1)!\) ways. So

$$\begin{aligned}G_{2k+1}&=\Bigl((k+1)!\,C_{k+1}-k!\,(M-k)\,C_k\Bigr)\\&\quad\cdot\frac{(M-k-1)!}{(M-2k-1)!}.\end{aligned}$$

Since \((k+1)!\,(M-k-1)!\,C_{k+1}=M!\,p_{k+1}\) and \(k!\,(M-k)!\,C_k=M!\,p_k\), this is

$$G_{2k+1}=\frac{M!}{(M-2k-1)!}\,\bigl(p_{k+1}-p_k\bigr).$$

The factor \(M!/(M-2k-1)!\) is the number of all sequences of \(2k+1\) distinct cells, and \(p_{k+1}-p_k\) is the fraction of them that are games won by White on move \(2k+1\).

4. Black's wins are twins of White's

A game that ends on move \(2k+2\), where \(2k+2\le M\), is won by Black with the \((k+1)\)-th black stone. Such games come from White's wins on move \(2k+1\) by reflecting the board and letting the players exchange their stones. Write \(x'\) for the image of a cell \(x\) under the reflection \((r,c)\mapsto(c,r)\) of Section 1. Take a game won by White on move \(2k+1\), with white stones \(w_1,\dots,w_{k+1}\) and black stones \(b_1,\dots,b_k\), and a cell \(z\) that is still empty at its end; there are \(M-2k-1\ge1\) choices of \(z\). Form the sequence

$$\begin{aligned}&b_1',\,w_1',\,b_2',\,w_2',\,\dots,\\&b_k',\,w_k',\,z',\,w_{k+1}'.\end{aligned}$$

In it White places \(b_1',\dots,b_k',z'\) and Black places \(w_1',\dots,w_{k+1}'\). Black's \(k+1\) stones contain a left–right path, because \(w_1,\dots,w_{k+1}\) contain a top–bottom path, and Black's first \(k\) stones do not. White's first \(k\) stones contain no top–bottom path, because \(b_1,\dots,b_k\) contain no left–right path, and White's \(k+1\) stones contain none either, since by Section 2 they share the final board with Black's path. So the sequence is a game won by Black on move \(2k+2\). Conversely, every game won by Black on move \(2k+2\) arises in this way from exactly one such pair: reflect it back, and Black's stones give \(w_1,\dots,w_{k+1}\), White's first \(k\) stones give \(b_1,\dots,b_k\), and White's last stone gives \(z\). By the criterion of Section 3 these lists form a game won by White on move \(2k+1\), since Black's \(k+1\) stones contain a left–right path and Black's first \(k\) stones do not, and \(z\) is still empty at its end. Hence

$$\begin{aligned}G_{2k+2}&=(M-2k-1)\,G_{2k+1}\\&=\frac{M!}{(M-2k-2)!}\,\bigl(p_{k+1}-p_k\bigr).\end{aligned}$$

In the language of probability: if all \(M\) cells were filled alternately in an order chosen uniformly at random, then for every \(k\) with \(2k+2\le M\) the first winning path would appear on move \(2k+1\) and on move \(2k+2\) with the same probability \(p_{k+1}-p_k\), because a game of length \(t\) is the beginning of exactly \((M-t)!\) of the \(M!\) orders. When \(M\) is odd, the last move \(M\) belongs to White and has no twin; the first winning path appears on it with probability \(p_{(M+1)/2}-p_{(M-1)/2}\), which is \(1/3\) for \(N=3\).

5. The number of games

Adding both cases for every \(k\), with \(1+(M-2k-1)=M-2k\), gives

$$H(N)=\sum_{k=0}^{\lceil M/2\rceil-1}\frac{M!\,(M-2k)}{(M-2k-1)!}\bigl(p_{k+1}-p_k\bigr).$$

The sum stops at \(k=\lceil M/2\rceil-1\) because the board is full after \(M\) moves. When \(M\) is odd, the last term, with \(2k+1=M\), counts the games of length \(M\), which White wins, and its factor \(M-2k=1\) counts each of them once, since there is no move \(M+1\). The formula needs \(p_0,\dots,p_{\lceil M/2\rceil}\) and nothing else.

For \(N=2\) the board has \(M=4\) cells, and a set of two cells contains a top–bottom path in \(3\) of the \(6\) cases, so \(p_2=1/2\), while \(p_0=p_1=0\). Only the term \(k=1\) contributes, \(24\cdot2\cdot\frac12=24\): \(12\) games won by White with the second white stone and \(12\) won by Black with the second black stone, so \(H(2)=24\). For \(N=3\), with \(M=9\), there are no paths with fewer than three cells, so \(p_0=p_1=p_2=0\), and the table shows \(C_k\) and \(p_k\) for \(k=3,4,5\), which is all that the formula still needs, together with the numbers of games won with the \(k\)-th white stone and with the \(k\)-th black stone. They add up to \(H(3)=257760\).

\(k\)\(C_k\)\(p_k\)\(G_{2k-1}\)\(G_{2k}\)
\(3\)\(8\)\(2/21\)\(1440\)\(5760\)
\(4\)\(42\)\(1/3\)\(43200\)\(86400\)
\(5\)\(84\)\(2/3\)\(120960\)\(0\)

For example \(G_5=9\cdot8\cdot7\cdot6\cdot5\cdot(p_3-p_2)=1440\) and \(G_6=(9-5)\,G_5=5760\); the row of \(k\) uses the difference \(p_k-p_{k-1}\), not \(p_k\) alone.

6. Half of the counts suffice

Colour a set \(S\) of \(k\) cells white and the other \(M-k\) cells black. By Section 2 exactly one of two things happens: \(S\) contains a top–bottom path, or its complement contains a left–right path. Counting over all \(\binom{M}{k}\) sets \(S\), and using that \(C_{M-k}\) also counts the sets of \(M-k\) cells with a left–right path, we get

$$C_k+C_{M-k}=\binom{M}{k},\qquad p_k+p_{M-k}=1.$$

So the counts for more than \(M/2\) cells follow from those for fewer. For \(N=11\) the number \(M=121\) is odd, the formula of Section 5 needs \(p_0,\dots,p_{61}\), and \(p_{61}=1-p_{60}\): only the sets of at most \(60\) cells have to be counted. In the table above, \(p_5=1-p_4\). Added up over all \(k\), the identity says that White wins on exactly half of the \(2^M\) colourings of the full board. When every cell is coloured white or black independently with probability \(1/2\), the white cells therefore contain a top–bottom path with probability exactly \(1/2\), on every \(N\times N\) board. In percolation theory, where the cells become the sites of the triangular lattice, this symmetry is one of the foundations of the theorem that the critical probability of site percolation on that lattice is \(1/2\), a statement about the infinite lattice that also needs further arguments. [4]

Two further facts give quick checks. A top–bottom path has a cell in every row, so \(C_k=0\) for \(k\lt N\). A set of exactly \(N\) cells contains such a path only if it is one, with one cell in each row, and then the cell in row \(r+1\) is \((r+1,c)\) or \((r+1,c-1)\) when the cell in row \(r\) is \((r,c)\). If the path steps to the left \(d\) times, its top cell can be in any of the columns \(d,\dots,N-1\), so for \(N\ge2\)

$$\begin{aligned}C_N&=\sum_{d=0}^{N-1}\binom{N-1}{d}(N-d)\\&=(N+1)\,2^{N-2}.\end{aligned}$$

This is \(8\) for \(N=3\), as in the table, and \(6144\) for \(N=11\). The values \(3,8,20,48,\dots\) for \(N=2,3,4,5,\dots\) are the terms \(a(N-1)\) of the sequence A001792 of the OEIS, \(a(n)=(n+2)\,2^{n-1}\). [5]

7. Sweeping the board

To count the sets, the sweep decides every cell as occupied or empty: first the whole top row at once, then the other cells one at a time in reading order, row by row and from left to right within a row; every set of cells is built exactly once. Just before a cell \((r,c)\) with \(r\ge1\) is decided, call the \(N\) most recently decided cells the frontier: the cells \((r,0),\dots,(r,c-1)\) of the current row followed by the cells \((r-1,c),\dots,(r-1,N-1)\) of the row above, so that position \(j\) of the frontier holds the cell of column \(j\). A decided cell off the frontier has only decided neighbours, and the decided neighbours of \((r,c)\), namely \((r,c-1)\), \((r-1,c)\) and \((r-1,c+1)\) as far as they lie on the board, sit at the frontier positions \(c-1\), \(c\) and \(c+1\).

Of the decided cells, the future only needs to know which frontier cells are occupied and how the occupied ones are connected, to each other and to the top side. For the top side, think of the top strip of Section 2 as an extra occupied region that touches every cell of row \(0\). The sweep stores a state, one label per frontier position: \(0\) for an empty cell, \(1\) for an occupied cell that is connected to the top strip through occupied decided cells, and \(2,3,\dots\) for the other pieces, the groups of occupied frontier cells that are connected to each other through occupied decided cells. These are numbered in the order in which they first appear from the left, so that each pattern of connections has exactly one name. The cells with the label \(1\) form the top piece: they need not be connected to each other inside the board, but the top strip joins them, and since a top–bottom path exists as soon as one of them reaches the bottom row, they need not be told apart. After the top row, every occupied cell has the label \(1\), since it touches the top strip. This is the transfer-matrix method in the form used to count lattice animals: a boundary line moves through the board one cell at a time, and the state records the pattern of occupied cells on it and their connections. [6]

Deciding the cell \((r,c)\) replaces the cell \((r-1,c)\) at position \(c\) by \((r,c)\); the labels at the other positions can only change by merging and renumbering. An empty cell gets the label \(0\). An occupied cell without occupied decided neighbours starts a new piece. Otherwise it joins every piece whose label appears at the positions \(c-1\), \(c\) and \(c+1\) before the change; the merged piece gets the label \(1\) if one of them had it, and every position of the merged piece receives the common label. In every case the labels are then numbered again, since a merge, or a piece that loses its last frontier cell, can leave a gap or change the order of first appearance. For example, let \(N=4\) and let the state be \((2,0,1,0)\) just before the cell \((r,1)\) is decided, so that the frontier holds the cells \((r,0)\), \((r-1,1)\), \((r-1,2)\) and \((r-1,3)\). If \((r,1)\) is occupied, it joins the piece \(2\) to the top piece, and the state becomes \((1,1,1,0)\). If \((r,1)\) is empty, the state stays \((2,0,1,0)\), and if the next cell \((r,2)\) is empty as well, it becomes \((2,0,0,0)\).

In the last case the label \(1\) has disappeared, so the top piece has left the frontier. Its cells now have only decided neighbours, so it can never grow again, and it contains no cell of the bottom row, because those cells stay on the frontier until the end. Every top–bottom path starts in the top piece, so no set built from this state contains one: the state is dropped, and so is every state in which the label \(1\) disappears. For the same reason an empty top row is dropped at once. When the last cell has been decided, the frontier is the bottom row, and a set contains a top–bottom path exactly when a cell of the bottom row carries the label \(1\), which holds for every state that has not been dropped. Each set of cells determines one sequence of states, so \(C_k\) is the number of ways to decide all cells, \(k\) of them occupied, without being dropped.

Two pieces never interleave along the frontier: if positions \(i_1\lt i_2\lt i_3\lt i_4\) held one piece at \(i_1\) and \(i_3\) and another at \(i_2\) and \(i_4\), paths joining \(i_1\) to \(i_3\) and \(i_2\) to \(i_4\) through occupied decided cells, where a path of the top piece may also run through the top strip, would cross. They cannot cross in the top strip, which the other piece does not reach, and on a board of hexagons crossing paths share a cell, so the two pieces would be one. This keeps the number of states small: the program finds \(1,3,9,25,69,189,\)\(518,1422,3915,10813\) and \(29964\) states for \(N=1,\dots,11\), and for \(N=11\) all \(29964\) of them occur together after the cell \((5,6)\). These eleven numbers are also the first terms of the sequence A005322 of the OEIS, a column of the Motzkin triangle; the computation does not depend on this agreement. [7]

8. Counting all sizes at once

The sweep has to count the sets by their number of occupied cells. Each state therefore carries a polynomial in a variable \(x\) whose coefficient of \(x^k\) is the number of ways to reach the state with \(k\) occupied cells. The top row gives the initial states: each non-empty set of \(j\) cells of row \(0\) is a state in which these positions have the label \(1\), and it carries \(x^j\). Deciding a cell passes the polynomial of a state on to the updated state, multiplied by \(x\) if the cell is occupied, and polynomials that arrive at the same state are added. After the last cell the polynomials of all states add up to

$$\sum_{k=0}^{M}C_k\,x^k.$$

By Section 6 only the coefficients of \(x^0,\dots,x^{\lfloor M/2\rfloor}\) are needed, and since a coefficient never moves to a lower power, all higher powers can be discarded at every step: for \(N=11\) each state carries \(61\) numbers. They are kept modulo the prime \(q=1\,000\,000\,007\). Since \(q\gt M\), every factorial up to \(M!\) is invertible modulo \(q\), so \(p_k=C_k\,k!\,(M-k)!/M!\) and the formula of Section 5 can be evaluated modulo \(q\). The result is \(H(N)\) modulo \(q\), because \(H(N)\) is an integer and the formula only divides by numbers prime to \(q\).

How the Code Works

The C++ program [8] stores a state in a \(64\)-bit integer, four bits per frontier position. Four bits are enough: neighbouring frontier positions always hold adjacent cells, so the pieces are separated by empty positions and at most \(\lceil N/2\rceil\) of them meet the frontier; since a live state always contains the top piece, no label exceeds \(6\) for \(N=11\). advance(state, column, side, occupied) performs the update of Section 7: it reads the labels at column - 1, column and column + 1, gives an occupied cell the smallest non-zero label among them, or the unused label \(15\) when all three are \(0\), carries this label to every position of a merged piece, renumbers the labels other than \(0\) and \(1\) as \(2,3,\dots\) from the left, and returns \(0\) when no label \(1\) is left. The value \(0\) marks a dropped state; no live state can be \(0\), since it contains a label \(1\).

The structure Frontier builds all states once. It starts from the \(2^N-1\) states of the first row, which initial keeps together with their numbers of occupied cells, applies advance for every position and both choices to every state found so far, numbers the new states through a hash map and records the results in edges, with NONE for a dropped state. For every column it then turns these transitions around into the lists incoming[column] (offsets and sources), so that each state can collect everything that flows into it.

crossing_counts(graph, threads) runs the sweep of Section 8 with stride \(=\lfloor M/2\rfloor+1\) residues per state. The array values starts with \(x^j\) at the initial states, and for each cell of the rows \(1\) to \(N-1\) every state adds up the polynomials of its sources, shifted by one degree when its own label at the current position shows an occupied cell. Only the degrees up to the number of cells decided so far are touched, since no more of them can be occupied. The states are split into nearly equal ranges, one per thread; after every cell the threads meet at a Barrier built from a pthread mutex and condition variable, and the arrays values and scratch exchange roles; since the number of steps, \(N(N-1)\), is even, the final polynomials end up in values. The calling thread holds the mutex while it starts the others, so that, if a thread cannot be started, the threads that did start share the work. Without an explicit count, or with the count \(0\), the program uses one thread for \(N\lt9\) and otherwise up to eight; explicit counts up to \(64\) are accepted, and no more threads than states are used. At the end the polynomials of all states are added.

game_count(cells, crossings) computes the factorials and their inverses modulo \(q\) with a single Fermat inversion of \(M!\), the residues of \(p_k=C_k\,k!\,(M-k)!/M!\) for \(k\le\lfloor M/2\rfloor\), for odd \(M\) also \(p_{(M+1)/2}=1-p_{(M-1)/2}\), and the sum of Section 5.

The rest serves the self-test. SmallBoard represents a board of at most \(4\times4\) cells by bit masks: wins(stones, player) grows the set of stones connected to the player's first side until it reaches the opposite side or stops growing, and games(white, black, player, memo) counts the games that can still be played from a position by trying every empty cell, with a memo over the pair of stone sets. evaluate_crossings(graph, weight) repeats the sweep with a number in place of \(x\) and without discarding any power, so it returns \(\sum_{k=0}^{M}C_k\,x^k\) modulo \(q\) for that number. run_tests, started with --self-test, takes every board from \(N=1\) to \(N=4\): it counts the sets with a top–bottom path among all \(2^M\) subsets, requiring on the way that every colouring of the full board has exactly one winner, compares the coefficients \(C_0,\dots,C_{\lfloor M/2\rfloor}\) with crossing_counts, checks \(H(1)=1\) and the three values of the statement, and for \(N\le3\) compares \(H(N)\) with a direct count of the game tree. For \(N=11\) it runs consistency checks: it requires that four threads give the same counts as one, that \(C_k=0\) for \(k\lt11\) and that \(C_{11}=12\cdot2^9\); it completes the polynomial with \(C_{121-k}=\binom{121}{k}-C_k\) and compares its values at \(x=1\) and \(x=2\) with evaluate_crossings, the first one also with \(2^{120}\). At \(x=1\) the completed polynomial equals \(\sum_{k=0}^{60}\binom{121}{k}=2^{120}\) whatever the counts are, so these two comparisons only confirm that the sweep without cut-off finds \(2^{120}\) sets with a top–bottom path, as Section 6 requires; the comparison at \(x=2\) is the one that tests the counts. Without arguments the program prints \(H(11)\) modulo \(q\); given a board side from \(1\) to \(11\), and optionally a number of threads, it prints \(H(N)\) modulo \(q\) for that side.

The Python and Java versions follow the C++ program function by function, with the same command line and the same self-test. Java keeps the residues in arrays of int and runs the sweep in threads that meet at a barrier after every cell, as the C++ program does. In Python a loop over single coefficients would be far too slow, so each state keeps its counts as the fields of one integer, each field \(M+1\) bits wide: every count of sets is below \(2^M\), so no field overflows into the next, adding two polynomials is one addition of integers, multiplying by \(x\) is a shift by one field, and discarding the powers above \(x^{\lfloor M/2\rfloor}\) is a mask. The Python counts are therefore exact and are reduced modulo \(q\) only at the end. Its threads follow the same scheme as in C++, but the interpreter runs them one at a time, so they bring no speed-up.

Complexity and Verification

The sweep has \(M-N\) steps, one for each cell below the first row. In each step every transition that is not dropped adds a polynomial of at most \(\lfloor M/2\rfloor+1\) coefficients, so the time is proportional to \(M-N\) times the number of transitions per position times \(M/2\). For \(N=11\) there are between \(57740\) and \(59487\) transitions per position and about \(3.2\cdot10^8\) additions in all, somewhat fewer than \(110\cdot59487\cdot61\approx4.0\cdot10^8\), because while fewer than \(60\) cells have been decided only the degrees up to their number are touched. The memory holds two arrays of \(29964\cdot61\) residues. From \(N=7\) to \(N=11\) the number of states grows by a factor between \(2.74\) and \(2.78\) each time \(N\) grows by one, so the method is exponential in the side of the board, but for \(N=11\) the work is small. Measured as the median wall-clock time of five runs of the whole program with the automatic thread count on an Apple M4 Max, the C++ program built with clang++ -std=c++20 -O3 prints \(H(11)\) in about \(0.14\) seconds, the Java version, compiled beforehand and run on OpenJDK 17, in about \(0.3\) seconds, and the Python version on CPython 3.11 in about \(3.3\) seconds; their self-tests take about \(0.2\), \(0.9\) and \(7\) seconds.

All three programs print the same value for every board side from \(1\) to \(11\) and pass their self-tests. Further checks, which are not part of the published solutions, test the formulas independently of the programs. A second program, written independently from the definition, counts the sets row by row instead of cell by cell: its states give every occupied cell of a row a component and record which components are connected to the top row, it tries all \(2^N\) choices for the next row, and it keeps the complete polynomials as exact integers, without cut-off and without using Section 6. It agrees with a direct count over all \(2^M\) sets for \(N\le5\), which also confirms that exactly one player wins on each of the \(2^{25}\) colourings of the \(5\times5\) board. For every \(N\le11\) its coefficients \(C_0,\dots,C_{\lfloor M/2\rfloor}\) agree with the residues of the sweep, all its coefficients satisfy \(C_k+C_{M-k}=\binom{M}{k}\), the identity that Section 6 derives from the Hex theorem, and they give \(C_k=0\) for \(k\lt N\) and \(C_N=(N+1)\,2^{N-2}\). When the components connected to the top row are merged, the states it reaches after complete rows, without the dropped ones, give exactly the set of states of the sweep. From these exact counts \(H(N)\) was computed in exact rational arithmetic for every \(N\le11\): every result is an integer, the values for \(N=2,3,4\) are those of the statement, and the residues agree with the programs. Finally, an exhaustive count over the positions of the game, which knows nothing of the formula, finds the number of games of each length on the boards up to \(4\times4\). For \(N=4\) it reaches about \(9.4\) million positions and the \(10663950528000\) games of the statement, and for every \(N\le4\) the numbers of games of each length agree with Sections 3 to 5, including \(G_{2k+2}=(M-2k-1)\,G_{2k+1}\).

Footnotes and References

The references below supply background facts. The game-count formula, the correspondence between the wins of the two players, the halving of the counts and the sweep are derived explicitly in this article.

  1. Project Euler 1018 — Hex. The official statement: the rules of Hex on a rhombus of hexagonal cells, the definition of \(H(N)\) as the number of possible sequences of moves until a player wins, the values \(H(2)=24\), \(H(3)=257760\) and \(H(4)=10663950528000\), and the task to find \(H(11)\) modulo \(1\,000\,000\,007\).
  2. Wikipedia — Hex (board game). Hex is played on a rhombus of hexagonal cells, traditionally of size \(11\times11\). The article reports that Piet Hein invented the game in 1942 and that John Nash rediscovered it at Princeton in 1948 or 1949, and it states the hex theorem: however the board is filled with stones, exactly one player has connected the two sides assigned to them, so the game cannot end in a draw. It also mentions Gale's proof that this property is equivalent to the Brouwer fixed-point theorem.
  3. David Gale — The Game of Hex and the Brouwer Fixed-Point Theorem, The American Mathematical Monthly 86 (1979), 818–827. Gale proves that every completely marked board contains a winning path, by walking from a corner of the board along the edges that separate the two players' regions, and shows that this property of Hex is equivalent to the Brouwer fixed-point theorem. He remarks that the stronger statement, that the two players cannot both have a winning path, is connected with the Jordan curve theorem; his arguments only use the first statement.
  4. Wendelin Werner — Lectures on two-dimensional critical percolation, IAS/Park City Mathematics Series 16 (2009), 297–360. Section 1.2 explains that colouring the sites of the triangular lattice is the same as colouring the cells of a board of hexagons, that in a rhombus exactly one of a white top-to-bottom crossing and a black left-to-right crossing occurs, so that both have probability \(1/2\) when each site is white with probability \(1/2\), and that the critical probability \(1/2\) of this lattice, a result of Kesten and Wierman, builds on such symmetry properties but also requires additional arguments.
  5. OEIS A001792. The sequence \(a(n)=(n+2)\,2^{n-1}\) for \(n\ge0\), which begins \(1,3,8,20,48,\dots\), so that \(C_N=a(N-1)\). A comment in the entry identifies these numbers with the numbers of shortest winning paths across a Hex board.
  6. Iwan Jensen — Enumerations of lattice animals and trees, Journal of Statistical Physics 102 (2001), 865–881. Section 2.1 describes the transfer-matrix algorithm for counting lattice animals: a boundary line cuts the lattice, a generating function of the partially built configurations is kept for every pattern of occupied and empty sites on the line, together with their connections through occupied sites behind it, and the line moves forward one site at a time.
  7. OEIS A005322. A column of the Motzkin triangle, beginning \(1,3,9,25,69,189,518,\dots\); one of its comments counts the paths of length \(n\) with up, down and level steps that never go below their starting level and end two levels higher. Its first eleven terms are the numbers of states found in Section 7 for \(N=1,\dots,11\); this article neither proves nor uses a connection between the two.
  8. C++ — Euler1018.cpp. The linked immutable C++ revision contains advance, Frontier, crossing_counts, game_count, SmallBoard, evaluate_crossings and the checks run by --self-test. It is the implementation record for this article.

Problem 1018 source code

C++

#include <algorithm>
#include <array>
#include <charconv>
#include <cstdint>
#include <cstdlib>
#include <iostream>
#include <limits>
#include <pthread.h>
#include <stdexcept>
#include <string>
#include <thread>
#include <unordered_map>
#include <utility>
#include <vector>

namespace {

using u32 = std::uint32_t;
using u64 = std::uint64_t;
constexpr u32 MOD = 1'000'000'007;
constexpr unsigned MAX_SIDE = 11;
constexpr u32 NONE = std::numeric_limits<u32>::max();

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

u32 add(const u32 a, const u32 b) {
    const u32 sum = a + b;
    return sum >= MOD ? sum - MOD : sum;
}

u32 multiply(const u32 a, const u32 b) {
    return static_cast<u32>(u64{a} * b % MOD);
}

u32 power(u32 base, u32 exponent) {
    u32 result = 1;
    for (; exponent != 0; exponent >>= 1, base = multiply(base, base)) {
        if (exponent & 1) result = multiply(result, base);
    }
    return result;
}

u64 advance(const u64 state, const unsigned column, const unsigned side, const bool occupied) {
    std::array<unsigned, MAX_SIDE> labels{};
    for (unsigned i = 0; i < side; ++i) labels[i] = static_cast<unsigned>((state >> (4 * i)) & 15);
    if (!occupied) {
        labels[column] = 0;
    } else {
        const std::array<unsigned, 3> adjacent{
            column == 0 ? 0 : labels[column - 1], labels[column],
            column + 1 == side ? 0 : labels[column + 1]
        };
        unsigned merged = 15;
        for (const unsigned label : adjacent) {
            if (label != 0) merged = std::min(merged, label);
        }
        for (unsigned i = 0; i < side; ++i) {
            for (const unsigned label : adjacent) {
                if (label != 0 && labels[i] == label) labels[i] = merged;
            }
        }
        labels[column] = merged;
    }

    // Label 1 connects to the top edge; other component labels are canonicalized.
    std::array<unsigned, 16> canonical{0, 1};
    unsigned next = 2;
    bool reaches_top = false;
    u64 result = 0;
    for (unsigned i = 0; i < side; ++i) {
        reaches_top |= labels[i] == 1;
        if (labels[i] != 0 && canonical[labels[i]] == 0) canonical[labels[i]] = next++;
        result |= u64{canonical[labels[i]]} << (4 * i);
    }
    return reaches_top ? result : 0;
}

struct Incoming {
    std::vector<u32> offsets;
    std::vector<u32> sources;
};

struct Frontier {
    unsigned side;
    std::vector<u64> states;
    std::vector<std::pair<u32, unsigned>> initial;
    std::vector<std::array<u32, 2 * MAX_SIDE>> edges;
    std::vector<Incoming> incoming;

    explicit Frontier(const unsigned width) : side(width) {
        require(side >= 1 && side <= MAX_SIDE, "board side between 1 and 11");
        incoming.resize(side);
        std::unordered_map<u64, u32> ids;
        ids.reserve(std::size_t{1} << (side + 4));
        const auto intern = [&](const u64 state) {
            const auto entry = ids.emplace(state, static_cast<u32>(states.size()));
            if (entry.second) states.push_back(state);
            return entry.first->second;
        };
        for (unsigned mask = 1; mask < (1U << side); ++mask) {
            u64 state = 0;
            unsigned count = 0;
            for (unsigned column = 0; column < side; ++column) {
                if (mask & (1U << column)) {
                    state |= u64{1} << (4 * column);
                    ++count;
                }
            }
            initial.emplace_back(intern(state), count);
        }
        for (std::size_t i = 0; i < states.size(); ++i) {
            const u64 state = states[i];
            std::array<u32, 2 * MAX_SIDE> transitions{};
            for (unsigned column = 0; column < side; ++column) {
                for (unsigned occupied = 0; occupied < 2; ++occupied) {
                    const u64 next = advance(state, column, side, occupied != 0);
                    transitions[2 * column + occupied] = next == 0 ? NONE : intern(next);
                }
            }
            edges.push_back(transitions);
        }
        for (unsigned column = 0; column < side; ++column) {
            auto& reverse = incoming[column];
            reverse.offsets.resize(states.size() + 1);
            for (const auto& transitions : edges) {
                for (unsigned occupied = 0; occupied < 2; ++occupied) {
                    const u32 target = transitions[2 * column + occupied];
                    if (target != NONE) ++reverse.offsets[target + 1];
                }
            }
            for (std::size_t i = 0; i < states.size(); ++i) reverse.offsets[i + 1] += reverse.offsets[i];
            reverse.sources.resize(reverse.offsets.back());
            auto positions = reverse.offsets;
            for (u32 source = 0; source < states.size(); ++source) {
                for (unsigned occupied = 0; occupied < 2; ++occupied) {
                    const u32 target = edges[source][2 * column + occupied];
                    if (target != NONE) reverse.sources[positions[target]++] = source;
                }
            }
        }
    }
};

struct Barrier {
    pthread_mutex_t mutex = PTHREAD_MUTEX_INITIALIZER;
    pthread_cond_t changed = PTHREAD_COND_INITIALIZER;
    unsigned arrived = 0;
    unsigned generation = 0;

    ~Barrier() {
        pthread_cond_destroy(&changed);
        pthread_mutex_destroy(&mutex);
    }

    void wait(const unsigned participants) {
        if (participants == 1) return;
        pthread_mutex_lock(&mutex);
        const unsigned previous = generation;
        if (++arrived == participants) {
            arrived = 0;
            ++generation;
            pthread_cond_broadcast(&changed);
        } else {
            while (generation == previous) pthread_cond_wait(&changed, &mutex);
        }
        pthread_mutex_unlock(&mutex);
    }
};

struct Work {
    const Frontier& frontier;
    u32* values;
    u32* scratch;
    unsigned stride;
    unsigned threads;
    Barrier barrier;

    Work(const Frontier& graph, u32* first, u32* second, const unsigned length, const unsigned count)
        : frontier(graph), values(first), scratch(second), stride(length), threads(count) {}
};

struct Task {
    Work* work;
    unsigned index;
};

void* count_worker(void* argument) {
    const auto& task = *static_cast<Task*>(argument);
    auto& work = *task.work;
    pthread_mutex_lock(&work.barrier.mutex);
    pthread_mutex_unlock(&work.barrier.mutex);
    const auto& graph = work.frontier;
    const std::size_t first = graph.states.size() * task.index / work.threads;
    const std::size_t last = graph.states.size() * (task.index + 1) / work.threads;
    u32* current = work.values;
    u32* next = work.scratch;
    for (unsigned cell = graph.side; cell < graph.side * graph.side; ++cell) {
        const unsigned column = cell % graph.side;
        const unsigned upper = std::min(work.stride - 1, cell + 1);
        const auto& incoming = graph.incoming[column];
        for (std::size_t target = first; target < last; ++target) {
            const unsigned occupied = ((graph.states[target] >> (4 * column)) & 15) != 0;
            u32* output = next + target * work.stride;
            std::fill(output, output + upper + 1, 0);
            for (u32 j = incoming.offsets[target]; j < incoming.offsets[target + 1]; ++j) {
                const u32* input = current + std::size_t{incoming.sources[j]} * work.stride;
                for (unsigned degree = occupied; degree <= upper; ++degree) {
                    output[degree] = add(output[degree], input[degree - occupied]);
                }
            }
        }
        work.barrier.wait(work.threads);
        std::swap(current, next);
    }
    return nullptr;
}

std::vector<u32> crossing_counts(const Frontier& graph, unsigned threads) {
    const unsigned stride = graph.side * graph.side / 2 + 1;
    std::vector<u32> values(graph.states.size() * stride), scratch(values.size());
    for (const auto& entry : graph.initial) {
        if (entry.second < stride) values[std::size_t{entry.first} * stride + entry.second] = 1;
    }
    if (threads == 0) {
        threads = graph.side < 9 ? 1 : std::min(8U, std::max(1U, std::thread::hardware_concurrency()));
    }
    threads = static_cast<unsigned>(std::min<std::size_t>(threads, graph.states.size()));
    Work work(graph, values.data(), scratch.data(), stride, threads);
    std::vector<Task> tasks;
    for (unsigned i = 0; i < threads; ++i) tasks.push_back({&work, i});
    std::vector<pthread_t> workers(threads - 1);
    pthread_mutex_lock(&work.barrier.mutex);
    unsigned started = 1;
    for (; started < threads; ++started) {
        if (pthread_create(&workers[started - 1], nullptr, count_worker, &tasks[started]) != 0) break;
    }
    work.threads = started;
    pthread_mutex_unlock(&work.barrier.mutex);
    count_worker(&tasks[0]);
    bool joined = true;
    for (unsigned i = 1; i < started; ++i) {
        if (pthread_join(workers[i - 1], nullptr) != 0) joined = false;
    }
    require(joined, "pthread join");
    std::vector<u32> result(stride);
    for (std::size_t state = 0; state < graph.states.size(); ++state) {
        for (unsigned degree = 0; degree < stride; ++degree) {
            result[degree] = add(result[degree], values[state * stride + degree]);
        }
    }
    return result;
}

u32 game_count(const unsigned cells, const std::vector<u32>& crossings) {
    std::vector<u32> factorial(cells + 1, 1), inverse(cells + 1, 1), probability(cells / 2 + 2);
    for (unsigned i = 1; i <= cells; ++i) factorial[i] = multiply(factorial[i - 1], i);
    inverse[cells] = power(factorial[cells], MOD - 2);
    for (unsigned i = cells; i != 0; --i) inverse[i - 1] = multiply(inverse[i], i);
    for (unsigned k = 0; k <= cells / 2; ++k) {
        probability[k] = multiply(crossings[k], multiply(inverse[cells],
            multiply(factorial[k], factorial[cells - k])));
    }
    if (cells % 2 != 0) probability[cells / 2 + 1] = (1 + MOD - probability[cells / 2]) % MOD;

    // With p_k the crossing probability, wins on moves 2k+1 and 2k+2 total
    // (cells)_(2k+1) * (cells - 2k) * (p_(k+1) - p_k).
    u32 result = 0;
    for (unsigned k = 0; k < (cells + 1) / 2; ++k) {
        const u32 difference = (probability[k + 1] + MOD - probability[k]) % MOD;
        const u32 weight = multiply(multiply(factorial[cells], inverse[cells - 2 * k - 1]), cells - 2 * k);
        result = add(result, multiply(weight, difference));
    }
    return result;
}

struct SmallBoard {
    unsigned side;
    unsigned cells;
    std::vector<u32> adjacent;
    std::array<u32, 2> first{}, last{};

    explicit SmallBoard(const unsigned n) : side(n), cells(n * n), adjacent(cells) {
        constexpr std::array<int, 6> DR{-1, -1, 0, 0, 1, 1};
        constexpr std::array<int, 6> DC{0, 1, -1, 1, -1, 0};
        for (unsigned cell = 0; cell < cells; ++cell) {
            const unsigned row = cell / side, column = cell % side;
            const u32 bit = u32{1} << cell;
            if (row == 0) first[0] |= bit;
            if (row + 1 == side) last[0] |= bit;
            if (column == 0) first[1] |= bit;
            if (column + 1 == side) last[1] |= bit;
            for (unsigned d = 0; d < DR.size(); ++d) {
                const int r = static_cast<int>(row) + DR[d], c = static_cast<int>(column) + DC[d];
                if (r >= 0 && c >= 0 && r < static_cast<int>(side) && c < static_cast<int>(side)) {
                    adjacent[cell] |= u32{1} << (static_cast<unsigned>(r) * side + static_cast<unsigned>(c));
                }
            }
        }
    }

    bool wins(const u32 stones, const unsigned player) const {
        u32 connected = stones & first[player];
        for (;;) {
            if (connected & last[player]) return true;
            u32 next = connected;
            for (unsigned cell = 0; cell < cells; ++cell) {
                if (connected & (u32{1} << cell)) next |= adjacent[cell] & stones;
            }
            if (next == connected) return false;
            connected = next;
        }
    }

    u64 games(const u32 white, const u32 black, const unsigned player, std::vector<u64>& memo) const {
        u64& result = memo[white | (black << cells)];
        if (result != std::numeric_limits<u64>::max()) return result;
        if (wins(white, 0) || wins(black, 1)) return result = 1;
        result = 0;
        for (unsigned cell = 0; cell < cells; ++cell) {
            const u32 bit = u32{1} << cell;
            if ((white | black) & bit) continue;
            result += player == 0 ? games(white | bit, black, 1, memo) : games(white, black | bit, 0, memo);
        }
        return result;
    }
};

u32 evaluate_crossings(const Frontier& graph, const u32 weight) {
    std::vector<u32> current(graph.states.size()), next(current.size());
    for (const auto& entry : graph.initial) current[entry.first] = power(weight, entry.second);
    for (unsigned cell = graph.side; cell < graph.side * graph.side; ++cell) {
        std::fill(next.begin(), next.end(), 0);
        const unsigned column = cell % graph.side;
        for (std::size_t source = 0; source < current.size(); ++source) {
            for (unsigned occupied = 0; occupied < 2; ++occupied) {
                const u32 target = graph.edges[source][2 * column + occupied];
                if (target == NONE) continue;
                next[target] = add(next[target], occupied == 0 ? current[source] : multiply(current[source], weight));
            }
        }
        current.swap(next);
    }
    u32 result = 0;
    for (const u32 value : current) result = add(result, value);
    return result;
}

void run_tests() {
    constexpr std::array<u64, 4> GIVEN{1, 24, 257'760, 10'663'950'528'000};
    for (unsigned side = 1; side <= 4; ++side) {
        const Frontier graph(side);
        const auto actual = crossing_counts(graph, 1);
        const SmallBoard board(side);
        std::vector<u32> expected(board.cells + 1);
        const u32 all = (u32{1} << board.cells) - 1;
        for (u32 stones = 0; stones <= all; ++stones) {
            const bool white_wins = board.wins(stones, 0);
            require(white_wins != board.wins(all ^ stones, 1), "exactly one winner on a full board");
            unsigned count = 0;
            for (u32 remaining = stones; remaining != 0; remaining &= remaining - 1) ++count;
            if (white_wins) ++expected[count];
        }
        require(std::equal(actual.begin(), actual.end(), expected.begin()), "all subsets through side 4");
        const u32 answer = game_count(board.cells, actual);
        require(answer == GIVEN[side - 1] % MOD, "given game count");
        if (side <= 3) {
            std::vector<u64> memo(std::size_t{1} << (2 * board.cells), std::numeric_limits<u64>::max());
            require(answer == board.games(0, 0, 0, memo) % MOD, "direct game tree through side 3");
        }
    }

    const Frontier graph(MAX_SIDE);
    const auto serial = crossing_counts(graph, 1);
    require(serial == crossing_counts(graph, 4), "serial and pthread agreement");
    for (unsigned k = 0; k < MAX_SIDE; ++k) require(serial[k] == 0, "at least one stone in every row");
    require(serial[MAX_SIDE] == (MAX_SIDE + 1) * (1U << (MAX_SIDE - 2)), "minimal crossing paths");

    const unsigned cells = MAX_SIDE * MAX_SIDE;
    std::vector<u32> binomial(cells + 1), full(cells + 1);
    binomial[0] = 1;
    for (unsigned row = 1; row <= cells; ++row) {
        for (unsigned k = row; k != 0; --k) binomial[k] = add(binomial[k], binomial[k - 1]);
    }
    for (unsigned k = 0; k <= cells / 2; ++k) {
        full[k] = serial[k];
        full[cells - k] = (binomial[k] + MOD - serial[k]) % MOD;
    }
    for (const u32 weight : {1U, 2U}) {
        u32 value = 0;
        for (unsigned k = cells + 1; k != 0; --k) value = add(multiply(value, weight), full[k - 1]);
        require(value == evaluate_crossings(graph, weight), "full polynomial and weighted forward transfer");
        if (weight == 1) require(value == power(2, cells - 1), "half of all full boards connect each side");
    }
    std::cout << "All checks passed; supplied examples, exhaustive subsets through side 4, "
                 "game trees through side 3, pthread agreement, and weighted transfers at side 11.\n";
}

unsigned parse_unsigned(const char* text) {
    const std::string value(text);
    unsigned result = 0;
    const auto parsed = std::from_chars(value.data(), value.data() + value.size(), result);
    if (parsed.ec != std::errc{} || parsed.ptr != value.data() + value.size()) {
        throw std::invalid_argument("Expected an unsigned integer");
    }
    return result;
}

}

int main(int argc, char* argv[]) {
    try {
        if (argc == 2 && std::string(argv[1]) == "--self-test") {
            run_tests();
        } else {
            if (argc > 3) throw std::invalid_argument("Usage: Euler1018 [side [threads]] | --self-test");
            const unsigned side = argc > 1 ? parse_unsigned(argv[1]) : MAX_SIDE;
            const unsigned threads = argc > 2 ? parse_unsigned(argv[2]) : 0;
            require(threads <= 64, "at most 64 threads, or 0 for automatic selection");
            const Frontier graph(side);
            std::cout << game_count(side * side, crossing_counts(graph, threads)) << '\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 1018 - Hex."""

import os
import sys
import threading

MOD = 1_000_000_007
MAX_SIDE = 11
NONE = (1 << 32) - 1


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


def add(a, b):
    total = a + b
    return total - MOD if total >= MOD else total


def multiply(a, b):
    return a * b % MOD


def power(base, exponent):
    result = 1
    while exponent != 0:
        if exponent & 1:
            result = multiply(result, base)
        exponent >>= 1
        base = multiply(base, base)
    return result


def advance(state, column, side, occupied):
    labels = [(state >> (4 * i)) & 15 for i in range(side)]
    if not occupied:
        labels[column] = 0
    else:
        adjacent = (0 if column == 0 else labels[column - 1], labels[column],
                    0 if column + 1 == side else labels[column + 1])
        merged = 15
        for label in adjacent:
            if label != 0:
                merged = min(merged, label)
        for i in range(side):
            if labels[i] != 0 and labels[i] in adjacent:
                labels[i] = merged
        labels[column] = merged

    # Label 1 connects to the top edge; other component labels are canonicalized.
    canonical = [0, 1] + [0] * 14
    next_label = 2
    reaches_top = False
    result = 0
    for i in range(side):
        reaches_top |= labels[i] == 1
        if labels[i] != 0 and canonical[labels[i]] == 0:
            canonical[labels[i]] = next_label
            next_label += 1
        result |= canonical[labels[i]] << (4 * i)
    return result if reaches_top else 0


class Incoming:
    __slots__ = ("offsets", "sources")

    def __init__(self):
        self.offsets = []
        self.sources = []


class Frontier:
    def __init__(self, width):
        self.side = width
        self.states = []
        self.initial = []
        self.edges = []
        require(1 <= self.side <= MAX_SIDE, "board side between 1 and 11")
        side = self.side
        states = self.states
        self.incoming = [Incoming() for _ in range(side)]
        ids = {}

        def intern(state):
            identifier = ids.setdefault(state, len(states))
            if identifier == len(states):
                states.append(state)
            return identifier

        for mask in range(1, 1 << side):
            state = 0
            count = 0
            for column in range(side):
                if mask & (1 << column):
                    state |= 1 << (4 * column)
                    count += 1
            self.initial.append((intern(state), count))
        i = 0
        while i < len(states):          # intern() appends to the list being scanned
            transitions = [0] * (2 * MAX_SIDE)
            for column in range(side):
                for occupied in range(2):
                    next_state = advance(states[i], column, side, occupied != 0)
                    transitions[2 * column + occupied] = (NONE if next_state == 0
                                                          else intern(next_state))
            self.edges.append(transitions)
            i += 1
        for column in range(side):
            reverse = self.incoming[column]
            reverse.offsets = [0] * (len(states) + 1)
            for transitions in self.edges:
                for occupied in range(2):
                    target = transitions[2 * column + occupied]
                    if target != NONE:
                        reverse.offsets[target + 1] += 1
            for i in range(len(states)):
                reverse.offsets[i + 1] += reverse.offsets[i]
            reverse.sources = [0] * reverse.offsets[-1]
            positions = reverse.offsets[:]
            for source in range(len(states)):
                for occupied in range(2):
                    target = self.edges[source][2 * column + occupied]
                    if target != NONE:
                        reverse.sources[positions[target]] = source
                        positions[target] += 1


class Work:
    __slots__ = ("frontier", "values", "scratch", "stride", "threads", "gate", "barrier")

    def __init__(self, graph, first, second, length, count):
        self.frontier = graph
        self.values = first
        self.scratch = second
        self.stride = length
        self.threads = count
        # The C++ barrier mutex doubles as this start gate; threading.Barrier needs its
        # party count up front, so it is made once the caller knows how many threads started.
        self.gate = threading.Lock()
        self.barrier = None


class Task:
    __slots__ = ("work", "index", "error")

    def __init__(self, work, index):
        self.work = work
        self.index = index
        self.error = None


def count_worker(task):
    work = task.work
    try:
        with work.gate:
            pass
        if work.barrier is None:        # the caller failed while starting the threads
            return
        graph = work.frontier
        first = len(graph.states) * task.index // work.threads
        last = len(graph.states) * (task.index + 1) // work.threads
        current = work.values
        following = work.scratch
        # Python keeps each truncated polynomial in one integer: field k, cells + 1 bits wide,
        # holds the exact count of degree k (every count of subsets is below 2^cells, so no
        # field overflows). An occupied cell shifts by one field; the mask drops degree stride.
        field = graph.side * graph.side + 1
        mask = (1 << (field * work.stride)) - 1
        states = graph.states
        for cell in range(graph.side, graph.side * graph.side):
            column = cell % graph.side
            offsets = graph.incoming[column].offsets
            sources = graph.incoming[column].sources
            for target in range(first, last):
                total = 0
                for j in range(offsets[target], offsets[target + 1]):
                    total += current[sources[j]]
                if (states[target] >> (4 * column)) & 15:
                    total = (total << field) & mask
                following[target] = total
            work.barrier.wait()
            current, following = following, current
    except BaseException as error:      # a thread cannot raise into its caller
        task.error = error
        work.barrier.abort()            # release the others instead of leaving them waiting


def crossing_counts(graph, threads):
    cells = graph.side * graph.side
    stride = cells // 2 + 1
    field = cells + 1
    values = [0] * len(graph.states)
    scratch = [0] * len(values)
    for state, count in graph.initial:
        if count < stride:
            values[state] |= 1 << (field * count)
    if threads == 0:
        threads = 1 if graph.side < 9 else min(8, max(1, os.cpu_count() or 0))
    threads = min(threads, len(graph.states))
    work = Work(graph, values, scratch, stride, threads)
    tasks = [Task(work, i) for i in range(threads)]
    workers = []
    with work.gate:
        for task in tasks[1:]:
            worker = threading.Thread(target=count_worker, args=(task,))
            try:
                worker.start()
            except RuntimeError:        # like a failed pthread_create: go on with those started
                break
            workers.append(worker)
        work.threads = 1 + len(workers)
        work.barrier = threading.Barrier(work.threads)
    count_worker(tasks[0])
    for worker in workers:
        worker.join()
    errors = [task.error for task in tasks if task.error is not None]
    if errors:
        # The first failure aborts the barrier; the others then stop with BrokenBarrierError.
        raise next((error for error in errors
                    if not isinstance(error, threading.BrokenBarrierError)), errors[0])
    total = sum(values)
    field_mask = (1 << field) - 1
    return [((total >> (field * degree)) & field_mask) % MOD for degree in range(stride)]


def game_count(cells, crossings):
    factorial = [1] * (cells + 1)
    inverse = [1] * (cells + 1)
    probability = [0] * (cells // 2 + 2)
    for i in range(1, cells + 1):
        factorial[i] = multiply(factorial[i - 1], i)
    inverse[cells] = power(factorial[cells], MOD - 2)
    for i in range(cells, 0, -1):
        inverse[i - 1] = multiply(inverse[i], i)
    for k in range(cells // 2 + 1):
        probability[k] = multiply(
            crossings[k], multiply(inverse[cells], multiply(factorial[k], factorial[cells - k])))
    if cells % 2 != 0:
        probability[cells // 2 + 1] = (1 + MOD - probability[cells // 2]) % MOD

    # With p_k the crossing probability, wins on moves 2k+1 and 2k+2 total
    # (cells)_(2k+1) * (cells - 2k) * (p_(k+1) - p_k).
    result = 0
    for k in range((cells + 1) // 2):
        difference = (probability[k + 1] + MOD - probability[k]) % MOD
        weight = multiply(multiply(factorial[cells], inverse[cells - 2 * k - 1]), cells - 2 * k)
        result = add(result, multiply(weight, difference))
    return result


class SmallBoard:
    def __init__(self, n):
        self.side = n
        self.cells = n * n
        self.adjacent = [0] * self.cells
        self.first = [0, 0]
        self.last = [0, 0]
        dr = (-1, -1, 0, 0, 1, 1)
        dc = (0, 1, -1, 1, -1, 0)
        for cell in range(self.cells):
            row, column = divmod(cell, self.side)
            bit = 1 << cell
            if row == 0:
                self.first[0] |= bit
            if row + 1 == self.side:
                self.last[0] |= bit
            if column == 0:
                self.first[1] |= bit
            if column + 1 == self.side:
                self.last[1] |= bit
            for d in range(len(dr)):
                r, c = row + dr[d], column + dc[d]
                if 0 <= r < self.side and 0 <= c < self.side:
                    self.adjacent[cell] |= 1 << (r * self.side + c)

    def wins(self, stones, player):
        connected = stones & self.first[player]
        while True:
            if connected & self.last[player]:
                return True
            reached = connected
            for cell in range(self.cells):
                if connected & (1 << cell):
                    reached |= self.adjacent[cell] & stones
            if reached == connected:
                return False
            connected = reached

    def games(self, white, black, player, memo):
        key = white | (black << self.cells)
        if memo[key] is not None:
            return memo[key]
        if self.wins(white, 0) or self.wins(black, 1):
            memo[key] = 1
            return 1
        result = 0
        for cell in range(self.cells):
            bit = 1 << cell
            if (white | black) & bit:
                continue
            result += (self.games(white | bit, black, 1, memo) if player == 0
                       else self.games(white, black | bit, 0, memo))
        memo[key] = result
        return result


def evaluate_crossings(graph, weight):
    current = [0] * len(graph.states)
    for state, count in graph.initial:
        current[state] = power(weight, count)
    for cell in range(graph.side, graph.side * graph.side):
        following = [0] * len(current)
        column = cell % graph.side
        for source in range(len(current)):
            for occupied in range(2):
                target = graph.edges[source][2 * column + occupied]
                if target == NONE:
                    continue
                following[target] = add(following[target], current[source] if occupied == 0
                                        else multiply(current[source], weight))
        current = following
    result = 0
    for value in current:
        result = add(result, value)
    return result


def run_tests():
    given = (1, 24, 257_760, 10_663_950_528_000)
    for side in range(1, 5):
        graph = Frontier(side)
        actual = crossing_counts(graph, 1)
        board = SmallBoard(side)
        expected = [0] * (board.cells + 1)
        all_cells = (1 << board.cells) - 1
        for stones in range(all_cells + 1):
            white_wins = board.wins(stones, 0)
            require(white_wins != board.wins(all_cells ^ stones, 1),
                    "exactly one winner on a full board")
            if white_wins:
                expected[stones.bit_count()] += 1
        require(actual == expected[:len(actual)], "all subsets through side 4")
        answer = game_count(board.cells, actual)
        require(answer == given[side - 1] % MOD, "given game count")
        if side <= 3:
            memo = [None] * (1 << (2 * board.cells))
            require(answer == board.games(0, 0, 0, memo) % MOD, "direct game tree through side 3")

    graph = Frontier(MAX_SIDE)
    serial = crossing_counts(graph, 1)
    require(serial == crossing_counts(graph, 4), "serial and thread agreement")
    for k in range(MAX_SIDE):
        require(serial[k] == 0, "at least one stone in every row")
    require(serial[MAX_SIDE] == (MAX_SIDE + 1) * (1 << (MAX_SIDE - 2)), "minimal crossing paths")

    cells = MAX_SIDE * MAX_SIDE
    binomial = [0] * (cells + 1)
    full = [0] * (cells + 1)
    binomial[0] = 1
    for row in range(1, cells + 1):
        for k in range(row, 0, -1):
            binomial[k] = add(binomial[k], binomial[k - 1])
    for k in range(cells // 2 + 1):
        full[k] = serial[k]
        full[cells - k] = (binomial[k] + MOD - serial[k]) % MOD
    for weight in (1, 2):
        value = 0
        for k in range(cells + 1, 0, -1):
            value = add(multiply(value, weight), full[k - 1])
        require(value == evaluate_crossings(graph, weight),
                "full polynomial and weighted forward transfer")
        if weight == 1:
            require(value == power(2, cells - 1), "half of all full boards connect each side")
    print("All checks passed; supplied examples, exhaustive subsets through side 4, "
          "game trees through side 3, thread agreement, and weighted transfers at side 11.")


# Like std::from_chars into an unsigned 32-bit value: ASCII digits only, no sign or spaces.
def parse_unsigned(text):
    digits = text.lstrip("0") or "0"
    if (not text or any(c not in "0123456789" for c in text)
            or len(digits) > 10 or int(digits) > (1 << 32) - 1):
        raise ValueError("Expected an unsigned integer")
    return int(digits)


def main(argv):
    try:
        if len(argv) == 2 and argv[1] == "--self-test":
            run_tests()
        else:
            if len(argv) > 3:
                raise ValueError("Usage: Euler1018.py [side [threads]] | --self-test")
            side = parse_unsigned(argv[1]) if len(argv) > 1 else MAX_SIDE
            threads = parse_unsigned(argv[2]) if len(argv) > 2 else 0
            require(threads <= 64, "at most 64 threads, or 0 for automatic selection")
            graph = Frontier(side)
            print(game_count(side * side, crossing_counts(graph, threads)))
    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.util.ArrayList;
import java.util.Arrays;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
import java.util.function.LongToIntFunction;

public class Euler1018 {
    // Residues below MOD live in an int, as the u32 values do in the C++ program: a sum of two
    // stays below 2^31 because 2 * (MOD - 1) < 2^31, and every product goes through a long.
    private static final int MOD = 1_000_000_007;
    private static final int MAX_SIDE = 11;
    private static final int NONE = -1;
    private static final long UNSIGNED_MAX = 4_294_967_295L;

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

    static int add(int a, int b) {
        int sum = a + b;
        return sum >= MOD ? sum - MOD : sum;
    }

    static int multiply(int a, int b) {
        return (int) ((long) a * b % MOD);
    }

    static int power(int base, int exponent) {
        int result = 1;
        for (; exponent != 0; exponent >>= 1, base = multiply(base, base)) {
            if ((exponent & 1) != 0) result = multiply(result, base);
        }
        return result;
    }

    static long advance(long state, int column, int side, boolean occupied) {
        int[] labels = new int[MAX_SIDE];
        for (int i = 0; i < side; ++i) labels[i] = (int) ((state >>> (4 * i)) & 15);
        if (!occupied) {
            labels[column] = 0;
        } else {
            int[] adjacent = {
                column == 0 ? 0 : labels[column - 1], labels[column],
                column + 1 == side ? 0 : labels[column + 1]
            };
            int merged = 15;
            for (int label : adjacent) {
                if (label != 0) merged = Math.min(merged, label);
            }
            for (int i = 0; i < side; ++i) {
                for (int label : adjacent) {
                    if (label != 0 && labels[i] == label) labels[i] = merged;
                }
            }
            labels[column] = merged;
        }

        // Label 1 connects to the top edge; other component labels are canonicalized.
        int[] canonical = new int[16];
        canonical[1] = 1;
        int next = 2;
        boolean reachesTop = false;
        long result = 0;
        for (int i = 0; i < side; ++i) {
            reachesTop |= labels[i] == 1;
            if (labels[i] != 0 && canonical[labels[i]] == 0) canonical[labels[i]] = next++;
            result |= (long) canonical[labels[i]] << (4 * i);
        }
        return reachesTop ? result : 0;
    }

    private static final class Incoming {
        final int[] offsets;
        final int[] sources;

        Incoming(int[] offsets, int[] sources) {
            this.offsets = offsets;
            this.sources = sources;
        }
    }

    private static final class Frontier {
        final int side;
        final long[] states;
        final int[][] initial; // {state, number of stones in the first row}
        final int[][] edges;
        final Incoming[] incoming;

        // The width is the parsed unsigned value, checked before it is narrowed to an int.
        Frontier(long width) {
            require(width >= 1 && width <= MAX_SIDE, "board side between 1 and 11");
            side = (int) width;
            incoming = new Incoming[side];
            Map<Long, Integer> ids = new HashMap<>(1 << (side + 4));
            List<Long> found = new ArrayList<>();
            LongToIntFunction intern = state -> {
                Integer id = ids.putIfAbsent(state, found.size());
                if (id != null) return id;
                found.add(state);
                return found.size() - 1;
            };
            initial = new int[(1 << side) - 1][];
            for (int mask = 1; mask < (1 << side); ++mask) {
                long state = 0;
                int count = 0;
                for (int column = 0; column < side; ++column) {
                    if ((mask & (1 << column)) != 0) {
                        state |= 1L << (4 * column);
                        ++count;
                    }
                }
                initial[mask - 1] = new int[] {intern.applyAsInt(state), count};
            }
            List<int[]> edgeList = new ArrayList<>();
            for (int i = 0; i < found.size(); ++i) {
                long state = found.get(i);
                int[] transitions = new int[2 * MAX_SIDE];
                for (int column = 0; column < side; ++column) {
                    for (int occupied = 0; occupied < 2; ++occupied) {
                        long next = advance(state, column, side, occupied != 0);
                        transitions[2 * column + occupied] = next == 0 ? NONE : intern.applyAsInt(next);
                    }
                }
                edgeList.add(transitions);
            }
            states = found.stream().mapToLong(Long::longValue).toArray();
            edges = edgeList.toArray(new int[0][]);
            for (int column = 0; column < side; ++column) {
                int[] offsets = new int[states.length + 1];
                for (int[] transitions : edges) {
                    for (int occupied = 0; occupied < 2; ++occupied) {
                        int target = transitions[2 * column + occupied];
                        if (target != NONE) ++offsets[target + 1];
                    }
                }
                for (int i = 0; i < states.length; ++i) offsets[i + 1] += offsets[i];
                int[] sources = new int[offsets[states.length]];
                int[] positions = offsets.clone();
                for (int source = 0; source < states.length; ++source) {
                    for (int occupied = 0; occupied < 2; ++occupied) {
                        int target = edges[source][2 * column + occupied];
                        if (target != NONE) sources[positions[target]++] = source;
                    }
                }
                incoming[column] = new Incoming(offsets, sources);
            }
        }
    }

    // The C++ barrier: one lock, one condition and a generation count. A participant that fails
    // also breaks it, so the others stop instead of waiting for an arrival that never comes.
    private static final class Barrier {
        private int arrived = 0;
        private int generation = 0;
        private Throwable failure = null;

        /** Waits for every participant; false once another participant has failed. */
        boolean await(int participants) throws InterruptedException {
            if (participants == 1) return true;
            synchronized (this) {
                if (failure != null) return false;
                int previous = generation;
                if (++arrived == participants) {
                    arrived = 0;
                    ++generation;
                    notifyAll();
                } else {
                    while (generation == previous && failure == null) wait();
                }
                return failure == null;
            }
        }

        /** Keeps the first failure and wakes every waiting participant. */
        synchronized void fail(Throwable error) {
            if (failure == null) failure = error;
            notifyAll();
        }

        synchronized Throwable failure() {
            return failure;
        }
    }

    private static final class Work {
        final Frontier frontier;
        final int[] values;
        final int[] scratch;
        final int stride;
        int threads; // final count, written under the barrier's lock before any worker reads it
        final Barrier barrier = new Barrier();

        Work(Frontier graph, int[] first, int[] second, int length, int count) {
            frontier = graph;
            values = first;
            scratch = second;
            stride = length;
            threads = count;
        }
    }

    private static final class Task {
        final Work work;
        final int index;

        Task(Work work, int index) {
            this.work = work;
            this.index = index;
        }
    }

    static void countWorker(Task task) {
        Work work = task.work;
        try {
            int threads;
            synchronized (work.barrier) {
                threads = work.threads; // the caller holds this lock until every worker has started
            }
            Frontier graph = work.frontier;
            int first = graph.states.length * task.index / threads;
            int last = graph.states.length * (task.index + 1) / threads;
            int[] current = work.values;
            int[] next = work.scratch;
            for (int cell = graph.side; cell < graph.side * graph.side; ++cell) {
                int column = cell % graph.side;
                int upper = Math.min(work.stride - 1, cell + 1);
                Incoming incoming = graph.incoming[column];
                for (int target = first; target < last; ++target) {
                    int occupied = ((graph.states[target] >>> (4 * column)) & 15) != 0 ? 1 : 0;
                    int output = target * work.stride;
                    Arrays.fill(next, output, output + upper + 1, 0);
                    for (int j = incoming.offsets[target]; j < incoming.offsets[target + 1]; ++j) {
                        int input = incoming.sources[j] * work.stride;
                        for (int degree = occupied; degree <= upper; ++degree) {
                            next[output + degree] = add(next[output + degree],
                                    current[input + degree - occupied]);
                        }
                    }
                }
                if (!work.barrier.await(threads)) return;
                int[] swap = current;
                current = next;
                next = swap;
            }
        } catch (Throwable error) {
            // The caller re-throws the first failure once every thread has stopped.
            work.barrier.fail(error);
        }
    }

    static int[] crossingCounts(Frontier graph, int threads) {
        int stride = graph.side * graph.side / 2 + 1;
        int[] values = new int[graph.states.length * stride], scratch = new int[values.length];
        for (int[] entry : graph.initial) {
            if (entry[1] < stride) values[entry[0] * stride + entry[1]] = 1;
        }
        if (threads == 0) {
            threads = graph.side < 9 ? 1
                    : Math.min(8, Math.max(1, Runtime.getRuntime().availableProcessors()));
        }
        threads = Math.min(threads, graph.states.length);
        Work work = new Work(graph, values, scratch, stride, threads);
        Task[] tasks = new Task[threads];
        for (int i = 0; i < threads; ++i) tasks[i] = new Task(work, i);
        Thread[] workers = new Thread[threads - 1];
        int started = 1;
        synchronized (work.barrier) {
            for (; started < threads; ++started) {
                Task task = tasks[started];
                try {
                    workers[started - 1] = new Thread(() -> countWorker(task));
                    workers[started - 1].start();
                } catch (Throwable notStarted) {
                    break; // like a failed pthread_create: the threads already running share the work
                }
            }
            work.threads = started;
        }
        countWorker(tasks[0]);
        boolean joined = true;
        for (int i = 1; i < started; ++i) {
            for (;;) {
                try {
                    workers[i - 1].join();
                    break;
                } catch (InterruptedException interrupted) {
                    joined = false; // keep waiting: the worker still writes into the shared arrays
                }
            }
        }
        if (!joined) Thread.currentThread().interrupt();
        Throwable failure = work.barrier.failure();
        if (failure instanceof RuntimeException runtime) throw runtime;
        if (failure instanceof Error error) throw error;
        if (failure != null) throw new IllegalStateException(failure);
        require(joined, "thread join");
        int[] result = new int[stride];
        for (int state = 0; state < graph.states.length; ++state) {
            for (int degree = 0; degree < stride; ++degree) {
                result[degree] = add(result[degree], values[state * stride + degree]);
            }
        }
        return result;
    }

    static int gameCount(int cells, int[] crossings) {
        int[] factorial = new int[cells + 1], inverse = new int[cells + 1];
        int[] probability = new int[cells / 2 + 2];
        Arrays.fill(factorial, 1);
        Arrays.fill(inverse, 1);
        for (int i = 1; i <= cells; ++i) factorial[i] = multiply(factorial[i - 1], i);
        inverse[cells] = power(factorial[cells], MOD - 2);
        for (int i = cells; i != 0; --i) inverse[i - 1] = multiply(inverse[i], i);
        for (int k = 0; k <= cells / 2; ++k) {
            probability[k] = multiply(crossings[k], multiply(inverse[cells],
                    multiply(factorial[k], factorial[cells - k])));
        }
        if (cells % 2 != 0) probability[cells / 2 + 1] = (1 + MOD - probability[cells / 2]) % MOD;

        // With p_k the crossing probability, wins on moves 2k+1 and 2k+2 total
        // (cells)_(2k+1) * (cells - 2k) * (p_(k+1) - p_k).
        int result = 0;
        for (int k = 0; k < (cells + 1) / 2; ++k) {
            int difference = (probability[k + 1] + MOD - probability[k]) % MOD;
            int weight = multiply(multiply(factorial[cells], inverse[cells - 2 * k - 1]), cells - 2 * k);
            result = add(result, multiply(weight, difference));
        }
        return result;
    }

    private static final class SmallBoard {
        private static final int[] DR = {-1, -1, 0, 0, 1, 1};
        private static final int[] DC = {0, 1, -1, 1, -1, 0};

        final int side;
        final int cells;
        final int[] adjacent;
        final int[] first = new int[2], last = new int[2];

        SmallBoard(int n) {
            side = n;
            cells = n * n;
            adjacent = new int[cells];
            for (int cell = 0; cell < cells; ++cell) {
                int row = cell / side, column = cell % side;
                int bit = 1 << cell;
                if (row == 0) first[0] |= bit;
                if (row + 1 == side) last[0] |= bit;
                if (column == 0) first[1] |= bit;
                if (column + 1 == side) last[1] |= bit;
                for (int d = 0; d < DR.length; ++d) {
                    int r = row + DR[d], c = column + DC[d];
                    if (r >= 0 && c >= 0 && r < side && c < side) adjacent[cell] |= 1 << (r * side + c);
                }
            }
        }

        boolean wins(int stones, int player) {
            int connected = stones & first[player];
            for (;;) {
                if ((connected & last[player]) != 0) return true;
                int next = connected;
                for (int cell = 0; cell < cells; ++cell) {
                    if ((connected & (1 << cell)) != 0) next |= adjacent[cell] & stones;
                }
                if (next == connected) return false;
                connected = next;
            }
        }

        long games(int white, int black, int player, long[] memo) {
            int key = white | (black << cells);
            if (memo[key] != -1) return memo[key];
            if (wins(white, 0) || wins(black, 1)) return memo[key] = 1;
            long result = 0;
            for (int cell = 0; cell < cells; ++cell) {
                int bit = 1 << cell;
                if (((white | black) & bit) != 0) continue;
                result += player == 0 ? games(white | bit, black, 1, memo)
                        : games(white, black | bit, 0, memo);
            }
            return memo[key] = result;
        }
    }

    static int evaluateCrossings(Frontier graph, int weight) {
        int[] current = new int[graph.states.length], next = new int[current.length];
        for (int[] entry : graph.initial) current[entry[0]] = power(weight, entry[1]);
        for (int cell = graph.side; cell < graph.side * graph.side; ++cell) {
            Arrays.fill(next, 0);
            int column = cell % graph.side;
            for (int source = 0; source < current.length; ++source) {
                for (int occupied = 0; occupied < 2; ++occupied) {
                    int target = graph.edges[source][2 * column + occupied];
                    if (target == NONE) continue;
                    next[target] = add(next[target],
                            occupied == 0 ? current[source] : multiply(current[source], weight));
                }
            }
            int[] swap = current;
            current = next;
            next = swap;
        }
        int result = 0;
        for (int value : current) result = add(result, value);
        return result;
    }

    static void runTests() {
        final long[] given = {1, 24, 257_760, 10_663_950_528_000L};
        for (int side = 1; side <= 4; ++side) {
            Frontier graph = new Frontier(side);
            int[] actual = crossingCounts(graph, 1);
            SmallBoard board = new SmallBoard(side);
            int[] expected = new int[board.cells + 1];
            int all = (1 << board.cells) - 1;
            for (int stones = 0; stones <= all; ++stones) {
                boolean whiteWins = board.wins(stones, 0);
                require(whiteWins != board.wins(all ^ stones, 1), "exactly one winner on a full board");
                if (whiteWins) ++expected[Integer.bitCount(stones)];
            }
            require(Arrays.equals(actual, 0, actual.length, expected, 0, actual.length),
                    "all subsets through side 4");
            int answer = gameCount(board.cells, actual);
            require(answer == given[side - 1] % MOD, "given game count");
            if (side <= 3) {
                long[] memo = new long[1 << (2 * board.cells)];
                Arrays.fill(memo, -1);
                require(answer == board.games(0, 0, 0, memo) % MOD, "direct game tree through side 3");
            }
        }

        Frontier graph = new Frontier(MAX_SIDE);
        int[] serial = crossingCounts(graph, 1);
        require(Arrays.equals(serial, crossingCounts(graph, 4)), "serial and thread agreement");
        for (int k = 0; k < MAX_SIDE; ++k) require(serial[k] == 0, "at least one stone in every row");
        require(serial[MAX_SIDE] == (MAX_SIDE + 1) * (1 << (MAX_SIDE - 2)), "minimal crossing paths");

        int cells = MAX_SIDE * MAX_SIDE;
        int[] binomial = new int[cells + 1], full = new int[cells + 1];
        binomial[0] = 1;
        for (int row = 1; row <= cells; ++row) {
            for (int k = row; k != 0; --k) binomial[k] = add(binomial[k], binomial[k - 1]);
        }
        for (int k = 0; k <= cells / 2; ++k) {
            full[k] = serial[k];
            full[cells - k] = (binomial[k] + MOD - serial[k]) % MOD;
        }
        for (int weight : new int[] {1, 2}) {
            int value = 0;
            for (int k = cells + 1; k != 0; --k) value = add(multiply(value, weight), full[k - 1]);
            require(value == evaluateCrossings(graph, weight),
                    "full polynomial and weighted forward transfer");
            if (weight == 1) {
                require(value == power(2, cells - 1), "half of all full boards connect each side");
            }
        }
        System.out.println("All checks passed; supplied examples, exhaustive subsets through side 4, "
                + "game trees through side 3, thread agreement, and weighted transfers at side 11.");
    }

    // Like std::from_chars into a 32-bit unsigned value: ASCII digits only, no sign or spaces.
    static long parseUnsigned(String text) {
        long value = 0;
        boolean digits = !text.isEmpty();
        for (int i = 0; i < text.length() && digits && value <= UNSIGNED_MAX; ++i) {
            char digit = text.charAt(i);
            digits = digit >= '0' && digit <= '9';
            value = 10 * value + (digit - '0');
        }
        if (!digits || value > UNSIGNED_MAX) {
            throw new IllegalArgumentException("Expected an unsigned integer");
        }
        return value;
    }

    public static void main(String[] args) {
        try {
            if (args.length == 1 && args[0].equals("--self-test")) {
                runTests();
            } else {
                if (args.length > 2) {
                    throw new IllegalArgumentException(
                            "Usage: java Euler1018 [side [threads]] | --self-test");
                }
                long side = args.length > 0 ? parseUnsigned(args[0]) : MAX_SIDE;
                long threads = args.length > 1 ? parseUnsigned(args[1]) : 0;
                require(threads <= 64, "at most 64 threads, or 0 for automatic selection");
                Frontier graph = new Frontier(side);
                System.out.println(gameCount(graph.side * graph.side, crossingCounts(graph, (int) threads)));
            }
        } catch (RuntimeException error) {
            System.err.println(error.getMessage());
            System.exit(1);
        }
    }
}