Problem 391: Hopping Game

View on Project Euler

Project Euler Problem 391 Solution

EulerSolve provides an optimized solution for Project Euler Problem 391, Hopping Game, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For each \(n\), the repository code computes a value \(m(n)\) and then adds \(m(n)^3\) to the global sum $$\sum_{n=1}^{1000} m(n)^3.$$ The important point is that the local solutions do not use a short closed form. Instead they model the game by repeatedly composing residue maps on \(R_n=\{0,1,\dots,n\}\) and stop as soon as the diagonal row becomes constant. The checkpoints hard-coded in all three implementations are \(m(2)=2\), \(m(7)=1\), \(m(20)=4\), and \(\sum_{n=1}^{20} m(n)^3=8150\). Mathematical Approach Step 1: Residues and the Overflow-to-Zero Shift Fix \(n\) and set \(M=n+1\). Every table entry is a function value on the residue set \(R_n=\{0,1,\dots,n\}\). The code does not use ordinary addition modulo \(M\). Instead it uses the truncated shift $$\sigma_k(x)=\begin{cases} x+k, & x+k \lt M,\\ 0, & x+k \ge M, \end{cases}\qquad 1 \le k \le M.$$ That distinction matters: ordinary modular addition would send \(M-1+k\) to \(k-1\), but the program sends every overflow directly to \(0\). This is exactly what the line if (value >= step) value = 0; implements in C++, Python, and Java. Step 2: The Clean Recurrence Behind the Two Tables Let \(F_s^{(m)}:R_n \to R_n\) denote the map stored in row \(m\) after stage \(s\) has finished....

Detailed mathematical approach

Problem Summary

Write the integers \(0,1,2,\dots,k\) in binary and count every \(1\) digit; call the total \(s_k\). For \(k=5\) the strings are \(0,1,10,11,100,101\), which contain seven \(1\) digits, so \(s_5=7\), and the set \(S=\{s_k:k\ge0\}\) begins \(0,1,2,4,5,7,9,12,\dots\) [2] Two players then share a counter \(c\) that starts at \(0\). A bound \(n\) is fixed before play; a move adds some integer between \(1\) and \(n\) to \(c\), and the resulting value of \(c\) must again belong to \(S\). A player with no legal move loses. Let \(M(n)\) be the highest opening move with which the first player can force a win, and \(M(n)=0\) when no opening move can. The statement gives \(M(2)=2\), \(M(7)=1\), \(M(20)=4\) and \(\sum M(n)^3=8150\) for \(1\le n\le20\), and asks for the same sum over \(1\le n\le1000\). [1]

The whole solution in one sentence: the losing positions can be listed by walking backwards along \(S\) from beyond the last position the game can reach, accumulating the gaps and resetting the running total to \(0\) whenever it would pass \(n\), and because the gaps of \(S\) are the binary digit sums, that astronomically long walk folds into a short recursion over blocks of length \(2^m\) whose outcome stops depending on where the walk began.

Mathematical Approach

The game looks like a game about a set of numbers, but all the information it uses sits in the gaps of that set. We turn the game into a walk along those gaps, show that the walk needs one number of memory, and then use the self-similarity of the binary digit sum to replace that walk, which is \(2^{1001}\) steps long at \(n=1000\), by a few hundred tabulated maps.

1. The gaps of \(S\) are binary digit sums

Let \(w(k)\) be the number of \(1\) digits of \(k\) in binary, the binary weight or population count of \(k\). [4] By definition \(s_k\) accumulates these weights,

$$s_k=\sum_{j=0}^{k}w(j),$$

so consecutive members of \(S\) differ by

$$s_k-s_{k-1}=w(k)\qquad(k\ge1).$$

Since \(w(k)\ge1\) for \(k\ge1\), the sequence \(s_0,s_1,s_2,\dots\) is strictly increasing: every element of \(S\) has exactly one index. The weights \(w(1),w(2),w(3),\dots\) are \(1,1,2,1,2,2,3,\dots\), and accumulating them reproduces \(0,1,2,4,5,7,9,12,\dots\) [2] [3] We shall need only the halving rule \(w(2r)=w(r)\), \(w(2r+1)=w(r)+1\): appending a binary digit adds that digit to the weight. [3]

2. The game as a walk on the indices of \(S\)

Because each value of \(c\) is some \(s_k\), we may describe a position by its index \(k\) instead of by the counter. A move from index \(k\) goes to an index \(j\gt k\) and costs \(s_j-s_k\), the sum of the weights \(w(k+1),\dots,w(j)\); it is legal exactly when that cost is at most \(n\). The starting position is \(k=0\). Index \(k\) has no legal move precisely when even the cheapest one, the step to \(k+1\), is too expensive:

$$w(k+1)\gt n.$$

3. Play stays in a finite part of the board

The number \(2^{n+1}-1\) is a string of \(n+1\) ones, so \(w(2^{n+1}-1)=n+1\gt n\). Hence the index

$$K_n=2^{n+1}-2$$

has no legal move, and no move can step over it either: leaving an index \(k\le K_n\) for an index \(j\gt K_n\) would have to pay that weight along the way. A game started at \(0\) therefore lives inside the indices \(0,1,\dots,K_n\), and it ends, because every move increases the index.

That bound is far too large for direct computation. Among the \(2^{n+1}\) integers below \(2^{n+1}\), each of the \(n+1\) bit positions carries a \(1\) exactly \(2^n\) times, so \(s_{2^{n+1}-1}=(n+1)2^n\); subtracting the weight \(n+1\) of the last of them leaves \(s_{K_n}=(n+1)(2^n-1)\). For \(n=20\) the board already holds \(2\,097\,151\) indices counting \(0\), and for \(n=1000\) about \(2^{1001}\). Everything below avoids that enumeration.

4. Losing and winning positions

The game is impartial and follows the normal play convention: both players have the same moves, and a player who cannot move loses. [5] Call an index a loss when the player who must move from it loses under perfect play, and a win otherwise. Unwinding the definition gives the usual backward recursion: an index is a loss exactly when every legal move from it lands on a win, and an index with no legal move is a loss. Write \(\mathcal L\) for the set of losing indices. The first player wins by moving to a loss, so

$$M(n)=\max\{s_j\le n:\ j\in\mathcal L\},$$

with \(M(n)=0\) if no such \(j\ge1\) exists. The opening move is just the counter value \(s_j\) it produces.

5. Only the next losing position matters

Let \(0\le k\lt K_n\) and let \(q\) be the smallest losing index with \(q\gt k\); it exists because \(K_n\) is a loss lying to the right of \(k\). We claim

$$k\in\mathcal L\iff s_q-s_k\gt n.$$

If \(s_q-s_k\le n\), then moving from \(k\) to \(q\) is legal and reaches a loss, so \(k\) is a win. If \(s_q-s_k\gt n\), then every legal target \(j\) satisfies \(s_j\le s_k+n\lt s_q\), hence \(j\lt q\); by the minimality of \(q\) all indices strictly between \(k\) and \(q\) are wins, so every legal move from \(k\) lands on a win and \(k\) is a loss. The terminal index \(K_n\) is the base case of this recursion: it is a loss because it has no legal move.

So the losing indices are determined from the right: after a loss, the next loss going left is the last index whose counter value is more than \(n\) below it, and two losing counter values therefore differ by at least \(n+1\).

6. One running total replaces the whole table

For \(0\le k\le K_n\), define the distance to the nearest loss at or after \(k\):

$$u_k=s_{q(k)}-s_k,\qquad q(k)=\min\{q\ge k:\ q\in\mathcal L\}.$$

Section 5 says \(u_k=0\) when \(k\) is a loss and \(0\lt u_k\le n\) when \(k\) is a win, so \(u_k\) always lies in \(\{0,1,\dots,n\}\). Introduce the capped shift on that set,

$$\sigma_j(u)=\begin{cases}u+j,&u+j\le n,\\0,&u+j\gt n,\end{cases}$$

which adds \(j\) but collapses to \(0\) on overflow instead of wrapping round. Then

$$u_{k-1}=\sigma_{w(k)}(u_k).$$

Indeed \(u_k+w(k)=s_{q(k)}-s_{k-1}\). If that is at most \(n\), the nearest loss has not changed, \(q(k-1)=q(k)\), and \(u_{k-1}\) is exactly this sum; if it exceeds \(n\), then by Section 5 the index \(k-1\) is itself a loss and \(u_{k-1}=0\). Finally \(q(0)\) is a loss whose counter value is \(u_0=s_{q(0)}\le n\), and by the spacing proved in Section 5 it is the only losing value in \(\{0,1,\dots,n\}\); together with Section 4 this gives

$$M(n)=u_0.$$

7. Worked example: \(n=2\) by hand

Here \(K_2=2^3-2=6\) and \(S\) up to that point is \(0,1,2,4,5,7,9\), with \(s_6=9\). Start the walk at \(k=7\) with any value, say \(u_7=0\): since \(w(7)=3\gt2\), the first step resets anyway. The weights \(w(7),\dots,w(1)\) are \(3,2,2,1,2,1,1\), giving

  • \(k=7\): \(0+3\gt2\), so \(u_6=0\) — index \(6\) (value \(9\)) is a loss;
  • \(k=6\): \(0+2\le2\), so \(u_5=2\);
  • \(k=5\): \(2+2\gt2\), so \(u_4=0\) — index \(4\) (value \(5\)) is a loss;
  • \(k=4\): \(0+1\le2\), so \(u_3=1\);
  • \(k=3\): \(1+2\gt2\), so \(u_2=0\) — index \(2\) (value \(2\)) is a loss;
  • \(k=2\): \(0+1\le2\), so \(u_1=1\);
  • \(k=1\): \(1+1\le2\), so \(u_0=2\).

Hence \(M(2)=2\), and the losing counter values are \(2\), \(5\) and \(9\). Each can be checked directly: from \(9\) neither \(10\) nor \(11\) lies in \(S\); from \(5\) the only move is to \(7\), from which one can move to the loss \(9\); from \(2\) the only move is to \(4\), from which one can move to the loss \(5\). Opening with \(2\) therefore wins, while \(1\) leaves the opponent a win.

The sample line of play in the statement, \(c=0\to4\to9\to12\) with \(n=5\), is an illustration rather than good play: the value \(4\) is index \(3\), which is a win for the player who moves next, while \(M(5)=2\).

8. At most one winning opening

Since two losing values differ by at least \(n+1\), exactly one of them lies in \(\{0,1,\dots,n\}\), namely \(u_0\). The first player therefore has at most one winning opening move, and the word "highest" in the definition of \(M(n)\) never has to choose between candidates. If that unique value is \(0\), index \(0\) is itself a loss and no opening move wins, so both readings give \(M(n)=0\); this happens for \(n=1\). For \(1\le n\le20\) the values \(M(n)\) are \(0,2,1,4,2,2,1,7,5,5,\)\({}7,4,13,7,15,2,9,1,7,4\), whose cubes sum to \(8150\), matching the statement.

9. The walk forgets where it started

Composing the capped shifts turns the walk into a function; composition feeds one output into the next, \((f\circ g)(u)=f(g(u))\), so \(g\) acts first. For \(L\ge1\) put

$$F_L=\sigma_{w(1)}\circ\sigma_{w(2)}\circ\cdots\circ\sigma_{w(L)},$$

the map that sends a value at index \(L\) to the value it becomes at index \(0\). Composition is associative, so any grouping of the factors computes the same map. [6] The map \(\sigma_{n+1}\) is constant \(0\), because \(u+n+1\gt n\) always. Taking \(L=2^{n+1}-1\) makes \(\sigma_{w(L)}=\sigma_{n+1}\) the first factor applied, so

$$F_L\equiv M(n)\qquad(L\ge2^{n+1}-1),$$

a constant map: the walk over that many indices wipes out its input and returns \(u_0\).

10. Halving the blocks

Now the self-similarity of \(w\) enters. For \(r=0,1,2,3\) the weights \(w(r)\) are \(0,1,1,2\); raised by \(1\) they read \(1,2,2,3\), whose upper half \(2,3\) is the lower half \(1,2\) with \(1\) added. For \(m\ge0\) and \(j\ge1\), let \(B_{m,j}\) be the map obtained by walking \(r=2^m-1,2^m-2,\dots,0\) and applying \(\sigma_{w(r)+j}\) at each step — a block of \(2^m\) consecutive indices whose weights are all raised by \(j\). Splitting the block at its midpoint and using \(w(2^{m-1}+r)=w(r)+1\) for \(0\le r\lt2^{m-1}\), which is the halving rule of Section 1, the upper half is a copy of the lower half with \(j\) replaced by \(j+1\). The upper half is walked first, so

$$B_{0,j}=\sigma_j,\qquad B_{m,j}=B_{m-1,j}\circ B_{m-1,j+1}.$$

Since \(w(p2^m+r)=w(p)+w(r)\) for \(0\le r\lt2^m\), the block of indices with binary prefix \(p\) is the one with \(j=w(p)\): the parameter \(j\) carries the weight of the prefix, while \(m\) counts the suffix bits. This is the kind of base-\(2\) recursion that makes sequences such as \(w\) tractable. [7]

11. One recurrence for prefixes

What we need is not a single block but the prefix \(F_{2^m-1}\), which omits the index \(0\). So let \(C_{m,j}\) be the map obtained by walking \(k=2^m-1,\dots,1\) and applying \(\sigma_{w(k)+j-1}\), with \(C_{0,j}\) the identity. Cutting that walk into the blocks \([2^i,2^{i+1})\) for \(i=m-1,\dots,0\) and applying the halving rule to each gives

$$C_{m,j}=B_{0,j}\circ B_{1,j}\circ\cdots\circ B_{m-1,j},$$

and comparing this with \(B_{m,j}\), which is the same walk with the index \(0\) appended at the end, yields

$$B_{m,j}=\sigma_j\circ C_{m,j+1}.$$

Peeling the last factor off that product gives \(C_{m,j}=C_{m-1,j}\circ B_{m-1,j}\) for \(m\ge1\), and substituting the identity above removes the block maps entirely:

$$\begin{aligned}C_{0,j}&=\mathrm{id},\\C_{m,j}&=C_{m-1,j}\circ\sigma_j\circ C_{m-1,j+1}.\end{aligned}$$

By construction \(C_{m,1}=F_{2^m-1}\), the map we are hunting.

12. The order of evaluation and when to stop

The recurrence needs \(C_{m-1,j}\) and \(C_{m-1,j+1}\), whose parameter sums are \(m+j-1\) and \(m+j\). Grouping the maps by \(m+j\) and taking the groups in increasing order therefore makes each group depend only on itself and its predecessor. With the group index \(s=m+j-1\), group \(s\) is \(C_{1,s},C_{2,s-1},\dots,C_{s,1}\), whose last member is \(F_{2^s-1}\), so two arrays of maps suffice.

Splitting the walk at the midpoint also gives \(F_{2^s-1}=F_{2^{s-1}-1}\circ B_{s-1,1}\), so the image of \(F_{2^s-1}\) is contained in the image of \(F_{2^{s-1}-1}\): the set of reachable values can only shrink as \(s\) grows. Once \(F_{2^s-1}\) is constant it stays constant with the same value, which by Section 9 must be \(M(n)\). We may therefore stop at the first \(s\) for which \(F_{2^s-1}\) is constant, and Section 9 guarantees that this happens by \(s=n+1\) at the latest.

In practice it happens far earlier, because each reset merges trajectories that started at different values. For \(n=1,2,20,100,500,1000\) the first constant group is \(s=2,3,8,14,18,21\), and over \(1\le n\le1000\) the largest is \(25\), at \(n=905\). For \(n=20\) the \(255\) indices \(k=1,\dots,2^8-1\) therefore already determine \(M(20)=4\), instead of the \(2\,097\,151\) that Section 9 asks for.

13. The recursion for \(n=2\)

With \(n=2\) the maps are triples \((f(0),f(1),f(2))\) and the capped shifts are \(\sigma_1=(1,2,0)\), \(\sigma_2=(2,0,0)\), \(\sigma_3=(0,0,0)\). Group \(s=1\) holds only \(C_{1,1}=\sigma_1=(1,2,0)\), which is not constant. Group \(s=2\) holds \(C_{1,2}=\sigma_2=(2,0,0)\) and then \(C_{2,1}=C_{1,1}\circ\sigma_1\circ C_{1,2}=(1,2,2)\), which is \(F_3\), the walk over the indices \(3,2,1\) of Section 7; it still depends on its input, so this prefix alone does not yet determine \(M(2)\). Group \(s=3\) holds \(C_{1,3}=\sigma_3=(0,0,0)\), \(C_{2,2}=(0,0,0)\) and \(C_{3,1}=(2,2,2)\). The last map is constant, so \(M(2)=2\), in agreement with the hand walk.

How the Code Works

All three programs implement Sections 6, 11 and 12 directly, and none of them ever computes a binary weight: the self-similarity is baked into the recurrence, so the weights survive only as shift amounts.

The C++ class Solver [9] holds two flat buffers table_a and table_b of \(1024\cdot1024\) integers with pointers a and b into them. A map is one row of \(n+1\) entries, placed by index(m, v) = (m << SHIFT) + v at a fixed stride of \(1024\). Row \(m\) of b holds \(C_{m,s-m+1}\) for the group being built, row \(m\) of a holds \(C_{m,s-m}\) from the previous group, and row \(0\) of both is the identity \(C_{0,j}\), written once in the constructor.

compute_m(n) sets step to \(n+1\), the overflow threshold of the capped shift, and loops over the groups s. For each row m it sets add = s - m + 1, the second parameter \(j\) of \(C_{m,j}\), and evaluates the recurrence of Section 11 pointwise: it reads the value of \(C_{m-1,j+1}\) at \(v\) from row \(m-1\) of b, adds add, replaces the result by \(0\) if it reaches step — the one line that is \(\sigma_j\) — and looks the outcome up in row \(m-1\) of a, which is \(C_{m-1,j}\).

Once the group is filled, the code scans row s for an entry differing from its first. Finding none means \(C_{s,1}=F_{2^s-1}\) is constant, and that constant is returned as \(M(n)\): the stopping test of Section 12. Otherwise swap(a, b) turns the finished group into the previous one. The return a[index(step, 0)] after the loop is a safety net that Section 9 makes unreachable, since group \(n+1\) is always constant.

Reusing the buffers across values of \(n\) is safe: group \(s\) reads only rows \(0,\dots,s-1\), written by the current or previous group of the same call, and row \(0\) is never overwritten. The stride also covers the worst case, since the highest row that could be needed, \(n+1=1001\), ends at offset \(1\,026\,024\) of \(1\,048\,576\); raising MAX_N past \(1022\) would demand a larger SHIFT.

The C++ main and the Java solve accumulate m * m * m over \(n=1,\dots,1000\) and, at \(n=20\), refuse to continue unless \(M(2)=2\), \(M(7)=1\), \(M(20)=4\) and the partial sum is \(8150\); C++ exits with an error and Java throws. The Python solve carries no such check. It uses the same stride, swaps with a, b = b, a, and puts the unreachable fallback in the else branch of its for loop; Java's computeM keeps tableA and tableB as fields and writes the swapped references back before returning. The accumulator is int64_t in C++ and long in Java, while Python's integers are unbounded.

Complexity and Verification

Let \(t(n)\) be the group at which the stopping test fires. Group \(s\) fills \(s\) rows of \(n+1\) entries, so one value of \(n\) costs

$$\sum_{s=1}^{t(n)}s\,(n+1)=\tfrac12(n+1)\,t(n)\bigl(t(n)+1\bigr)$$

constant-time updates, plus \(O(n)\) per group for the constancy scan. One value of \(n\) therefore costs \(O(n\,t(n)^2)\), which the proven bound \(t(n)\le n+1\) turns into \(O(n^3)\) in the worst case. Over the tested range \(1\le n\le1000\) we measured \(t(n)\le25\), capping one value at \(325(n+1)\) updates: \(n=1000\) has \(t=21\) and \(231\,231\) updates, and all of them together need \(99\,769\,950\). The C++ and Java buffers hold \(2\cdot1024\cdot1024\) four-byte entries, \(8\) MiB of payload, while Python's two lists hold that many object references, \(16\) MiB of slots alone; the allocation is fixed for \(n\le1000\) even though only rows \(0,\dots,t(n)\) are used.

Every \(M(n)\) is at most \(n\), so the final sum is below \(\sum_{n\le1000}n^3=250\,500\,250\,000\), far inside the signed 64-bit range of C++ and Java. [8] The three programs print identical output, and C++ and Java finish in a fraction of a second while Python takes a few seconds.

Beyond the checkpoints built into the C++ and Java programs, the following were verified for this article. A direct solver that labels every index \(0,\dots,2^{n+1}-2\) by the backward recursion of Section 4 alone reproduces \(M(n)\) for all \(1\le n\le20\), including the four checkpoints, and confirms that exactly one losing counter value lies in \(\{0,\dots,n\}\) and equals \(M(n)\). The running total of Section 6, walked over the whole range \(k=2^{n+1}-1,\dots,1\), agrees with the group recursion for all \(1\le n\le28\); at \(n=28\) that walk visits \(536\,870\,911\) indices. For \(n\le10\) the rows produced by the code match \(C_{m,s-m+1}\) computed from the definition entry by entry; the identities of Sections 10 and 11 were checked as explicit tables of maps for \(1\le n\le8\), and the image of \(F_{2^s-1}\) was confirmed to shrink and, once a singleton, to stay equal to \(M(n)\). These checks are not part of the published solutions.

Footnotes and References

The references below supply the background facts used above; the reduction to a capped running total, the block identities and the group recursion are derived in this article.

  1. Project Euler 391 — Hopping Game. The official statement defines s_k as the number of 1 digits used in writing 0 to k in binary, notes s_5 = 7, describes the game in which each move adds between 1 and n to the counter and must keep it inside S, defines M(n) as the highest winning opening move or 0, and gives M(2) = 2, M(7) = 1, M(20) = 4 with the cube sum 8150 for n up to 20 before asking for the same sum up to 1000.
  2. OEIS A000788 — "Total number of 1's in binary expansions of 0, ..., n." Its terms 0, 1, 2, 4, 5, 7, 9, 12, … are the set S of the problem. The entry gives a(2^m − 1) = m·2^(m−1), used in Section 3, and the leading-bit split a(2^m + r) = m·2^(m−1) + r + 1 + a(r) for r < 2^m behind the self-similarity of Section 10.
  3. OEIS A000120 — "Ones-counting sequence: number of 1's in binary expansion of n (or the binary weight of n)." Its first term is w(0) = 0, and its terms from index 1, namely 1, 1, 2, 1, 2, 2, 3, …, are the gaps of S; the entry states the halving rule a(2n) = a(n), a(2n+1) = a(n) + 1 used in Sections 1 and 10.
  4. Wikipedia — Hamming weight. The article defines the weight of a binary number as "the number of bits set to 1, or the digit sum of the binary representation", lists the names population count and popcount, and notes that the weight of a number can differ from that of its half only in the last bit, the halving rule of Section 1.
  5. Wikipedia — Sprague–Grundy theorem. The article states the normal play convention, "a player who cannot move loses", describes an impartial game as one in which "each player is allowed exactly the same set of moves", and distinguishes positions won by the next player from those won by the previous player. Section 4 unwinds these definitions into the backward recursion.
  6. Wikipedia — Function composition. The article states that composition of functions "is always associative", so the long product of capped shifts may be regrouped into blocks as in Sections 10 to 12, and that chains of transformations of a set into itself "have the algebraic structure of a monoid".
  7. Wikipedia — k-regular sequence. A k-regular sequence is described as "a sequence satisfying linear recurrence equations that reflect the base-k representations of the integers", and the article lists the ruler, Thue–Morse and Stern–Brocot sequences as 2-regular examples. The halving rule of Section 10 is a recursion of this base-2 kind, applied to maps rather than numbers.
  8. Java Language Specification, §4.2.1. The section gives the range of long as −9223372036854775808 to 9223372036854775807, which contains the bound on the final sum computed above, and the range of int as −2147483648 to 2147483647, enough for every table entry.
  9. C++ — Euler391.cpp. The linked immutable revision contains the Solver class with MAX_N = 1000 and SHIFT = 10, the group loop of compute_m with its capped shift and constancy test, and the validation at n = 20. It is the implementation record for this article; the extra comparisons above are not packaged in it.

Problem 391 source code

C++

#include <cstdint>
#include <iostream>
#include <utility>
#include <vector>

using namespace std;

class Solver {
public:
  static constexpr int MAX_N = 1000;

  Solver() : table_a(STRIDE * STRIDE, 0), table_b(STRIDE * STRIDE, 0) {
    a = table_a.data();
    b = table_b.data();
    for (int v = 0; v <= MAX_N; ++v) {
      a[index(0, v)] = v;
      b[index(0, v)] = v;
    }
  }

  int compute_m(int n) {
    const int step = n + 1;
    for (int s = 1; s <= step; ++s) {
      for (int m = 1; m <= s; ++m) {
        const int add = s - m + 1;
        const int prev_row = (m - 1) << SHIFT;
        const int cur_row = m << SHIFT;
        for (int v = 0; v <= n; ++v) {
          int value = b[prev_row + v] + add;
          if (value >= step) {
            value = 0;
          }
          b[cur_row + v] = a[prev_row + value];
        }
      }

      const int row = s << SHIFT;
      const int value = b[row];
      int v = 1;
      while (v <= n && b[row + v] == value) {
        ++v;
      }
      if (v > n) {
        return value;
      }
      swap(a, b);
    }
    return a[index(step, 0)];
  }

private:
  static constexpr int SHIFT = 10;
  static constexpr int STRIDE = 1 << SHIFT;

  vector<int> table_a;
  vector<int> table_b;
  int *a;
  int *b;

  static constexpr int index(int m, int v) { return (m << SHIFT) + v; }
};

int main() {
  Solver solver;
  int m2 = -1;
  int m7 = -1;
  int m20 = -1;
  int64_t sum = 0;

  for (int n = 1; n <= Solver::MAX_N; ++n) {
    const int m = solver.compute_m(n);
    if (n == 2) {
      m2 = m;
    }
    if (n == 7) {
      m7 = m;
    }
    if (n == 20) {
      m20 = m;
    }
    sum += static_cast<int64_t>(m) * m * m;

    if (n == 20) {
      if (m2 != 2 || m7 != 1 || m20 != 4 || sum != 8150) {
        cerr << "Validation failed.\n";
        return 1;
      }
    }
  }

  cout << sum << '\n';
  return 0;
}

Python

def solve():
    MAX_N = 1000
    SHIFT = 10
    STRIDE = 1 << SHIFT

    a = [0] * (STRIDE * STRIDE)
    b = [0] * (STRIDE * STRIDE)

    for v in range(MAX_N + 1):
        a[(0 << SHIFT) + v] = v
        b[(0 << SHIFT) + v] = v

    total = 0
    for n in range(1, MAX_N + 1):
        step = n + 1
        for s in range(1, step + 1):
            for m in range(1, s + 1):
                add = s - m + 1
                prev_row = (m - 1) << SHIFT
                cur_row = m << SHIFT
                for v in range(n + 1):
                    value = b[prev_row + v] + add
                    if value >= step:
                        value = 0
                    b[cur_row + v] = a[prev_row + value]

            row = s << SHIFT
            val = b[row]
            v = 1
            while v <= n and b[row + v] == val:
                v += 1
            if v > n:
                m_val = val
                break
            a, b = b, a
        else:
            m_val = a[(step << SHIFT)]

        total += m_val * m_val * m_val

    return str(total)

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

Java

public class Euler391 {
    static class Solver {
        static final int SHIFT = 10;
        static final int STRIDE = 1 << SHIFT;

        int[] tableA = new int[STRIDE * STRIDE];
        int[] tableB = new int[STRIDE * STRIDE];

        Solver() {
            for (int v = 0; v <= 1000; v++) {
                int idx = v; // index(0, v)
                tableA[idx] = v;
                tableB[idx] = v;
            }
        }

        int computeM(int n) {
            int[] a = tableA;
            int[] b = tableB;
            int step = n + 1;

            for (int s = 1; s <= step; ++s) {
                for (int m = 1; m <= s; ++m) {
                    int add = s - m + 1;
                    int prevRow = (m - 1) << SHIFT;
                    int curRow = m << SHIFT;
                    for (int v = 0; v <= n; ++v) {
                        int value = b[prevRow + v] + add;
                        if (value >= step) {
                            value = 0;
                        }
                        b[curRow + v] = a[prevRow + value];
                    }
                }

                int row = s << SHIFT;
                int value = b[row];
                int v = 1;
                while (v <= n && b[row + v] == value) {
                    ++v;
                }
                if (v > n) {
                    tableA = a;
                    tableB = b;
                    return value;
                }

                int[] temp = a;
                a = b;
                b = temp;
            }

            tableA = a;
            tableB = b;
            return a[step << SHIFT];
        }
    }

    static String solve() {
        Solver solver = new Solver();
        long sum = 0;
        int m2 = -1;
        int m7 = -1;
        int m20 = -1;

        for (int n = 1; n <= 1000; n++) {
            int m = solver.computeM(n);
            if (n == 2)
                m2 = m;
            if (n == 7)
                m7 = m;
            if (n == 20)
                m20 = m;

            sum += (long) m * m * m;

            if (n == 20) {
                if (m2 != 2 || m7 != 1 || m20 != 4 || sum != 8150) {
                    throw new RuntimeException("Validation failed.");
                }
            }
        }
        return Long.toString(sum);
    }

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