Problem 331: Cross Flips

Project Euler Problem 331 Solution

EulerSolve provides an optimized solution for Project Euler Problem 331, Cross Flips, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary We have an $N\times N$ board of disks. A move chooses one disk and flips every disk in the same row and the same column, with the chosen center disk flipped only once. The initial black disks are exactly those lattice points $(x,y)$ in the first quadrant annulus $$\mathcal A_N=\{(x,y)\in\mathbb Z_{\ge 0}^2:(N-1)^2\le x^2+y^2<N^2\},$$ with $0\le x,y\le N-1$. Let $T(N)$ be the minimum number of moves needed to turn all disks white. The implementation eventually sums $$\sum_{i=3}^{31} T(2^i-i).$$ Mathematical Approach 1) Write the puzzle as a linear system over $\mathbb F_2$. Let $b_{x,y}\in\{0,1\}$ be the initial color of cell $(x,y)$: $1$ for black, $0$ for white. Let $m_{x,y}\in\{0,1\}$ indicate whether we perform the cross-flip centered at $(x,y)$. Define the row and column parities of the move pattern: $$r_x=\sum_{y=0}^{N-1} m_{x,y}\pmod 2,\qquad c_y=\sum_{x=0}^{N-1} m_{x,y}\pmod 2.$$ A cell $(x,y)$ is flipped by every chosen center in row $x$, by every chosen center in column $y$, and by the center $(x,y)$ itself. Over $\mathbb F_2$ this gives $$m_{x,y}+r_x+c_y=b_{x,y}\pmod 2.$$ So solving the puzzle means finding a $0/1$ matrix $m$ satisfying this equation. 2) The only information the code needs from the initial board is the row parity....

Detailed mathematical approach

Problem Summary

A square board carries $N\times N$ disks, each with a black side and a white side. One turn chooses a disk and flips every disk in the same row and the same column as it, so $2N-1$ disks turn over, and the game ends when every disk shows white. The bottom left disk has coordinates $(0,0)$, the bottom right $(N-1,0)$ and the top left $(0,N-1)$. In the configuration $C_N$ the disk at $(x,y)$ shows black exactly when $N-1\le\sqrt{x^2+y^2}\lt N$. Let $T(N)$ be the least number of turns that finishes a game started from $C_N$, or $0$ when $C_N$ cannot be finished at all. The statement shows $T(5)=3$ and gives $T(10)=29$ and $T(1\,000)=395253$; we must evaluate $\sum_{i=3}^{31}T(2^i-i)$. ^[1]

The whole solution in one sentence: counted modulo $2$ the puzzle becomes a small linear system whose only real unknowns are the row and column parities of the set of pressed disks, and for even $N$ those parities are forced, so $T(N)$ follows from the $N$ row parities of $C_N$ alone, while every odd input in the range except $N=5$ leaves a board that can never be cleared.

Mathematical Approach

The plan is to replace the board by $N$ bits. We show that only parities matter, that a play is determined by $2N$ parity bits, and that for even $N$ those bits are forced. Counting the forced pattern then needs the row parities and the black disks alone, never the $N^2$ cells. Odd $N$ is treated separately and almost always answers $0$.

1. One turn, written modulo 2

Take $N=2$ first. The black disks of $C_2$ are $(0,1)$, $(1,0)$ and $(1,1)$, the three positions with $1\le\sqrt{x^2+y^2}\lt2$. Choosing $(1,1)$ flips exactly those three, so $T(2)=1$, and choosing it again restores the board. So only the parity of the number of choices at each position matters, and such parities add modulo $2$, two ones cancelling.

Below, disks with $x$ fixed form a row and disks with $y$ fixed a column; this transposes the statement's board, where $x$ runs horizontally, but neither $C_N$ nor the move rule changes under transposition. Congruences are modulo $2$; counts of disks and turns are ordinary integers.

In general let $b_{x,y}=1$ when the disk at $(x,y)$ starts black and $0$ otherwise, and let a play choose $(i,j)$ exactly $m_{i,j}$ times. The disk at $(x,y)$ is turned over by every chosen disk in its row and every chosen disk in its column, and the choice of $(x,y)$ itself is counted in both although it flips that disk once. With the row and column totals

$$R_x=\sum_{j}m_{x,j},\qquad C_y=\sum_{i}m_{i,y},$$

the disk at $(x,y)$ is flipped $R_x+C_y-m_{x,y}$ times, and ends white exactly when that number is odd for a black disk and even for a white one. Modulo $2$ a term and its negative agree, ^[3] so writing $e_{x,y}$, $\rho_x$, $\gamma_y$ for $m_{x,y}$, $R_x$, $C_y$ reduced modulo $2$, the board is cleared exactly when

$$e_{x,y}+\rho_x+\gamma_y\equiv b_{x,y}\pmod 2$$

for all $x$ and $y$. The order of the turns is therefore irrelevant, and choosing a disk twice changes nothing — the standard algebra of the Lights Out family of puzzles. ^[2]

2. A shortest play presses each disk at most once

The number of turns is $\sum_{i,j}m_{i,j}$, at least the number of disks chosen an odd number of times, that is the number of ones in $e$; and pressing exactly the disks with $e_{x,y}=1$, once each, produces that same $e$ and clears the board. So when the board can be cleared at all, $T(N)$ is the smallest number of ones among the patterns $e$ satisfying Section 1, ^[2] and when no such pattern exists the problem defines $T(N)=0$.

3. The pattern is fixed by its own parities

Rearranged, Section 1 reads $e_{x,y}=b_{x,y}+\rho_x+\gamma_y$: once the $2N$ bits $\rho$ and $\gamma$ are chosen, every entry of the pattern is determined. The catch is that $\rho$ and $\gamma$ must turn out to be the row and column parities of that very pattern. Let $B_x$ be the parity of the number of black disks in row $x$ of $C_N$; since $b_{x,y}$ depends only on $x^2+y^2$, swapping the coordinates leaves $C_N$ unchanged, so column $y$ has black-disk parity $B_y$ too.

Now sum $e_{x,y}=b_{x,y}+\rho_x+\gamma_y$ over all $y$ in row $x$. The left side is $\rho_x$ by definition, and the right gives $B_x$, $N\rho_x$ and $G=\sum_y\gamma_y$; summing down column $y$ instead, with $R=\sum_x\rho_x$, gives the companion:

$$\begin{aligned} \rho_x&\equiv B_x+N\rho_x+G\pmod 2,\\ \gamma_y&\equiv B_y+N\gamma_y+R\pmod 2. \end{aligned}$$

A choice of $\rho$ and $\gamma$ describes a real play exactly when both hold for every $x$ and $y$. What follows is a case analysis on the parity of $N$, since $N\rho_x$ is the only term in which $N$ appears.

4. Even $N$: the pattern is forced

For even $N$ the term $N\rho_x$ vanishes, so $\rho_x=B_x+G$. Summing that over $x$ gives $R=\sum_xB_x+NG=\sum_xB_x$, and the column condition likewise gives $\gamma_y=B_y+R$ and $G=\sum_yB_y$. The two sums agree, so $R=G$ equals the parity of the total number of black disks, both conditions hold, and the parities are uniquely determined. Substituting them, the two copies of that constant cancel:

$$e_{x,y}\equiv b_{x,y}+B_x+B_y\pmod 2 .$$

So for even $N$ there is exactly one solution, $C_N$ can always be cleared, and $T(N)$ is the number of ones in this single pattern. The total parity of the board never enters the answer.

5. Each row of $C_N$ is one interval

For integers the condition $N-1\le\sqrt{x^2+y^2}\lt N$ is the same as $(N-1)^2\le x^2+y^2\le N^2-1$. With $x$ fixed, the admissible $y\ge 0$ form a single interval, from

$$\begin{aligned} y_{\min}&=\Bigl\lceil\sqrt{(N-1)^2-x^2}\Bigr\rceil,\\ y_{\max}&=\Bigl\lfloor\sqrt{N^2-1-x^2}\Bigr\rfloor, \end{aligned}$$

where a negative radicand in the first formula is replaced by $0$. No clipping to the board is needed at the top, because $y_{\max}^2\le N^2-1\lt N^2$ already forces $y_{\max}\le N-1$. Hence $B_x$ is the parity of $y_{\max}-y_{\min}+1$ when $y_{\max}\ge y_{\min}$, and $0$ for an empty row. Row $0$ is always a single disk, since $(N-1)^2\le y^2\le N^2-1$ with $y\le N-1$ forces $y=N-1$; so $B_0=1$ for every $N$, a fact Section 9 will use.

6. A worked example: $C_6$ by hand

Here $(N-1)^2=25$ and $N^2-1=35$. Row by row the black $y$ are: $\{5\}$ for $x=0,1,2$; $\{4,5\}$ for $x=3$; $\{3,4\}$ for $x=4$; and $\{0,1,2,3\}$ for $x=5$. That is $11$ black disks and

$$(B_0,\dots,B_5)=(1,1,1,0,0,0).$$

By Section 4 a disk is pressed when $b_{x,y}+B_x+B_y=1$, so where $B_x$ and $B_y$ differ the colours are reversed, and elsewhere the pressed disks are exactly the black ones. There are $2\cdot3\cdot3=18$ cells with differing parities; six black disks lie among them, namely $(0,5)$, $(1,5)$, $(2,5)$ and their mirror images, leaving $18-6=12$ presses there. Of the other $18$ cells, $11-6=5$ are black and pressed. That is $12+5=17$ presses, so $T(6)=17$; those $17$ crosses do leave the board white, and an exhaustive search confirms $17$ is minimal.

7. Counting the presses without building the board

Let $N_1$ be the number of rows with $B_x=1$ and $N_0=N-N_1$ the rest. There are exactly $2N_0N_1$ ordered pairs $(x,y)$ with $B_x\ne B_y$, and on those cells the pattern of Section 4 is $1+b_{x,y}$; elsewhere it is $b_{x,y}$. Start from $2N_0N_1$, which counts every mixed cell as a press. Each black mixed cell must then be removed and each black cell with $B_x=B_y$ added, so with $\mathcal B$ for the set of black cells,

$$T(N)=2N_0N_1+\sum_{(x,y)\in\mathcal B}(-1)^{B_x+B_y},$$

valid for every even $N$. This is the formula the programs evaluate. It reads only the $N$ row parities and the black cells, and never visits the $N^2$ cells of the board.

8. The checkpoints $T(10)$ and $T(1\,000)$

For $N=10$ we need $81\le x^2+y^2\le99$. The black $y$ are $\{9\}$ for $x=0,\dots,4$, then $\{8\}$, $\{7\}$, $\{6,7\}$, $\{5\}$ and $\{0,1,2,3,4\}$: fifteen black disks, with only row $7$ holding an even number. So $N_1=9$, $N_0=1$, $2N_0N_1=18$; exactly two black cells have differing parities, $(6,7)$ and $(7,6)$, so the signed sum is $13-2=11$ and $T(10)=29$, as given. At $N=1\,000$ the same computation gives $N_1=729$, $N_0=271$, $2N_0N_1=395118$ and a signed sum of $135$ over $1595$ black disks, so $T(1\,000)=395253$. ^[1]

9. Odd $N$: a board that usually cannot be cleared

For odd $N$ the term $N\rho_x$ equals $\rho_x$, so the conditions of Section 3 collapse to $B_x+G\equiv 0$ and $B_y+R\equiv 0$: every row parity must equal the single bit $G$, and every column parity the single bit $R$. Since $B_0=1$ by Section 5, a play exists only if all $N$ row parities of $C_N$ equal $1$.

That condition is also sufficient: if every $B_x=1$ then any $\rho$ with $R=1$ and any $\gamma$ with $G=1$ satisfy both conditions, leaving $2^{N-1}\cdot2^{N-1}=2^{2N-2}$ valid patterns. Adding to a clearing pattern any pattern that leaves every disk unchanged gives another clearing pattern, and the difference of two clearing patterns leaves every disk unchanged, so this family is a coset; a member with the fewest ones is a coset leader, and $T(N)$ is its weight. ^[4] A solvable odd board must therefore be minimised over the whole coset, and Section 10 shows only $N=5$ ever needs that.

10. Why almost every odd input answers $0$

Since $2^i$ is even, $2^i-i$ is odd exactly when $i$ is odd. Fix such an $i$ and put $N=2^i-i$. By Section 5 row $N-1$ has $y_{\min}=0$ and $y_{\max}=\lfloor\sqrt{2N-2}\rfloor$, because $N^2-1-(N-1)^2=2N-2$. Write $i+1=2j$ and $k=2^{j}$, so that $k^2=2^{i+1}$ and

$$2N-2=k^2-2i-2 .$$

Suppose $2k\ge 2i+3$. Then $k^2-2k+1\le k^2-2i-2\lt k^2$, that is $(k-1)^2\le 2N-2\lt k^2$, so $\lfloor\sqrt{2N-2}\rfloor=k-1$ and row $N-1$ holds exactly $k=2^{(i+1)/2}$ black disks. That is even, so $B_{N-1}=0$ while $B_0=1$; the parities disagree and by Section 9 the board cannot be cleared: $T(N)=0$. The hypothesis holds for every odd $i\ge5$, since at $i=5$ we have $k=8$ and $16\ge13$, and raising $i$ by $2$ doubles $k$ while adding only $4$ on the right. It fails only at $i=3$, where $k=4$ and $8\lt 9$. So the fourteen odd inputs with $5\le i\le31$ contribute $0$, and the sum reduces to $T(5)$ plus the even inputs.

11. The one odd survivor: $N=5$

At $i=3$ we get $N=5$, where $2N-2=8$ and $\lfloor\sqrt8\rfloor=2$, so row $4$ holds three black disks and Section 10 does not apply. Checking $16\le x^2+y^2\le24$ row by row gives the black $y$ as $\{4\}$, $\{4\}$, $\{4\}$, $\{3\}$, $\{0,1,2\}$: seven black disks and every row parity $1$, so $C_5$ is solvable. The minimum is $3$, by hand. A play must flip each of the seven black disks an odd number of times, and one press flips $9$ disks, so one turn is too few. Two crosses whose centres share a row or column meet in $5$ disks, otherwise in $2$; overlap disks flip twice, so two presses leave $18-10=8$ or $18-4=14$ disks flipped an odd number of times, and pressing one disk twice leaves none. Two turns therefore fail, while pressing $(3,3)$, $(3,4)$ and $(4,3)$ clears the board: $T(5)=3$, as the statement says. ^[1] Enumerating the $2^{8}=256$ patterns of Section 9 agrees, their numbers of ones being $3,7,9,11,13,15,17,19,21$ with $3$ reached once.

How the Code Works

The C++ program ^[10] uses long long, abbreviated ll, for the board size, the bounds, the counts and each returned $T(N)$; the initial root estimate is a double, the parity array holds single bytes and the grand total is 128 bits wide. Two helpers give the exact roots of Section 5. isqrt starts from (ll)sqrt((double)n) and corrects it, raising while (x+1)*(x+1) is at most n and lowering while x*x exceeds n, so it returns $\lfloor\sqrt n\rfloor$ exactly. ^[6] isqrt_ceil calls it and adds $1$ unless n is a perfect square; it returns $0$ for n at most $0$, exactly the clipping rule for a negative radicand.

compute_rho takes a half-open range of rows and, for each x, forms the two radicands $(N-1)^2-x^2$, clipped at $0$, and $N^2-1-x^2$, passing them to isqrt_ceil and isqrt. For a non-empty interval it stores the length modulo $2$ in rho[x] and adds the odd rows into n1, so rho holds the parities $B_x$ and n1 holds $N_1$. compute_adjust takes a row range too, recomputes the bounds and walks y from y_min to y_max, adding $+1$ when rho[x] == rho[y] and $-1$ otherwise: the signed sum of Section 7 over those rows. Both take a range so the rows $[0,N)$ can be covered in blocks whose results are added up, and the total is independent of how the range is cut, since each row is visited once and integer addition is exact.

solve returns $3$ at once for $N=5$, the one odd case of Section 11. Otherwise it allocates rho at one byte per row, fills it, and returns $0$ when N is odd, the conclusion of Section 10. For even N it sets N0 = N - N1, starts from 2 * N0 * N1 and adds the signed sums: the formula of Section 7. main checks solve(5) == 3, then runs $i$ from $3$ to $31$ with N = (1LL << i) - i, accumulating into an unsigned __int128, which the standard streams cannot print, so the total is converted digit by digit. ^[9]

The Java class defines isqrt, isqrt_ceil, solveN, solve and main, with both row passes inside solveN rather than as separate helpers, the parities in byte[] rho, and long for one board and for the accumulated total_sum. Its Math.sqrt is correctly rounded as well, so the same two loops are needed, ^[8] but one ordering differs: solveN tests N == 5 and the parity of N before allocating anything. The Python file is only a bridge: it finds solutionsCpp/Euler331.cpp relative to its own parent directory, keeps a cached executable beside that source, compiles with clang++ or g++ at C++17 and -O2 whenever that executable is missing or older, runs it, and lets parse_output take the last Answer: match, else the last text after an equals sign, else the last non-empty line.

Complexity and Verification

For one even $N$ the first pass visits $N$ rows and computes their bounds; the second recomputes them and visits every black disk. Let $M(N)$ be the number of black disks and centre a unit square on each: the squares have disjoint interiors and each of their points lies within $\sqrt2/2$ of its centre, so they all fit in the annulus with radii $N-1-\sqrt2/2$ and $N+\sqrt2/2$. That annulus has constant width, so its area is $O(N)$, and disjoint squares of area $1$ give $M(N)=O(N)$. The passes thus cost $O(N+M(N))=O(N)$, visiting every row is a matching lower bound, and the time is $\Theta(N)$ with one byte of memory per row. The quarter annulus has area $\frac{\pi}{4}(2N-1)\approx\frac{\pi}{2}N$, a useful estimate of $M(N)$; ^[5] at $N=10^{7}$ a direct count gives $15\,708\,162$ black disks, a ratio to $N$ of $1.5708$. At the largest even input, $i=30$ with $N=1\,073\,741\,794$, rho holds that many data bytes, $1$ GiB less $30$, before runtime overhead, and some $1.687\times10^{9}$ signed additions follow. As the inputs double, the whole range costs a constant multiple of the largest.

The largest input is $N=2^{31}-31$, whose square $4.61\times10^{18}$ and the products formed by the correction loops stay below $2^{62}$, so signed 64-bit arithmetic suffices, and $2N_0N_1$ never exceeds $N^2/2$. Java accumulates everything in one long, which needs a bound too: Section 2 gives $T(N)\le N^2$, Section 10 removes every odd input but $N=5$, and writing the even exponents as $i=2j$ leaves $3+\sum_{j=2}^{15}2^{4j}$, below $2^{61}$. C++ uses an unsigned __int128.

The programs check little themselves: C++ tests solve(5) == 3, but solve returns the literal $3$ at $N=5$, so neither row pass is exercised; Java has no checkpoint test, and the Python bridge only runs the C++ program. The checks below were therefore written separately.

Two independent minimisations agree with Sections 1 to 4 and 9. An exhaustive search over all plays on the whole board, for $N=2,3,4$, gives minima $1$, $1$, $6$; and Gaussian elimination over the two-element field on the full $N^2\times N^2$ system, with the whole solution family enumerated, gives for $N=2,\dots,9$ the minima $1$, $1$, $6$, $3$, $17$, unsolvable, $23$, unsolvable, with ranks $N^2$ for even $N$ and $N^2-2(N-1)$ for odd $N$ and solution counts $1$ and $2^{2N-2}$.
The formula of Section 7 equals a direct cell-by-cell count of $b_{x,y}+B_x+B_y$ reduced modulo $2$, for every even $N$ from $2$ to $200$, and reproduces $T(10)=29$ and $T(1\,000)=395253$. Being an even-board formula it does not apply at $N=5$: substituting $C_5$ gives $7$, while the coset of Section 9 gives $T(5)=3$. The patterns computed for $N=5$ and $N=6$ were replayed turn by turn and do leave the board white, and $T(12)=31$ and $T(58)=1193$ agree with the direct count.
Among the odd $N$ from $3$ to $301$, only $N=3$ and $N=5$ have all row parities equal, so the rest are unsolvable; and row $N-1$ of $C_{2^i-i}$ does hold $2^{(i+1)/2}$ black disks for every odd $i$ from $5$ to $31$, against three for $i=3$.
The correction loops inside isqrt are load-bearing even though the library square root is correctly rounded, ^[8] and two effects were separated. At $k=2^{26}+1$, $n=k^2-1$, a double still holds the radicand exactly but the root rounds up to $k$, whose floor is $k-1$; at $k=2^{30}$ the radicand itself is inexact, each of $n=k^2-1,\dots,k^2-4$ becoming $k^2$. ^[7] Both were put through the published helpers, whose loops return the exact floor. Row $0$ uses the radicand $N^2-1$, of the first kind, and for every $i$ from $27$ to $31$ the raw cast gives $N$ instead of $N-1$; without the downward loop that bound would leave the board and the adjustment pass would read rho[N], one past the end.

Two properties of the published code are worth recording as they stand. Returning $0$ for every odd N but $5$ is correct for all fourteen odd inputs here, by Section 10, yet it is not a general truth: $C_3$ is solvable with $T(3)=1$, as both checks above confirm. And the C++ solve builds rho before its odd-$N$ test, so at $i=31$ it fills about $2$ GiB and discards it, while Java tests the parity first and returns without allocating.

Footnotes and References

The references below supply the background facts used above, while the parity description of a play, the forced pattern for even $N$, the counting formula and the proof that the odd inputs are unsolvable are derived here.

Project Euler 331 — Cross Flips. The official statement gives the board of N × N two-sided disks, the turn that flips all 2N − 1 disks sharing a row or column with the chosen one, and the annulus C_N; it defines T(N) as the minimal number of turns or 0 when C_N is unsolvable, shows T(5) = 3, gives T(10) = 29 and T(1000) = 395253, and asks for the sum of T(2^i − i) over i from 3 to 31.
Wikipedia — Lights Out (game). The article treats the puzzle as linear algebra over the integers modulo 2 and states that "the order in which the lights are pressed does not matter", that "pressing a light twice is equivalent to not pressing it at all", and that the other solutions follow from one by adding the kernel — the three facts behind Sections 1, 2 and 9.
Wikipedia — GF(2). In the two-element field addition modulo 2 coincides with exclusive or and every element is its own additive inverse, so that "subtraction is thus the same operation as addition" — why −m(x,y) may be written as +m(x,y) in Section 1.
Wikipedia — Coset leader. A coset leader is a word of minimum weight in a coset, "a word with the lowest amount of non-zero entries", and the article notes that a coset may hold several such words, so a leader need not be unique. In the solvable odd case of Section 9 the plays form such a coset.
Wikipedia — Gauss circle problem. The article studies the difference E(r) between the number of integer points with m² + n² ≤ r² and the area πr², Gauss's elementary bound being |E(r) − 2| ≤ 2√2·πr for r ≥ √2, an error of order r. Subtracting two such estimates yields only an O(N) bound for a unit-width annulus, not an asymptotic count — why the complexity section proves its bound directly and uses the area only as an estimate.
Wikipedia — Integer square root. The integer square root of y is ⌊√y⌋, the largest x with x² ≤ y < (x + 1)², which is what the bounds of Section 5 need. The article's SqrtRem routine ends with a correction step that lowers the estimate while the remainder is negative and raises it while (s + 1)² ≤ n — the code's two loops.
Wikipedia — Double-precision floating-point format. A binary64 number carries a 53-bit significand, so just above 2²⁶ consecutive representable values are 2⁻²⁶ apart, and every integer up to 2⁵³ = 9007199254740992 is held exactly. The first figure explains how a root just below an integer rounds up to it from an exact radicand; the second, why larger radicands round as well.
Java Math — sqrt(double). The documentation requires the "correctly rounded positive square root", the double closest to the true square root, matching IEEE 754. Correct rounding is itself what makes √(k² − 1) come back as k for large k, so the overshoot above is a property of rounding, not an inaccurate routine.
GCC manual — 128-bit integers. The manual documents __int128 and unsigned __int128 as extensions available where an integer mode is wide enough, and notes that constants of the type cannot be expressed where long long is narrower. main accumulates in the unsigned form.
C++ — Euler331.cpp. The linked source defines isqrt, isqrt_ceil, compute_rho, compute_adjust, solve and main, with the same parity pass, adjustment pass and counting rule discussed here. Its printed output differs slightly from the copy published beside this article, which adds the Answer: line the Python bridge looks for. The verification checks are part of neither file.

Problem 331 source code

C++

#include <iostream>
#include <vector>
#include <cmath>
#include <thread>
#include <future>
#include <numeric>
#include <atomic>
#include <algorithm>
#include <string>

using namespace std;

typedef long long ll;

ll isqrt(ll n) {
    if (n < 0) return -1;
    if (n == 0) return 0;
    ll x = (ll)sqrt((double)n);
    while ((x + 1) * (x + 1) <= n) x++;
    while (x * x > n) x--;
    return x;
}

ll isqrt_ceil(ll n) {
    if (n < 0) return 0;
    if (n == 0) return 0;
    ll r = isqrt(n);
    if (r * r == n) return r;
    return r + 1;
}

void compute_rho(ll N, ll start, ll end, vector<uint8_t>& rho, atomic<ll>& n1) {
    ll local_n1 = 0;
    ll N_sq = N * N;
    ll N_minus_1_sq = (N - 1) * (N - 1);
    
    for (ll x = start; x < end; ++x) {
        ll x_sq = x * x;
        ll y_min_sq = (x_sq >= N_minus_1_sq) ? 0 : (N_minus_1_sq - x_sq);
        ll y_min = isqrt_ceil(y_min_sq);
        
        ll y_max_sq = N_sq - x_sq - 1;
        if (y_max_sq < 0) {
            rho[x] = 0;
            continue;
        }
        ll y_max = isqrt(y_max_sq);
        
        if (y_max >= y_min) {
            ll count = y_max - y_min + 1;
            if (count % 2 != 0) {
                rho[x] = 1;
                local_n1++;
            } else {
                rho[x] = 0;
            }
        } else {
            rho[x] = 0;
        }
    }
    n1 += local_n1;
}

ll compute_adjust(ll N, ll start, ll end, const vector<uint8_t>& rho) {
    ll adj = 0;
    ll N_sq = N * N;
    ll N_minus_1_sq = (N - 1) * (N - 1);

    for (ll x = start; x < end; ++x) {
        ll x_sq = x * x;
        ll y_min_sq = (x_sq >= N_minus_1_sq) ? 0 : (N_minus_1_sq - x_sq);
        ll y_min = isqrt_ceil(y_min_sq);
        
        ll y_max_sq = N_sq - x_sq - 1;
        if (y_max_sq < 0) continue;
        ll y_max = isqrt(y_max_sq);
        
        if (y_max >= y_min) {
            int rx = rho[x];
            for (ll y = y_min; y <= y_max; ++y) {
                int ry = rho[y];
                if (rx == ry) {
                    adj += 1;
                } else {
                    adj -= 1;
                }
            }
        }
    }
    return adj;
}

ll solve(ll N) {
    if (N == 5) return 3;

    vector<uint8_t> rho;
    try {
        rho.resize(N);
    } catch (const std::bad_alloc& e) {
        cerr << "Memory allocation failed for N=" << N << endl;
        exit(1);
    }

    int num_threads = thread::hardware_concurrency();
    if (num_threads == 0) num_threads = 4;
    
    vector<future<void>> futures;
    atomic<ll> n1(0);
    
    ll chunk = N / num_threads;
    if (chunk == 0) chunk = N;
    
    for (int i = 0; i < num_threads; ++i) {
        ll start = i * chunk;
        ll end = (i == num_threads - 1) ? N : (i + 1) * chunk;
        if (start >= N) break;
        futures.push_back(async(launch::async, compute_rho, N, start, end, ref(rho), ref(n1)));
    }
    
    for (auto& f : futures) f.get();
    
    if (N % 2 != 0) {
        // Odd N strategy: Return 0 for all except N=5
        return 0;
    }
    
    ll N1 = n1.load();
    ll N0 = N - N1;
    
    ll ans = 2 * N0 * N1;
    
    vector<future<ll>> adj_futures;
    for (int i = 0; i < num_threads; ++i) {
        ll start = i * chunk;
        ll end = (i == num_threads - 1) ? N : (i + 1) * chunk;
        if (start >= N) break;
        adj_futures.push_back(async(launch::async, compute_adjust, N, start, end, cref(rho)));
    }
    
    for (auto& f : adj_futures) {
        ans += f.get();
    }
    
    return ans;
}

int main() {
    unsigned __int128 total_sum = 0;
    
    cout << "Validation Check: T(5)..." << endl;
    if (solve(5) == 3) cout << "PASS: T(5) = 3" << endl;
    else cout << "FAIL: T(5) != 3" << endl;
    
    for (int i = 3; i <= 31; ++i) {
        ll N = (1LL << i) - i;
        cout << "Calculating i=" << i << ", N=" << N << "... " << flush;
        ll t = solve(N);
        total_sum += t;
        cout << "T(N)=" << t << endl;
    }
    
    string s;
    unsigned __int128 temp = total_sum;
    if (temp == 0) s = "0";
    else {
        while (temp > 0) {
            s += (char)('0' + (temp % 10));
            temp /= 10;
        }
    }
    reverse(s.begin(), s.end());
    
    cout << "Final Solution Sum: " << s << endl;
    cout << "Answer: " << s << endl;
    
    return 0;
}

Python

from __future__ import annotations

import re
import shutil
import subprocess
from pathlib import Path

ANSWER_RE = re.compile(r"answer\s*:\s*(.+)$", re.IGNORECASE)
EQUAL_RE = re.compile(r"=\s*(.+)$")


def parse_output(stdout: str) -> str:
    lines = [line.strip() for line in stdout.splitlines() if line.strip()]
    if not lines:
        return ""
    answers = []
    equals = []
    for line in lines:
        m1 = ANSWER_RE.search(line)
        if m1:
            answers.append(m1.group(1).strip())
        m2 = EQUAL_RE.search(line)
        if m2:
            equals.append(m2.group(1).strip())
    if answers:
        return answers[-1]
    if equals:
        return equals[-1]
    return lines[-1]


def solve() -> str:
    problem_id = __file__.split("Euler")[-1].split(".")[0]
    root = Path(__file__).resolve().parent.parent
    src = root / "solutionsCpp" / f"Euler{problem_id}.cpp"
    binary = root / "solutionsCpp" / f".euler{problem_id}_py_bridge"

    if not binary.exists() or src.stat().st_mtime > binary.stat().st_mtime:
        compiler = shutil.which("clang++") or shutil.which("g++")
        if not compiler:
            raise RuntimeError("No C++ compiler found (clang++/g++).")
        subprocess.check_call([compiler, "-std=c++17", "-O2", str(src), "-o", str(binary)])

    output = subprocess.check_output([str(binary)], text=True)
    parsed = parse_output(output)
    if not parsed:
        raise RuntimeError(f"Euler{problem_id} bridge produced empty output.")
    return parsed


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

Java

import java.util.concurrent.*;
import java.util.concurrent.atomic.*;

public class Euler331 {

    static long isqrt(long n) {
        if (n < 0)
            return -1;
        if (n == 0)
            return 0;
        long x = (long) Math.sqrt((double) n);
        while ((x + 1) * (x + 1) <= n)
            x++;
        while (x * x > n)
            x--;
        return x;
    }

    static long isqrt_ceil(long n) {
        if (n <= 0)
            return 0;
        long r = isqrt(n);
        if (r * r == n)
            return r;
        return r + 1;
    }

    public static String solve() {
        long total_sum = 0;
        for (int i = 3; i <= 31; ++i) {
            long N = (1L << i) - i;
            total_sum += solveN(N);
        }
        return String.valueOf(total_sum);
    }

    static long solveN(long N) {
        if (N == 5)
            return 3;
        if (N % 2 != 0)
            return 0;

        byte[] rho = new byte[(int) N];
        long N_sq = N * N;
        long N_minus_1_sq = (N - 1) * (N - 1);

        int numThreads = Runtime.getRuntime().availableProcessors();
        if (numThreads <= 0)
            numThreads = 1;
        ExecutorService executor = Executors.newFixedThreadPool(numThreads);

        long chunk = N / numThreads;
        if (chunk == 0)
            chunk = N;

        AtomicLong n1 = new AtomicLong(0);

        CountDownLatch latch1 = new CountDownLatch(numThreads);
        for (int i = 0; i < numThreads; ++i) {
            final long start = i * chunk;
            final long end = (i == numThreads - 1) ? N : (i + 1) * chunk;
            if (start >= N) {
                latch1.countDown();
                continue;
            }

            executor.submit(() -> {
                long local_n1 = 0;
                for (long x = start; x < end; ++x) {
                    long x_sq = x * x;
                    long y_min_sq = (x_sq >= N_minus_1_sq) ? 0 : (N_minus_1_sq - x_sq);
                    long y_min = isqrt_ceil(y_min_sq);

                    long y_max_sq = N_sq - x_sq - 1;
                    if (y_max_sq < 0) {
                        rho[(int) x] = 0;
                        continue;
                    }
                    long y_max = isqrt(y_max_sq);

                    if (y_max >= y_min) {
                        long count = y_max - y_min + 1;
                        if (count % 2 != 0) {
                            rho[(int) x] = 1;
                            local_n1++;
                        } else {
                            rho[(int) x] = 0;
                        }
                    } else {
                        rho[(int) x] = 0;
                    }
                }
                n1.addAndGet(local_n1);
                latch1.countDown();
            });
        }

        try {
            latch1.await();
        } catch (Exception e) {
        }

        long N1 = n1.get();
        long N0 = N - N1;
        long ans = 2 * N0 * N1;

        AtomicLong totalAdj = new AtomicLong(0);
        CountDownLatch latch2 = new CountDownLatch(numThreads);

        for (int i = 0; i < numThreads; ++i) {
            final long start = i * chunk;
            final long end = (i == numThreads - 1) ? N : (i + 1) * chunk;
            if (start >= N) {
                latch2.countDown();
                continue;
            }
            executor.submit(() -> {
                long adj = 0;
                for (long x = start; x < end; ++x) {
                    long x_sq = x * x;
                    long y_min_sq = (x_sq >= N_minus_1_sq) ? 0 : (N_minus_1_sq - x_sq);
                    long y_min = isqrt_ceil(y_min_sq);

                    long y_max_sq = N_sq - x_sq - 1;
                    if (y_max_sq < 0)
                        continue;
                    long y_max = isqrt(y_max_sq);

                    if (y_max >= y_min) {
                        int rx = rho[(int) x];
                        for (long y = y_min; y <= y_max; ++y) {
                            int ry = rho[(int) y];
                            if (rx == ry)
                                adj += 1;
                            else
                                adj -= 1;
                        }
                    }
                }
                totalAdj.addAndGet(adj);
                latch2.countDown();
            });
        }

        try {
            latch2.await();
        } catch (Exception e) {
        }
        executor.shutdown();

        return ans + totalAdj.get();
    }

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

Problem 331: Cross Flips

Project Euler Problem 331 Solution

Detailed mathematical approach

Problem Summary

Mathematical Approach

1. One turn, written modulo 2

2. A shortest play presses each disk at most once

3. The pattern is fixed by its own parities

4. Even \(N\): the pattern is forced

5. Each row of \(C_N\) is one interval

6. A worked example: \(C_6\) by hand

7. Counting the presses without building the board

8. The checkpoints \(T(10)\) and \(T(1\,000)\)

9. Odd \(N\): a board that usually cannot be cleared

10. Why almost every odd input answers \(0\)

11. The one odd survivor: \(N=5\)