Problem 310: Nim Square

Project Euler Problem 310 Solution

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

Problem Summary In Nim Square, a move on one heap of size $n$ removes any positive square $k^2\le n$. The position has three heaps, and we must count how many ordered-by-size triples $$0\le a\le b\le c\le N$$ are losing positions. Mathematical Approach 1) Single-heap Grundy numbers For one heap, the legal moves are $$n\to n-1,\ n-4,\ n-9,\dots$$ for every square not exceeding $n$. So the Grundy value of a single heap is $$g(0)=0,$$ $$g(n)=\operatorname{mex}\{g(n-k^2):k^2\le n\}.$$ This is the usual subtract-a-square recurrence. 2) Why three heaps reduce to xor The full game is the disjoint sum of three impartial subgames. Therefore Sprague-Grundy theory says that the triple $(a,b,c)$ is losing exactly when the nim-sum vanishes: $$g(a)\oplus g(b)\oplus g(c)=0.$$ So once the array $g(0),\dots,g(N)$ is known, the entire counting problem becomes a combinatorial xor-counting problem. 3) Small worked example The first few single-heap values are $$g(0..4)=(0,1,0,1,2).$$ For instance: $$g(1)=\operatorname{mex}\{g(0)\}=\operatorname{mex}\{0\}=1,$$ $$g(2)=\operatorname{mex}\{g(1)\}=\operatorname{mex}\{1\}=0,$$ $$g(4)=\operatorname{mex}\{g(3),g(0)\}=\operatorname{mex}\{1,0\}=2.$$ If $N=4$, the losing triples are $$\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}$$ so the count is $12$....

Detailed mathematical approach

Problem Summary

In Nim Square, a move on one heap of size $n$ removes any positive square $k^2\le n$. The position has three heaps, and we must count how many ordered-by-size triples

$$0\le a\le b\le c\le N$$

are losing positions.

Mathematical Approach

1) Single-heap Grundy numbers

For one heap, the legal moves are

$$n\to n-1,\ n-4,\ n-9,\dots$$

for every square not exceeding $n$. So the Grundy value of a single heap is

$$g(0)=0,$$

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

This is the usual subtract-a-square recurrence.

2) Why three heaps reduce to xor

The full game is the disjoint sum of three impartial subgames. Therefore Sprague-Grundy theory says that the triple $(a,b,c)$ is losing exactly when the nim-sum vanishes:

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

So once the array $g(0),\dots,g(N)$ is known, the entire counting problem becomes a combinatorial xor-counting problem.

3) Small worked example

The first few single-heap values are

$$g(0..4)=(0,1,0,1,2).$$

For instance:

$$g(1)=\operatorname{mex}\{g(0)\}=\operatorname{mex}\{0\}=1,$$

$$g(2)=\operatorname{mex}\{g(1)\}=\operatorname{mex}\{1\}=0,$$

$$g(4)=\operatorname{mex}\{g(3),g(0)\}=\operatorname{mex}\{1,0\}=2.$$

If $N=4$, the losing triples are

$$\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}$$

so the count is $12$.

4) The naive count is cubic

A direct scan over all triples $(a,b,c)$ would take

$$O(N^3)$$

time, which is impossible for $N=100000$. The code instead reorganizes the count by fixing the middle heap $b$.

5) Prefix/suffix counting for a fixed middle heap

Fix $b$. Then:

Left side: $a$ may range over $0\le a\le b$.

Right side: $c$ may range over $b\le c\le N$.

Let

$$\text{pref}[x]=\#\{a\le b:g(a)=x\},$$

$$\text{suf}[y]=\#\{c\ge b:g(c)=y\}.$$

For a fixed Grundy value $x=g(a)$, the losing condition

$$x\oplus g(b)\oplus g(c)=0$$

forces

$$g(c)=x\oplus g(b).$$

Hence the number of losing triples with this middle value $b$ is

$$\text{add}_b=\sum_x \text{pref}[x]\cdot \text{suf}[x\oplus g(b)].$$

Summing $\text{add}_b$ over all $b$ gives the final answer.

6) Why this automatically respects $a\le b\le c$

The implementation maintains the arrays in a very specific order:

Before counting for $b$: it adds $g(b)$ to the prefix array.

During counting: prefix therefore represents $a\in[0,b]$, while suffix still represents $c\in[b,N]$.

After counting: it removes $g(b)$ from suffix, so the next step uses $c\ge b+1$.

This means the ordering constraint is built directly into the data structure update order; no extra combinatorial correction is needed.

How the Code Works

1) Build Grundy table. build_grundy(limit) computes the subtract-a-square Grundy values using a timestamp-based mex array.

2) Determine Grundy alphabet size. The code records

$$\max g(n)$$

so it knows how large the prefix and suffix frequency arrays must be.

3) Initialize suffix. Initially all heaps $0,\dots,N$ are available on the right, so suffix counts every Grundy value in the full range.

4) Sweep the middle index. For each $b$, the code updates prefix, applies the xor-count formula above, then removes $b$ from suffix.

5) Checkpoints. The C++ version verifies

$$\texttt{solve(29)}=1160$$

and also checks that $\texttt{solve(40)}$ matches a literal brute-force triple scan.

6) Final Project Euler value. For the required bound $N=100000$, the implementation returns

$$2586528661783.$$

Complexity Analysis

Computing the Grundy table costs about $O(N\sqrt N)$, because each $g(n)$ examines all squares up to $n$. The counting sweep costs $O(N\cdot G)$, where $G$ is the number of distinct Grundy values encountered, which stays small in practice. Memory usage is $O(N+G)$.

Problem 310 source code

C++

#include <algorithm>
#include <cstdint>
#include <iostream>
#include <string>
#include <vector>

namespace {

using i64 = long long;

struct Options {
    int limit = 100000;
    bool run_checkpoints = true;
};

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }
    int parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<int>(c - '0');
    }
    value = parsed;
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_int_after_prefix(arg, "--limit=", options.limit)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.limit >= 0;
}

std::vector<int> build_grundy(const int limit) {
    std::vector<int> g(static_cast<std::size_t>(limit + 1), 0);
    std::vector<int> seen(512, -1);

    for (int n = 1; n <= limit; ++n) {
        for (int k = 1; k * k <= n; ++k) {
            const int val = g[static_cast<std::size_t>(n - k * k)];
            if (val >= static_cast<int>(seen.size())) {
                seen.resize(static_cast<std::size_t>(val + 64), -1);
            }
            seen[static_cast<std::size_t>(val)] = n;
        }
        while (seen[static_cast<std::size_t>(g[static_cast<std::size_t>(n)])] == n) {
            ++g[static_cast<std::size_t>(n)];
        }
    }

    return g;
}

i64 solve(const int limit) {
    const std::vector<int> g = build_grundy(limit);
    int max_g = 0;
    for (int x : g) {
        max_g = std::max(max_g, x);
    }

    std::vector<i64> prefix(static_cast<std::size_t>(max_g + 1), 0);
    std::vector<i64> suffix(static_cast<std::size_t>(max_g + 1), 0);
    for (int n = 0; n <= limit; ++n) {
        ++suffix[static_cast<std::size_t>(g[static_cast<std::size_t>(n)])];
    }

    i64 answer = 0;
    for (int b = 0; b <= limit; ++b) {
        const int gb = g[static_cast<std::size_t>(b)];
        ++prefix[static_cast<std::size_t>(gb)];

        for (int ga = 0; ga <= max_g; ++ga) {
            const int gc = ga ^ gb;
            if (gc > max_g) {
                continue;
            }
            answer += prefix[static_cast<std::size_t>(ga)] * suffix[static_cast<std::size_t>(gc)];
        }

        --suffix[static_cast<std::size_t>(gb)];
    }

    return answer;
}

i64 brute_small(const int limit) {
    const std::vector<int> g = build_grundy(limit);
    i64 count = 0;
    for (int a = 0; a <= limit; ++a) {
        for (int b = a; b <= limit; ++b) {
            for (int c = b; c <= limit; ++c) {
                if ((g[static_cast<std::size_t>(a)] ^ g[static_cast<std::size_t>(b)] ^ g[static_cast<std::size_t>(c)]) == 0) {
                    ++count;
                }
            }
        }
    }
    return count;
}

bool run_checkpoints() {
    if (solve(29) != 1160LL) {
        std::cerr << "Checkpoint failed for limit 29" << '\n';
        return false;
    }
    if (solve(40) != brute_small(40)) {
        std::cerr << "Checkpoint failed for brute cross-check at limit 40" << '\n';
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }
    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }

    std::cout << solve(options.limit) << '\n';
    return 0;
}

Python

def build_grundy(limit):
    g = [0] * (limit + 1)
    seen = [-1] * 512
    
    for n in range(1, limit + 1):
        k = 1
        while k * k <= n:
            val = g[n - k * k]
            if val >= len(seen):
                seen.extend([-1] * (val + 64))
            seen[val] = n
            k += 1
            
        while seen[g[n]] == n:
            g[n] += 1
            
    return g

def solve(limit=100000):
    g = build_grundy(limit)
    max_g = max(g)
    
    prefix = [0] * (max_g + 1)
    suffix = [0] * (max_g + 1)
    
    for n in range(limit + 1):
        suffix[g[n]] += 1
        
    answer = 0
    for b in range(limit + 1):
        gb = g[b]
        prefix[gb] += 1
        
        for ga in range(max_g + 1):
            gc = ga ^ gb
            if gc > max_g:
                continue
            answer += prefix[ga] * suffix[gc]
            
        suffix[gb] -= 1
        
    return str(answer)

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

Java

import java.util.Arrays;

public class Euler310 {
    static int[] buildGrundy(int limit) {
        int[] g = new int[limit + 1];
        int[] seen = new int[512];
        Arrays.fill(seen, -1);

        for (int n = 1; n <= limit; ++n) {
            for (int k = 1; k * k <= n; ++k) {
                int val = g[n - k * k];
                if (val >= seen.length) {
                    int[] newSeen = new int[val + 64];
                    System.arraycopy(seen, 0, newSeen, 0, seen.length);
                    Arrays.fill(newSeen, seen.length, newSeen.length, -1);
                    seen = newSeen;
                }
                seen[val] = n;
            }
            while (seen[g[n]] == n) {
                ++g[n];
            }
        }
        return g;
    }

    public static String solve() {
        int limit = 100000;
        int[] g = buildGrundy(limit);
        int maxG = 0;
        for (int x : g) {
            if (x > maxG) {
                maxG = x;
            }
        }

        long[] prefix = new long[maxG + 1];
        long[] suffix = new long[maxG + 1];

        for (int n = 0; n <= limit; ++n) {
            ++suffix[g[n]];
        }

        long answer = 0;
        for (int b = 0; b <= limit; ++b) {
            int gb = g[b];
            ++prefix[gb];

            for (int ga = 0; ga <= maxG; ++ga) {
                int gc = ga ^ gb;
                if (gc > maxG) {
                    continue;
                }
                answer += prefix[ga] * suffix[gc];
            }

            --suffix[gb];
        }

        return String.valueOf(answer);
    }

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