Problem 310: Nim Square

Detailed mathematical approach

Problem Summary

Alice and Bob play with three heaps of stones whose sizes form the ordered triple \((a,b,c)\). A move takes a positive square number of stones — \(1\), \(4\), \(9\), \(16\), and so on — from one heap, never more than it holds, leaving the other two alone. Play is normal: a player with no legal move loses, so whoever takes the last stone wins. Nim Square is thus ordinary three-heap Nim, restricted to square takes. ^[1]

A position is losing for the next player when the player to move loses against perfect play. We must count the triples with \(0\le a\le b\le c\le N\) that are losing, for \(N=100\,000\). The statement gives one checkpoint: for \(N=29\) exactly \(1160\) of the \(\binom{32}{3}=4960\) triples in range are losing. ^[1]

The whole solution in one sentence: give every heap size the Grundy value of the one-heap game "remove a square", which turns Nim Square into ordinary Nim, so a triple is losing exactly when its three values have bitwise exclusive or zero, and count those triples with one sweep over the middle heap carrying a frequency table of values on each side.

Mathematical Approach

1. The game ends, and only one position is stuck

Because \(1=1^2\), a heap with at least one stone always offers a move, and every move removes at least one stone, so the total strictly decreases and play ends within \(a+b+c\) moves. The only position with no move is \((0,0,0)\); whoever faces it has lost, so \((0,0,0)\) is itself losing and belongs in the count. That is why the ranges start at \(0\).

2. One heap first: the subtract-a-square game

Strip the game to a single heap. Call a size cold when the player to move from it loses, hot otherwise; a size is hot exactly when some move reaches a cold size. Working upwards: \(0\) is cold, having no move; \(1\) is hot, reaching \(0\); \(2\) is cold, its only move reaching the hot \(1\); \(3\) and \(4\) are hot, reaching \(2\) and \(0\); \(5\) is cold, its moves reaching \(4\) and \(1\), both hot. Continuing, the cold sizes below \(100\) are

$$\begin{aligned} &0,\ 2,\ 5,\ 7,\ 10,\ 12,\ 15,\ 17,\ 20,\ 22,\\ &34,\ 39,\ 44,\ 52,\ 57,\ 62,\ 65,\ 67,\ 72,\\ &85,\ 95. \end{aligned}$$

This one-heap game is subtract-a-square, invented by Epstein and analysed by Golomb in 1966; its cold sizes are OEIS A030193, and up to \(100\,000\) there are \(2781\) of them. ^{[2] [6]} But cold sizes alone do not settle three heaps: that each heap separately is hot says nothing about the combination. We need a finer measurement.

3. Grundy values and the mex rule

That measurement is the Grundy value. The minimum excludant \(\operatorname{mex}\) of a set of non-negative integers is the smallest one it omits, so \(\operatorname{mex}\{0,1,3\}=2\) and \(\operatorname{mex}\varnothing=0\). ^[3] Put

$$g(n)=\operatorname{mex}\{\,g(n-k^2)\;:\;k\ge1,\ k^2\le n\,\},$$

so \(g(0)=0\), the empty case. Each value depends only on smaller ones, so the table fills from the bottom up. For \(n=0,1,\dots,29\) it reads

$$\begin{aligned} &0,1,0,1,2,0,1,0,1,2,\\ &0,1,0,1,2,0,1,0,1,2,\\ &0,1,0,1,2,3,2,3,4,5. \end{aligned}$$

Here \(g(1)=\operatorname{mex}\{0\}=1\), \(g(2)=\operatorname{mex}\{1\}=0\) and \(g(4)=\operatorname{mex}\{1,0\}=2\). The block \(0,1,0,1,2\) repeats up to \(n=24\) and then breaks, when the square \(25\) first becomes available: from \(25\) the options \(24,21,16,9,0\) have values \(2,1,1,2,0\), so \(g(25)=\operatorname{mex}\{0,1,2\}=3\), where without the new move the values \(\{1,2\}\) would give mex \(0\). The sequence is OEIS A014586. ^[5]

Two properties of the mex are all we use. (A) Every \(u\lt g(n)\) is the value of some option: were such a \(u\) omitted, the mex would be at most \(u\), not \(g(n)\). (B) No option has value \(g(n)\), omitted by definition. So a heap of value \(v\) can be moved to any smaller value and never to its own, exactly like a Nim heap of \(v\) stones. A move can also raise the value — one stone off a heap of \(2\) takes it from \(0\) to \(1\) — so the equivalence to Nim is one of outcomes, not of move sets. Property (A) with \(u=0\) also recovers Section 2, since \(g(n)\gt0\) means some option is cold and \(g(n)=0\) means none is, so by induction the cold sizes are those with \(g(n)=0\).

One bound is needed twice. A heap of size \(n\) has \(\lfloor\sqrt n\rfloor\) moves, and a set of at most \(m\) numbers cannot contain all \(m+1\) integers \(0,1,\dots,m\), so its mex is at most \(m\). Hence \(g(n)\le\lfloor\sqrt n\rfloor\le316\) for \(n\le100\,000\). The true maximum is far smaller, \(74\), first reached at \(n=99\,078\).

4. Three heaps become ordinary Nim

Write \(x\oplus y\) for the Nim-sum: write both numbers in binary and put a \(1\) in each position where the digits differ, a \(0\) where they agree. For instance \(5=101_2\) and \(3=011_2\) give \(5\oplus3=110_2=6\). Order and grouping do not matter, and because each digit obeys the same rule, \(x\oplus x=0\) and \(x\oplus0=x\). ^[4]

Claim. The position \((a,b,c)\) is losing for the next player exactly when

$$g(a)\oplus g(b)\oplus g(c)=0 .$$

Proof. Call that Nim-sum \(S\). If \(S=0\), a move changes one heap, say \(a\) to \(a'\); by (B) \(g(a')\ne g(a)\), so the new Nim-sum \(S\oplus g(a)\oplus g(a')=g(a)\oplus g(a')\) is not \(0\). Hence \(S=0\) leads only to \(S\ne0\). If \(S\ne0\), take the highest bit position in which \(S\) has a \(1\); at least one of the three values has a \(1\) there, say \(g(a)\). Then \(u=g(a)\oplus S\) has a \(0\) in that position while all higher bits are unchanged, so \(u\lt g(a)\); by (A) some move takes \(a\) to an \(a'\) with \(g(a')=u\), and the new Nim-sum is \(S\oplus g(a)\oplus u=0\). Hence \(S\ne0\) leads to \(S=0\).

The terminal \((0,0,0)\) has \(S=0\) and is losing, and every move strictly decreases the total number of stones, so induction on that total finishes the argument: from \(S=0\) the mover hands over only positions with \(S\ne0\), from each of which the opponent returns a position with \(S=0\), and the descent must end at \((0,0,0)\) with the first player to move.

This is the Sprague–Grundy theorem for our game: an impartial normal-play game with a finite game tree is equivalent to a Nim heap of its Grundy value, the value of a sum is the Nim-sum of the parts, and the mover loses exactly when it is \(0\); ^[3] the argument above is Bouton's classical proof for Nim. ^[4]

5. A completely worked example: heaps of at most four stones

Take \(N=4\), where \(g(0),\dots,g(4)=0,1,0,1,2\). Of the \(\binom73=35\) triples with \(0\le a\le b\le c\le4\), testing the Nim-sum leaves exactly \(12\) losing ones:

$$\begin{aligned} &(0,0,0),\ (0,0,2),\ (0,1,1),\ (0,1,3),\\ &(0,2,2),\ (0,3,3),\ (0,4,4),\ (1,1,2),\\ &(1,2,3),\ (2,2,2),\ (2,3,3),\ (2,4,4). \end{aligned}$$

Note that \((1,2,3)\) is losing although no two heaps are equal. The triple \((0,1,2)\) is absent: its values \(0,1,0\) have Nim-sum \(1\), and the proof's recipe moves the heap carrying the top bit — the heap of one stone — to value \(1\oplus1=0\), giving \((0,0,2)\), which is in the list.

6. Why the triples cannot be visited one at a time

The map \((a,b,c)\mapsto(a,b+1,c+2)\) sends our triples bijectively onto the strictly increasing triples from \(\{0,1,\dots,N+2\}\), so there are \(\binom{N+3}{3}\) of them: for \(N=100\,000\) that is \(166\,676\,666\,850\,001\), far too many to enumerate. What saves us is that the losing condition reads only the three Grundy values, and by Section 3 no value in the range exceeds \(74\). So we count in bulk, grouped by value.

7. The middle-heap sweep

Write \(L(N)\) for the number we want. A sorted triple \(a\le b\le c\) singles out its middle entry \(b\), and the losing condition rearranges to \(g(c)=g(a)\oplus g(b)\). Fix \(b\) and define two frequency tables,

$$P_b(x)=\#\{a: 0\le a\le b,\ g(a)=x\},$$

$$S_b(y)=\#\{c: b\le c\le N,\ g(c)=y\}.$$

Every losing triple with middle entry \(b\) arises by choosing \(a\le b\) of some value \(x\) and then \(c\ge b\) of the forced value \(x\oplus g(b)\), and distinct choices give distinct triples, so

$$L(N)=\sum_{b=0}^{N}\ \sum_{x}P_b(x)\,S_b\bigl(x\oplus g(b)\bigr).$$

Each losing triple is counted once, at the single \(b\) that is its middle entry. Both ranges are closed, \(a\le b\) and \(c\ge b\), so triples with repeated entries need no separate treatment and no symmetry factor appears. Moving from \(b\) to \(b+1\) adds one entry on the left and removes one on the right, and the order matters: when \(b\) is processed, \(g(b)\) must already be in \(P_b\) and still be in \(S_b\).

8. The same example through the sweep

For \(N=4\) the values \(0,1,2\) occur \(2,2,1\) times, so the right table starts at \((2,2,1)\). Both tables are written as counts for the values \(0,1,2\):

• \(b=0\), \(g(b)=0\): \(P=(1,0,0)\), \(S=(2,2,1)\), giving \(1\cdot2=2\): \((0,0,0)\), \((0,0,2)\).
• \(b=1\), \(g(b)=1\): \(P=(1,1,0)\), \(S=(1,2,1)\), giving \(1\cdot2+1\cdot1=3\): \((0,1,1)\), \((0,1,3)\), \((1,1,2)\).
• \(b=2\), \(g(b)=0\): \(P=(2,1,0)\), \(S=(1,1,1)\), giving \(2\cdot1+1\cdot1=3\): \((0,2,2)\), \((2,2,2)\), \((1,2,3)\).
• \(b=3\), \(g(b)=1\): \(P=(2,2,0)\), \(S=(0,1,1)\), giving \(2\cdot1=2\): \((0,3,3)\), \((2,3,3)\).
• \(b=4\), \(g(b)=2\): \(P=(2,2,1)\), \(S=(0,0,1)\), giving \(2\cdot1=2\): \((0,4,4)\), \((2,4,4)\).

The contributions add to \(12\), as in Section 5, and every triple there appears once.

9. One index that leaves the table

In the step \(b=1\) above, the value \(x=2\) calls for \(S_b(2\oplus1)=S_b(3)\), yet the tables reach only the largest value that occurs. Such a term is zero, but its index lies outside the table and must be skipped rather than read — and not rarely: for \(N=100\,000\) the largest value is \(74\), and \(75\) is not a power of two, so \(x\oplus g(b)\) can reach \(127\), and \(883\,354\) of the \(7\,500\,075\) inner steps, \(11.8\%\), have an index above \(74\). Only \(S_b\) is indexed this way, so padding it alone to \(128\) entries — zero, since those values never occur — removes the test and leaves the inner loop at \(75\) values.

10. How wide is the alphabet of values

Put \(G=1+\max_{0\le n\le N}g(n)\). By property (A) a heap of maximal value has options of every smaller value, so the values occurring are exactly \(0,1,\dots,G-1\) — no gaps — and the sweep makes \(G\) inner steps per middle heap. Section 3 bounds \(G\) by \(\lfloor\sqrt N\rfloor+1\), here \(317\), but the observed maxima are far smaller: \(5\) for \(n\le40\), \(14\) for \(n\le1000\), \(33\) for \(n\le10\,000\), \(74\) for \(n\le100\,000\). Eppstein observes the same and computes the table by another algorithm in \(\tilde O(mn)\) time, the tilde hiding logarithmic factors and \(m\) being the largest value rather than the size of the move set. ^[7]

How the Code Works

In C++ ^[9] main parses the options, runs run_checkpoints unless --skip-checkpoints was given, and prints solve(limit); each solve builds the value table, then sweeps. build_grundy(limit) is Section 3 from the bottom up: an array g of limit + 1 entries starting at \(0\), and a scratchpad seen of \(512\) entries starting at \(-1\). For the current \(n\) the loop for (int k = 1; k * k <= n; ++k) stamps every reachable value with seen[val] = n, and while (seen[g[n]] == n) ++g[n]; then walks up from \(0\) to the first value not stamped with this \(n\): the mex. Stamping rather than clearing flags saves a pass per heap size. The code enlarges seen if a value reaches its length; by Section 3 no value exceeds \(316\), so that branch never runs here.

solve(limit) is Section 7. It scans the table once for max_g, allocates prefix and suffix with max_g + 1 counters, and fills suffix with the value counts of the whole range. The sweep then runs b upwards in the order Section 7 demands: ++prefix[gb] first, so \(b\) counts as a possible \(a\); then the inner loop over ga, forming gc = ga ^ gb, applying the guard if (gc > max_g) continue; of Section 9 and adding prefix[ga] * suffix[gc]; and only afterwards --suffix[gb], so \(b\) counts as a possible \(c\) during its own step and not later.

The C++ file carries its own checks. brute_small(limit) is the definition itself, three nested loops over \(a\le b\le c\) testing the Nim-sum. run_checkpoints demands solve(29) == 1160, the statement's value, and solve(40) == brute_small(40), over all \(\binom{43}{3}=12\,341\) triples of that range. parse_arguments accepts --limit= and --skip-checkpoints, which is how the small cases above were produced.

The three languages differ only in machinery. C++ accumulates in a long long and keeps the tables in std::vector, so growing seen is one resize; Java mirrors it with int[] and long[] and grows seen by hand, with System.arraycopy then Arrays.fill; Python, whose integers are unbounded, needs no overflow analysis and grows seen with extend, appending val + 64 entries instead of resizing to that length — harmless, as the branch is unreachable here. Neither port runs checkpoints; Java's parameterless solve() sets the bound to \(100\,000\), while Python's solve(limit=100000) only defaults to it and accepts another.

Complexity and Verification

Building the value table costs one stamping step per move, so \(\lfloor\sqrt n\rfloor\) steps for the size \(n\). Comparing the sum with \(\int_0^N\sqrt x\,\mathrm dx\) and absorbing the dropped fractional parts into the error term,

$$\sum_{n=1}^{N}\lfloor\sqrt n\rfloor=\tfrac23N^{3/2}+O(N)$$

which is exactly \(21\,032\,170\) steps for \(N=100\,000\) against \(\tfrac23N^{3/2}\approx21\,081\,851\). The mex loops add \(2\,469\,832\) increments, one per unit of value, plus the \(N\) failed tests that end them: \(2\,569\,832\) tests in all. The sweep costs \(G\) inner steps per heap size, \((N+1)\cdot75=7\,500\,075\) together. The method is thus \(O(N\sqrt N)\) in time, dominated by the table, with an \(O(NG)\) counting stage, and \(O(N+G)\) stored entries: \(100\,001\) values, about \(0.4\) MB where int is 32 bits, beside \(512\) scratchpad entries and two tables of \(75\) counters. Python holds references to integer objects, so the same counts cost it more.

Fixed-width accumulators are safe, and Section 6 proves it: \(L(N)\) cannot exceed \(\binom{N+3}{3}=166\,676\,666\,850\,001\), far below \(2^{63}-1\), and no product \(P_b(x)\,S_b(y)\) exceeds \((N+1)^2=10\,000\,200\,001\). On the machine used here the C++ program takes about \(0.015\) s including its checkpoints, Java about \(0.05\) s and Python about \(1.35\) s.

Checked by computation for this article: a brute force over all sorted triples agrees with the sweep for every \(N\le60\), among them \(12\) at \(N=4\) as in Section 5, \(100\) at \(N=10\), \(1160\) at \(N=29\) — the statement's checkpoint — \(2632\) at \(N=40\) and \(7160\) at \(N=60\). The first thirty values of \(g\) match the opening terms of A014586, ^[5] the \(21\) cold sizes below \(100\) match those of A030193, ^[6] and \(g(n)\le\lfloor\sqrt n\rfloor\) holds for every \(n\le100\,000\), where the maximum is \(74\).

The count was then recomputed by a different route, to test the sweep rather than repeat it. Let \(X\) hold the ordered triples in \([0,N]^3\) of Nim-sum \(0\). Sorting collects them into the classes we want, but dividing \(|X|\) by six would be wrong, since \((0,1,1)\) has three orderings and \((0,0,0)\) one; Burnside's lemma instead averages, over the six permutations of the coordinates, the triples each leaves fixed. ^[8] Let \(f(x)\) count the sizes in \([0,N]\) of value \(x\), zero for a value that never occurs, and put \(z=f(0)\). Then

$$|X|=\sum_{x,y}f(x)\,f(y)\,f(x\oplus y),$$

a sum of only \(75^2=5625\) terms, so no fast transform is worth the trouble. A transposition fixes the triples \((a,a,c)\) with \(g(c)=0\), so \((N+1)z\) of them, and a three-cycle fixes \((a,a,a)\) with \(g(a)=0\), so \(z\); the three transpositions and two three-cycles then give

$$L(N)=\frac{|X|+3(N+1)z+2z}{6}.$$

Each bound has its own table, so \(z=2\) at \(N=4\) and \(z=2781\) at \(N=100\,000\). With those, the formula agrees with the sweep for \(N=4,10,29,40,100,1000,10\,000\) and \(100\,000\). Finally the C++, Python and Java programs were each run to completion and produce byte-identical output.

Footnotes and References

The references below supply the background facts used above; the middle-heap sweep, the Nim-sum proof, the worked examples, the operation counts and the Burnside recount are derived or computed here.

1. Project Euler 310 — Nim Square. The official statement describes three-heap normal-play Nim in which only a square number of stones may be removed, writes a position as the triple (a, b, c), gives 1160 losing positions for the next player when 0 ≤ a ≤ b ≤ c ≤ 29, and asks for that count with the bound 100000.
2. Wikipedia — Subtract a square. The one-heap game removes a non-zero square number of coins; the article credits its invention to Richard A. Epstein and its analysis to Solomon W. Golomb (1966), lists its cold positions 0, 2, 5, 7, 10, 12, 15, … as A030193, and analyses the several-heap version with the Sprague–Grundy theorem, whose values begin 0, 1, 0, 1, 2, 0, 1, 0, 1, 2, ….
3. Wikipedia — Sprague–Grundy theorem. Every impartial normal-play game with a finite game tree is equivalent to a one-heap Nim game of its Grundy value (an infinite game needs a generalisation of Nim), that value is zero exactly for a loss for the next player, the value of a sum is the Nim-sum of the summands, and the proof defines the value as the mex of the option values.
4. Wikipedia — Nim. The Nim-sum is the binary digital sum of the heap sizes without carries, their bitwise exclusive or; Bouton's theorem of 1901 says the player about to move wins in normal play exactly when that Nim-sum is not zero.
5. OEIS A014586 — Nim-Grundy function for the subtract-a-square game. The terms are the mex of the values at n − k² for k² ≤ n; the first thirty are those tabulated in Section 3, and the value vanishes exactly at the terms of A030193.
6. OEIS A030193 — losing positions of the subtract-a-square game. It lists the cold positions of the one-heap game, those from which the player to move cannot force a win; its opening terms are the twenty-one values below 100 given in Section 2.
7. David Eppstein, Faster Evaluation of Subtraction Games (arXiv:1804.06515). The paper computes the Nim-values of all heap sizes up to n in time Õ(mn) with m the largest Nim-value occurring, and records that for subtract-a-square those values are experimentally far smaller than the subtraction set.
8. Wikipedia — Burnside's lemma. For a finite group G acting on a finite set X the number of orbits is the average number of fixed points, |X/G| = (1/|G|) Σ |X^g|, applied above to the symmetric group on three letters acting on the ordered Nim-sum-zero triples.
9. C++ — Euler310.cpp. The linked immutable C++ revision contains build_grundy, solve, brute_small and run_checkpoints, with the --limit= and --skip-checkpoints options. It is the implementation record for this article.