Problem 649: Low-Prime Chessboard Nim

Project Euler Problem 649 Solution

EulerSolve provides an optimized solution for Project Euler Problem 649, Low-Prime Chessboard Nim, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Consider an $n \times n$ board, indexed by coordinates $0,1,\dots,n-1$ on each axis. A move decreases exactly one coordinate by one of the low primes $2,3,5,7$, provided the result stays nonnegative. With $c$ independent tokens, the whole position is an impartial normal-play sum, so it is losing exactly when the XOR of the token Grundy values is $0$. The task is to count winning starting placements modulo $10^9$. Mathematical Approach The direct search space has size $(n^2)^c$, which is hopeless for the actual parameters. The key simplification is to solve the one-coordinate game first, then lift that information to one square, and finally combine $c$ independent squares by XOR convolution. Step 1: Solve the one-dimensional subtraction game Let $s(i)$ be the Sprague-Grundy value of a single coordinate at position $i$. Legal moves go to $i-2$, $i-3$, $i-5$, and $i-7$ whenever those positions exist, so $$s(i)=\operatorname{mex}\{s(i-2),s(i-3),s(i-5),s(i-7)\},$$ where invalid predecessors are ignored. The base cases are immediate: $$s(0)=s(1)=0.$$ The next values are obtained mechanically from the recurrence: $$0,0,1,1,2,2,3,3,4,0,0,1,1,\dots$$ Only four follower values are ever inspected, so the mex can never exceed $4$....

Detailed mathematical approach

Problem Summary

Consider an $n \times n$ board, indexed by coordinates $0,1,\dots,n-1$ on each axis. A move decreases exactly one coordinate by one of the low primes $2,3,5,7$, provided the result stays nonnegative. With $c$ independent tokens, the whole position is an impartial normal-play sum, so it is losing exactly when the XOR of the token Grundy values is $0$. The task is to count winning starting placements modulo $10^9$.

Mathematical Approach

The direct search space has size $(n^2)^c$, which is hopeless for the actual parameters. The key simplification is to solve the one-coordinate game first, then lift that information to one square, and finally combine $c$ independent squares by XOR convolution.

Step 1: Solve the one-dimensional subtraction game

Let $s(i)$ be the Sprague-Grundy value of a single coordinate at position $i$. Legal moves go to $i-2$, $i-3$, $i-5$, and $i-7$ whenever those positions exist, so

$$s(i)=\operatorname{mex}\{s(i-2),s(i-3),s(i-5),s(i-7)\},$$

where invalid predecessors are ignored. The base cases are immediate:

$$s(0)=s(1)=0.$$

The next values are obtained mechanically from the recurrence:

$$0,0,1,1,2,2,3,3,4,0,0,1,1,\dots$$

Only four follower values are ever inspected, so the mex can never exceed $4$. Therefore

$$s(i)\in\{0,1,2,3,4\}\qquad (\forall i\ge 0).$$

Step 2: Lift one coordinate pair to one board square

A token on square $(x,y)$ is the disjunctive sum of two independent coordinate games, one on $x$ and one on $y$. By the Sprague-Grundy theorem, the square value is therefore

$$G(x,y)=s(x)\oplus s(y).$$

Now define the one-dimensional frequency table

$$A_r=\#\{\,i\in\{0,\dots,n-1\}:s(i)=r\,\}.$$

Then the number of board squares with value $t$ is

$$B_t=\sum_{a\oplus b=t} A_aA_b.$$

This is just the XOR convolution of the coordinate distribution with itself, and it satisfies

$$\sum_{t=0}^7 B_t=n^2.$$

Step 3: Combine $c$ tokens by XOR convolution

Let $D_j(u)$ denote the number of ways to place $j$ tokens so that the XOR of their square values is $u$. The empty placement has XOR $0$, so

$$D_0(0)=1,\qquad D_0(u)=0\quad (u\ne 0).$$

When one more token is added on a square of value $t$, the total XOR changes from $v$ to $v\oplus t$. Hence

$$D_{j+1}(u)=\sum_{t=0}^7 D_j(u\oplus t)\,B_t.$$

After $c$ steps, the losing placements are exactly those with zero nim-sum:

$$L=D_c(0).$$

Step 4: Convert losing placements into winning placements

Each token independently chooses one of the $n^2$ squares, so the total number of placements is

$$T=(n^2)^c.$$

Therefore the required answer is

$$W(n,c)\equiv T-L\equiv (n^2)^c-D_c(0)\pmod{10^9}.$$

The entire problem is now reduced to a small fixed-state dynamic program after the initial one-dimensional scan.

Step 5: Why the state space stays tiny

From Step 1, each coordinate Grundy value lies in $\{0,1,2,3,4\}$. Therefore a square value $s(x)\oplus s(y)$ can only lie in

$$\{0,1,2,3,4,5,6,7\}.$$

So both $B_t$ and $D_j(u)$ live in only eight states, regardless of how large $n$ becomes. That is why the huge board size does not produce a huge state space.

Worked Example: $n=3,\ c=2$

For coordinates $0,1,2$, the one-dimensional values are

$$s(0)=0,\qquad s(1)=0,\qquad s(2)=1,$$

so the frequency table is

$$A_0=2,\qquad A_1=1.$$

The board distribution becomes

$$B_0=A_0^2+A_1^2=2^2+1^2=5,$$

$$B_1=A_0A_1+A_1A_0=2\cdot 1+1\cdot 2=4,$$

and all other $B_t$ are $0$. For one token this already gives

$$D_1(0)=5,\qquad D_1(1)=4.$$

For two tokens, the zero-XOR count is

$$D_2(0)=D_1(0)B_0+D_1(1)B_1=5\cdot 5+4\cdot 4=41.$$

The total number of placements is

$$T=(3^2)^2=81,$$

so the winning count is

$$W(3,2)=81-41=40,$$

which matches the implementation checkpoint.

How the Code Works

The C++, Python, and Java implementations first scan the coordinates $0,1,\dots,n-1$ and compute the one-dimensional Grundy values directly from the mex recurrence. Because the recurrence only looks back $2$, $3$, $5$, and $7$ steps, the implementation keeps a rolling window of recent values instead of storing the whole sequence. During the same pass it accumulates the frequency table $A_r$.

Next, the implementation forms the board distribution $B_t$ by combining every pair of one-dimensional values with XOR. This is a tiny constant-size computation because only eight square states are relevant. It then performs a $c$-step dynamic program on those eight XOR states, reducing after every multiplication and addition modulo $10^9$.

Finally, the implementation computes $(n^2)^c \bmod 10^9$ with fast modular exponentiation, subtracts the zero-XOR count, and formats the result as a nine-digit string. The small consistency checks $W(3,1)=4$, $W(3,2)=40$, and $W(9,3)=450304$ agree with the derivation above.

Complexity Analysis

Computing the one-dimensional Grundy sequence for $n$ coordinates takes $O(n)$ time. Building the board distribution is constant work because the set of square values has fixed size $8$. The $c$-token dynamic program also runs over those same eight states, so it costs $O(c)$ time with a small constant factor. Fast modular exponentiation adds $O(\log c)$. For this fixed move set, the overall method is $O(n+c)$ time and $O(1)$ auxiliary memory.

Footnotes and References

Problem page: https://projecteuler.net/problem=649
Sprague-Grundy theorem: Wikipedia - Sprague-Grundy theorem
Nim: Wikipedia - Nim
mex (minimum excluded value): Wikipedia - mex
Combinatorial game theory: Wikipedia - Combinatorial game theory

Problem 649 source code

C++

#include <algorithm>
#include <cassert>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <vector>

namespace {

using u64 = std::uint64_t;
using u128 = __uint128_t;

constexpr u64 kMod = 1'000'000'000ULL;

u64 mod_mul(u64 a, u64 b) { return static_cast<u64>((u128)a * b % (u128)kMod); }

u64 mod_pow(u64 a, u64 e) {
    u64 r = 1 % kMod;
    while (e) {
        if (e & 1) r = mod_mul(r, a);
        a = mod_mul(a, a);
        e >>= 1;
    }
    return r;
}

std::vector<u64> sg_counts_1d(const int n) {
    std::vector<u64> cnt(8, 0);
    std::vector<std::uint8_t> buf(8, 0);
    const int moves[4] = {2, 3, 5, 7};
    for (int i = 0; i < n; ++i) {
        std::uint8_t g = 0;
        if (i >= 2) {
            bool seen[6] = {false, false, false, false, false, false};
            for (int d : moves) {
                if (i >= d) seen[buf[static_cast<std::size_t>((i - d) & 7)]] = true;
            }
            while (g < 6 && seen[g]) ++g;
        }
        buf[static_cast<std::size_t>(i & 7)] = g;
        ++cnt[g];
    }
    return cnt;
}

u64 M(const int n, const int c) {
    const auto cnt1 = sg_counts_1d(n);
    const int SZ = 8;
    std::vector<u64> cnt2(SZ, 0);
    for (int a = 0; a < SZ; ++a) {
        if (cnt1[a] == 0) continue;
        for (int b = 0; b < SZ; ++b) {
            if (cnt1[b] == 0) continue;
            const int g = a ^ b;
            cnt2[g] += cnt1[a] * cnt1[b];
        }
    }
    std::vector<u64> cnt2m(SZ, 0);
    for (int i = 0; i < SZ; ++i) cnt2m[i] = cnt2[i] % kMod;

    std::vector<u64> dp(SZ, 0), ndp(SZ, 0);
    dp[0] = 1;
    for (int t = 0; t < c; ++t) {
        std::fill(ndp.begin(), ndp.end(), 0);
        for (int x = 0; x < SZ; ++x) {
            const u64 vx = dp[x];
            if (vx == 0) continue;
            for (int g = 0; g < SZ; ++g) {
                ndp[x ^ g] += mod_mul(vx, cnt2m[g]);
                ndp[x ^ g] %= kMod;
            }
        }
        dp.swap(ndp);
    }

    const u64 n2m = mod_mul(static_cast<u64>(n) % kMod, static_cast<u64>(n) % kMod);
    const u64 total = mod_pow(n2m, static_cast<u64>(c));
    return (total + kMod - dp[0]) % kMod;
}

}  // namespace

int main() {
    assert(M(3, 1) == 4);
    assert(M(3, 2) == 40);
    assert(M(9, 3) == 450304);

    const u64 ans = M(10'000'019, 100);
    std::cout << std::setw(9) << std::setfill('0') << ans << "\n";
    return 0;
}

Python

def solve():
    n = 10000019
    c = 100
    MOD = 1000000000

    cnt = [0] * 8
    buf = [0] * 8
    moves = [2, 3, 5, 7]
    
    for i in range(n):
        g = 0
        if i >= 2:
            seen = [False] * 6
            for d in moves:
                if i >= d:
                    seen[buf[(i - d) & 7]] = True
            while g < 6 and seen[g]:
                g += 1
        buf[i & 7] = g
        cnt[g] += 1

    cnt2 = [0] * 8
    for a in range(8):
        if cnt[a] == 0: continue
        for b in range(8):
            if cnt[b] == 0: continue
            cnt2[a ^ b] += cnt[a] * cnt[b]
            
    cnt2m = [x % MOD for x in cnt2]
    
    dp = [0] * 8
    dp[0] = 1
    
    for _ in range(c):
        ndp = [0] * 8
        for x in range(8):
            vx = dp[x]
            if vx == 0: continue
            for g in range(8):
                ndp[x ^ g] = (ndp[x ^ g] + vx * cnt2m[g]) % MOD
        dp = ndp
        
    n2m = ((n % MOD) * (n % MOD)) % MOD
    total = pow(n2m, c, MOD)
    ans = (total + MOD - dp[0]) % MOD
    return "{:09d}".format(ans)

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

Java

public class Euler649 {
    static final long kMod = 1000000000L;

    static long modPow(long a, long e) {
        long r = 1 % kMod;
        while (e > 0) {
            if ((e & 1) != 0)
                r = (r * a) % kMod;
            a = (a * a) % kMod;
            e >>= 1;
        }
        return r;
    }

    public static String solve() {
        int n = 10000019;
        int c = 100;

        long[] cnt = new long[8];
        int[] buf = new int[8];
        int[] moves = { 2, 3, 5, 7 };

        for (int i = 0; i < n; ++i) {
            int g = 0;
            if (i >= 2) {
                boolean[] seen = new boolean[6];
                for (int d : moves) {
                    if (i >= d)
                        seen[buf[(i - d) & 7]] = true;
                }
                while (g < 6 && seen[g])
                    ++g;
            }
            buf[i & 7] = g;
            cnt[g]++;
        }

        long[] cnt2 = new long[8];
        for (int a = 0; a < 8; ++a) {
            if (cnt[a] == 0)
                continue;
            for (int b = 0; b < 8; ++b) {
                if (cnt[b] == 0)
                    continue;
                cnt2[a ^ b] += cnt[a] * cnt[b];
            }
        }

        long[] cnt2m = new long[8];
        for (int i = 0; i < 8; ++i)
            cnt2m[i] = cnt2[i] % kMod;

        long[] dp = new long[8];
        long[] ndp = new long[8];
        dp[0] = 1;

        for (int t = 0; t < c; ++t) {
            for (int i = 0; i < 8; i++)
                ndp[i] = 0;
            for (int x = 0; x < 8; ++x) {
                long vx = dp[x];
                if (vx == 0)
                    continue;
                for (int g = 0; g < 8; ++g) {
                    ndp[x ^ g] = (ndp[x ^ g] + vx * cnt2m[g]) % kMod;
                }
            }
            System.arraycopy(ndp, 0, dp, 0, 8);
        }

        long n2m = ((long) (n % kMod) * (long) (n % kMod)) % kMod;
        long total = modPow(n2m, c);
        long ans = (total + kMod - dp[0]) % kMod;

        return String.format("%09d", ans);
    }

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