Problem 260: Stone Game

Project Euler Problem 260 Solution

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

Problem Summary We work with three heap sizes, written canonically as $$0\le x\le y\le z.$$ A legal move chooses a positive integer $n$ and subtracts that same $n$ from exactly one heap, from any pair of heaps, or from all three heaps, provided no coordinate becomes negative. A position is losing (a P-position) if there is no legal move from it to another losing position. The task is to sum $$x+y+z$$ over all losing triples with coordinates at most a fixed limit. Mathematical Approach 1. Normal-Play Interpretation This is a finite impartial normal-play game. Every move strictly decreases at least one coordinate, so the game graph is acyclic. Therefore every position is either: 1. losing: no move reaches a losing state, 2. winning: at least one move reaches a losing state. The whole problem is therefore to characterize the losing triples efficiently. 2. Canonical Ordering and Why a Forward Scan Works The code stores each position in sorted order $x\le y\le z$. After any legal move we sort again, so each game state appears exactly once. If we move from $(x,y,z)$ to a successor and then sort, the new triple is always lexicographically smaller than $(x,y,z)$: at least one coordinate decreases, and no coordinate increases....

Detailed mathematical approach

Problem Summary

Two players take turns removing stones from three piles. A move chooses a number $N\gt0$ and removes $N$ stones from one pile, or $N$ stones from each of two piles, or $N$ stones from each of the three piles; the player who takes the last stone or stones wins. A configuration is winning if the player about to move can force a win, and losing if the other player can force a win whatever the first player does. For example, $(0,0,13)$, $(0,11,11)$ and $(5,5,5)$ are winning, since the player to move takes everything at once, while $(0,1,2)$ and $(1,3,3)$ are losing. Over the losing configurations $(x,y,z)$ with $x\le y\le z\le100$, the sum of $x+y+z$ is $173895$; we must find this sum over the losing configurations with $x\le y\le z\le1000$. ^[1] With one pile empty, the game is Wythoff's classical two-pile game. ^[2]

The whole solution in one sentence: each kind of move keeps a pair of numbers fixed (two pile sizes, or one pile size and the difference of the other two, or the two gaps between the sorted piles), so we scan the configurations from small to large and call one losing exactly when none of its seven such pairs belongs to an earlier losing configuration, using at most seven table look-ups.

Mathematical Approach

We classify positions backwards from the end of the game, order them so that every position comes after the positions it can move to, and replace the moves by seven pairs of numbers.

1. Winning and losing positions

Every move removes at least one stone, so a game that starts at $(x,y,z)$ ends after at most $x+y+z$ moves, without draws. Such a game is solved backwards from its end, by backward induction: a position is losing exactly when every move leads to a winning position, and winning exactly when some move leads to a losing position. ^[3] The empty position $(0,0,0)$ is losing, because the player to move has no move left. Losing positions are also called P-positions, since the previous player wins, and winning positions N-positions, since the next player wins.

From $(0,1,2)$ the moves lead to $(0,0,2)$, $(0,1,1)$ and $(0,0,1)$, and the opponent empties each of these in one move, so $(0,1,2)$ is losing. Similarly, each of the $9$ positions reachable from $(1,3,3)$ has a move to $(0,0,0)$ or to $(0,1,2)$; for example, $(1,1,3)$ becomes $(0,1,2)$ when one stone is taken from a pile of $1$ and one from the pile of $3$. So $(1,3,3)$ is losing too.

2. Sorted triples and the order of the scan

The piles are interchangeable, so we describe a position by its pile sizes from smallest to largest, a sorted triple $(x,y,z)$ with $x\le y\le z$, and compare sorted triples in lexicographic order, as in a dictionary: by $x$ first, then by $y$, then by $z$. ^[4] For example, $(0,4,4)$ comes before $(1,1,1)$.

Claim: every move turns a sorted triple into a lexicographically smaller one. For $k=1,2,3$, the $k$ piles that held the $k$ smallest sizes still hold at most the old $k$-th smallest size after the move, so the new $k$-th smallest size is at most the old one. No entry grows while the total drops, so at the first entry where the two triples differ, the new one is smaller.

Hence three nested loops over $x\le y\le z$, with $x$ outermost and $z$ innermost, visit every position after all positions it can move to. Those positions have no pile larger than $z$, so a scan of the triples with $z\le L$ never leaves its range. There are $\binom{L+3}{3}$ such triples, because $(x,y,z)\mapsto(x,y+1,z+2)$ maps them one-to-one onto the sets of three different numbers from $0,1,\dots,L+2$, and subtracting $0$, $1$ and $2$ undoes it; this is a case of stars and bars. ^[5] There are $35$ triples for $L=4$ and $167\,668\,501$ for $L=1000$.

3. Seven keys

For a sorted triple $(x,y,z)$, each of the seven kinds of moves keeps a pair of numbers fixed, which we call its key; the index of $K$ is the number of piles that the move touches:

one pile: the other two piles stay, which gives the keys $K_1(y,z)$, $K_1(x,z)$ and $K_1(x,y)$ for moves on $x$, on $y$ and on $z$;
two piles: the third pile and the difference of the two reduced piles stay, which gives $K_2(y-x,z)$, $K_2(z-y,x)$ and $K_2(z-x,y)$ for moves on $x,y$, on $y,z$ and on $x,z$;
all three piles: the gaps stay, which gives $K_3(y-x,z-x)$.

A key is a property that does not depend on the order of the piles, so it survives a move that reorders them: $K_1(a,b)$ with $a\le b$ says that two piles have the sizes $a$ and $b$; $K_2(d,c)$ says that one pile has size $c$ and the other two differ by $d$; $K_3(d,e)$ says that the sorted piles are $t$, $t+d$ and $t+e$ for some $t$. For example, $(1,3,3)$ has the five different keys $K_1(3,3)$, $K_1(1,3)$, $K_2(2,3)$, $K_2(0,1)$ and $K_3(2,2)$; the second and third occur twice, because two of its piles are equal.

4. Two losing positions never share a key

Claim: if two different positions share a key, the one that comes later in lexicographic order can move to the other.

Shared piles $a$ and $b$: the third piles differ, say $c\gt c'$, and taking $c-c'$ stones from the pile $c$ turns the first position into the second.
A shared pile $c$, with other piles $p,p+d$ and $p',p'+d$, say $p\gt p'$: taking $p-p'$ stones from each of the piles $p$ and $p+d$ turns the first into the second.
Shared gaps, $(t,t+d,t+e)$ and $(t',t'+d,t'+e)$ with $t\gt t'$: taking $t-t'$ stones from every pile turns the first into the second.

By Section 2, the position that moves is the later one.

Theorem. A sorted triple $T$ is losing exactly when no losing triple that comes before $T$ has one of the seven keys of $T$. If one does, $T$ can move to it by the claim, so $T$ is winning. If none does, no move of $T$ reaches a losing position, because such a position would come before $T$ by Section 2 and would keep the key of the move; so $T$ is losing.

In particular two losing positions never share a key; for example, two given sizes $a\le b$ occur together in at most one losing position.

5. The scan

The theorem gives the algorithm: three tables of marks, one for each kind of key, all empty at the start, and a scan of the triples with $z\le L$ in lexicographic order. If a key of the current triple is marked, the triple is winning; otherwise it is losing, $x+y+z$ is added to the sum, and its seven keys are marked. By induction on the scan order, the marks are exactly the keys of the losing triples visited so far, and these suffice: the losing triple reached by a winning move has no pile larger than $z$, so it has been visited. The entries of a key lie between $0$ and $L$, so each table has $(L+1)^2$ cells, and a triple costs at most seven look-ups instead of a search through all of its moves.

6. Worked example: at most four stones per pile

For $L=4$ the scan visits $35$ triples. $(0,0,0)$ is losing and marks $K_1(0,0)$, $K_2(0,0)$ and $K_3(0,0)$. Then $(0,0,1)$ to $(0,0,4)$ are winning through $K_1(0,0)$, by emptying the largest pile, and $(0,1,1)$ through $K_2(0,0)$, by taking one stone from each non-empty pile. $(0,1,2)$ has no marked key and is losing; it marks $K_1(1,2)$, $K_1(0,2)$, $K_1(0,1)$, $K_2(1,2)$, $K_2(1,0)$, $K_2(2,1)$ and $K_3(1,2)$. Further on:

$(0,3,4)$ is winning through $K_2(1,0)$: beside the empty pile, its piles $3$ and $4$ differ by $1$, and taking $2$ stones from each gives $(0,1,2)$;
$(1,1,3)$ is winning through $K_2(2,1)$: beside a pile $1$, its piles $1$ and $3$ differ by $2$, like the piles $0$ and $2$ of $(0,1,2)$, and taking one stone from each gives $(0,1,2)$;
$(1,1,4)$ is losing and marks $K_1(1,4)$, $K_1(1,1)$, $K_2(0,4)$, $K_2(3,1)$ and $K_3(0,3)$;
$(1,3,3)$ is losing, since none of its five keys from Section 3 is marked;
$(2,2,4)$ is winning through $K_2(0,4)$, moving to $(1,1,4)$, and $(2,4,4)$ through $K_3(2,2)$, moving to $(1,3,3)$;
$(3,4,4)$ is the last losing triple.

So the losing triples with $z\le4$ are $(0,0,0)$, $(0,1,2)$, $(1,1,4)$, $(1,3,3)$ and $(3,4,4)$, with the sum $0+3+6+7+11=27$. For $L=100$ the scan gives the sum $173895$ of the statement.

7. One empty pile: Wythoff's game

If $x=0$, no move touches the empty pile, and the game on the other two piles is Wythoff's game. W. A. Wythoff, who analysed it in 1907, found its losing positions:

$$\bigl(\lfloor k\varphi\rfloor,\lfloor k\varphi^2\rfloor\bigr),\qquad\varphi=\frac{1+\sqrt5}{2},$$

for $k=0,1,2,\dots$, where $\varphi$ is the golden ratio and $\lfloor\cdot\rfloor$ rounds down: $(0,0)$, $(1,2)$, $(3,5)$, $(4,7)$, $(6,10)$, $(8,13)$ and so on. ^[2] The scan agrees: for $L=1000$ its losing triples with $x=0$ are exactly the $383$ triples $(0,\lfloor k\varphi\rfloor,\lfloor k\varphi^2\rfloor)$ with $0\le k\le382$. Two classical facts about these pairs reflect Section 4 in this slice. Every positive integer lies in exactly one pair, ^[2] so no two pairs share a key $K_1(0,b)$; and since $\varphi^2=\varphi+1$ and $k$ is an integer, $\lfloor k\varphi^2\rfloor=\lfloor k\varphi\rfloor+k$, so the differences $k$ of the pairs are all distinct, and no two pairs share a key $K_2(d,0)$.

8. Every pair of piles has exactly one partner

By Section 4, at most one losing position has two piles of given sizes $a\le b$; in fact exactly one has, so every pair has a unique partner $c$, which need not be the largest pile: by Section 6, the sizes $1$ and $1$ have the partner $4$, and the sizes $1$ and $4$ the partner $1$. To prove existence, fix two labelled piles of sizes $a$ and $b$, let the third pile have any size $c$, and call $c$ covered if some move that touches at least one of the two fixed piles reaches a losing position. For a given kind of move and a given $n$, that losing position contains the two piles $a-n$ and $b$, or $a$ and $b-n$, or $a-n$ and $b-n$, so by Section 4 it is unique, and it fixes $c$: its third pile is $c$, or $c-n$ if pile $c$ was reduced as well. The six kinds of such moves, on $a$, on $b$, on $a,b$, on $a,c$, on $b,c$ and on all three, allow at most $a$, $b$, $\min(a,b)$, $a$, $b$ and $\min(a,b)$ values of $n$, so at most $2a+2b+2\min(a,b)$ sizes are covered. For the smallest uncovered size $c^*$, the moves that reduce only the third pile reach covered, hence winning, positions, and the other moves reach no losing position; so the position with piles $a$, $b$, $c^*$ is losing, and $c^*\le2(a+b+\min(a,b))$.

By this uniqueness, a losing sorted triple is determined by its two smallest piles: at most $(L+1)(L+2)/2$ triples with $z\le L$ are losing, $501\,501$ for $L=1000$, and the required sum is at most $3L$ times this number, about $1.5\cdot10^9$. The scan in fact finds $1264$ losing triples for $L=100$ and $121\,060$ for $L=1000$, about $0.12\,L^2$, so for $L=1000$ only about one sorted triple in $1400$ is losing.

How the Code Works

The C++ program ^[6] reads --limit= (default $1000$), --check-limit= (default $35$) and --skip-checkpoints in parse_arguments, runs run_checkpoints unless told not to, and prints the result of solve_fast for the limit $L$.

solve_fast is the scan of Section 5. Its tables one, two and all hold the marks of $K_1$, $K_2$ and $K_3$ in zero-filled byte vectors of $(L+1)^2$ cells, and the lambda id(a, b) computes the cell $a(L+1)+b$. Three nested loops visit $x\le y\le z\le L$ in lexicographic order, and one condition tests the cells of the seven keys of Section 3, from one[id(y, z)] for $K_1(y,z)$ to all[id(y - x, z - x)] for $K_3(y-x,z-x)$. If none is set, $x+y+z$ is added to a 64-bit unsigned sum and the seven cells are set to $1$.

The checks use a brute force that follows Section 1 literally and knows nothing about keys. For the check limit $M$ it keeps one byte for each of the $(M+1)^3$ triples and visits the sorted triples layer by layer in order of their total, which every move lowers. has_move_to_losing_state tries the moves one after another until one reaches a losing triple, at most $4x+2y+z$ of them: $n\le x$, $n\le y$ or $n\le z$ stones from one pile, $n\le x$ from $x,y$ or from $x,z$, $n\le y$ from $y,z$, and $n\le x$ from all three. It sorts each result with sort_triplet before the look-up, and a triple without a move to a losing triple is losing. run_checkpoints refuses $M\lt13$, so that $(0,0,13)$ fits, and then requires that the brute-force sum equals solve_fast for $L=M$, that query_is_losing of the returned BruteResult finds $(0,1,2)$ and $(1,3,3)$ losing and $(0,0,13)$, $(0,11,11)$ and $(5,5,5)$ winning, and that solve_fast returns $173895$ for $L=100$. A failed check prints a message and ends the program with exit code $1$.

The Python solve is the same scan with bytearray tables, without options or checks; Python integers are unbounded. The Java solve uses byte arrays, a long sum and two shortcuts: it skips a pair $(x,y)$ whose key $K_1(x,y)$ is already marked, and it leaves the $z$-loop with break after the first losing triple. Both are safe: while $K_1(x,y)$ is marked, no triple $(x,y,z)$ can pass the test, and a losing triple found by the loop marks it; by Section 4, no second losing triple has the smallest piles $x$ and $y$ anyway. For $L=1000$ the sum is below $1.51\cdot10^9$ by Section 8, so it fits easily into the integer types of all three programs.

Complexity and Verification

The C++ and Python programs test all $\binom{L+3}{3}$ sorted triples, $167\,668\,501$ for $L=1000$, with at most seven look-ups each, so the time is $O(L^3)$; the tables take $3(L+1)^2$ bytes, $3\,006\,003$ for $L=1000$, so the memory is $O(L^2)$. The Java shortcuts skip $241\,370$ of the $501\,501$ pairs $(x,y)$ and run the inner loop $46\,529\,505$ times, about $3.6$ times fewer; a pair whose partner exceeds $L$ still runs its $z$-loop to the end, and the bound $O(L^3)$ still holds. On our test machine the C++ program, checks included, and the Java program finish in well under a second, the Python program in about $30$ seconds. The brute force of the checks needs $(M+1)^3$ bytes and up to $4x+2y+z$ moves per triple: negligible for $M=35$, with $8436$ sorted triples, but about $10^9$ bytes and up to $4.6\cdot10^{11}$ moves for $L=1000$, against at most $7\binom{1003}{3}\approx1.2\cdot10^9$ look-ups for the scan.

For this article we also ran checks that are not part of the published programs. An independent brute force in C, which tries every move of every position in order of the total, agrees with the scan on all $4\,590\,551$ sorted triples with $z\le300$, of which $10\,943$ are losing. A scan up to $L=4000$ finds exactly one partner, never above $1896$, for every pair $a\le b\le1000$, as Section 8 proves, and these partners reproduce both Java counts above exactly. With --check-limit=100 the brute force of the C++ program agrees with solve_fast, and the three published programs print the same number for $L=1000$ as an independent C implementation of the scan.

Footnotes and References

The references below supply the background facts used above; the seven keys, the order of the scan, the partner of a pair of piles and the operation counts are derived in this article.

Project Euler 260 — Stone Game. The official statement allows a player to remove N > 0 stones from one pile, or N stones from each of two or of all three piles; the player taking the last stones wins. It gives the sum 173895 over the losing configurations with x ≤ y ≤ z ≤ 100 and asks for the sum for x ≤ y ≤ z ≤ 1000.
Wikipedia — Wythoff's game. Players remove counters from one pile, or the same number from both piles, and the player who removes the last counters wins; W. A. Wythoff published an analysis in 1907. The k-th cold position, lost by the player to move, is (⌊kφ⌋, ⌊kφ²⌋) with the golden ratio φ, and each positive integer appears in exactly one of the two sequences.
T. S. Ferguson — Game Theory, Part I: Impartial Combinatorial Games. Sections 1.1 to 1.3 solve a take-away game by backward induction and characterise P-positions, won by the previous player, and N-positions, won by the next player, in games that always end under normal play: terminal positions are P-positions, every N-position has a move to a P-position, and every move from a P-position leads to an N-position.
Wikipedia — Lexicographic order. Words of the same length are ordered by their symbols at the first place where they differ, and the order extends to Cartesian products of ordered sets.
Wikipedia — Stars and bars (combinatorics). There are C(n + k − 1, n) multisets of size n with elements from a set of size k; n = 3 and k = L + 1 count the sorted triples.
C++ — Euler260.cpp. The linked immutable C++ revision contains solve_fast, has_move_to_losing_state and run_checkpoints. It is the implementation record for this article.

Problem 260 source code

C++

#include <algorithm>
#include <array>
#include <cstddef>
#include <cstdint>
#include <exception>
#include <iostream>
#include <limits>
#include <stdexcept>
#include <string>
#include <thread>
#include <vector>

namespace {

using u8 = std::uint8_t;
using u64 = std::uint64_t;

struct Options {
    int limit = 1000;
    bool run_checkpoints = true;
    int checkpoint_bruteforce_limit = 35;
    unsigned requested_threads = 0U;
};

bool parse_unsigned_after_prefix(const std::string& arg,
                                 const std::string& prefix,
                                 unsigned& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }

    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }

    u64 parsed = 0ULL;
    for (const char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        const u64 digit = static_cast<u64>(c - '0');
        if (parsed > (static_cast<u64>(std::numeric_limits<unsigned>::max()) - digit) /
                          10ULL) {
            return false;
        }
        parsed = parsed * 10ULL + digit;
    }

    value = static_cast<unsigned>(parsed);
    return 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;
    }

    std::int64_t parsed = 0;
    for (const char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<std::int64_t>(c - '0');
        if (parsed > static_cast<std::int64_t>(std::numeric_limits<int>::max())) {
            return false;
        }
    }

    value = static_cast<int>(parsed);
    return true;
}

bool parse_arguments(const 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;
        }
        if (parse_int_after_prefix(arg, "--check-limit=",
                                   options.checkpoint_bruteforce_limit)) {
            continue;
        }
        if (parse_unsigned_after_prefix(arg, "--threads=", options.requested_threads)) {
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }

    return true;
}

unsigned pick_thread_count(const unsigned requested) {
    if (requested > 0U) {
        return requested;
    }

    unsigned hw = std::thread::hardware_concurrency();
    if (hw == 0U) {
        hw = 4U;
    }
    return hw;
}

std::size_t checked_square_size(const int side) {
    if (side <= 0) {
        throw std::runtime_error("Side length must be positive");
    }

    const std::size_t s = static_cast<std::size_t>(side);
    if (s > std::numeric_limits<std::size_t>::max() / s) {
        throw std::overflow_error("2D table size overflow");
    }
    return s * s;
}

std::size_t checked_cube_size(const int side) {
    const std::size_t square = checked_square_size(side);
    const std::size_t s = static_cast<std::size_t>(side);
    if (square > std::numeric_limits<std::size_t>::max() / s) {
        throw std::overflow_error("3D table size overflow");
    }
    return square * s;
}

u64 solve_fast(const int limit) {
    if (limit < 0) {
        throw std::invalid_argument("Limit must be non-negative");
    }

    const int side = limit + 1;
    const std::size_t table_size = checked_square_size(side);

    std::vector<u8> one(table_size, 0U);
    std::vector<u8> two(table_size, 0U);
    std::vector<u8> all(table_size, 0U);

    const auto id = [side](const int a, const int b) -> std::size_t {
        return static_cast<std::size_t>(a) * static_cast<std::size_t>(side) +
               static_cast<std::size_t>(b);
    };

    u64 sum = 0ULL;

    for (int x = 0; x <= limit; ++x) {
        for (int y = x; y <= limit; ++y) {
            for (int z = y; z <= limit; ++z) {
                if (one[id(y, z)] != 0U || one[id(x, z)] != 0U ||
                    one[id(x, y)] != 0U || two[id(y - x, z)] != 0U ||
                    two[id(z - y, x)] != 0U || two[id(z - x, y)] != 0U ||
                    all[id(y - x, z - x)] != 0U) {
                    continue;
                }

                sum += static_cast<u64>(x + y + z);

                one[id(y, z)] = 1U;
                one[id(x, z)] = 1U;
                one[id(x, y)] = 1U;

                two[id(y - x, z)] = 1U;
                two[id(z - y, x)] = 1U;
                two[id(z - x, y)] = 1U;

                all[id(y - x, z - x)] = 1U;
            }
        }
    }

    return sum;
}

inline void sort_triplet(int& a, int& b, int& c) {
    if (a > b) {
        std::swap(a, b);
    }
    if (b > c) {
        std::swap(b, c);
    }
    if (a > b) {
        std::swap(a, b);
    }
}

struct BruteResult {
    int limit = 0;
    int side = 0;
    u64 losing_sum = 0ULL;
    std::vector<u8> is_losing;

    bool query_is_losing(int x, int y, int z) const {
        sort_triplet(x, y, z);
        const std::size_t index =
            (static_cast<std::size_t>(x) * static_cast<std::size_t>(side) +
             static_cast<std::size_t>(y)) *
                static_cast<std::size_t>(side) +
            static_cast<std::size_t>(z);
        return is_losing[index] != 0U;
    }
};

bool has_move_to_losing_state(const std::vector<u8>& is_losing,
                              const int side,
                              const int x,
                              const int y,
                              const int z) {
    const auto index3 = [side](const int a, const int b, const int c) -> std::size_t {
        return (static_cast<std::size_t>(a) * static_cast<std::size_t>(side) +
                static_cast<std::size_t>(b)) *
                   static_cast<std::size_t>(side) +
               static_cast<std::size_t>(c);
    };

    const auto can_reach_losing = [&is_losing, &index3](int a, int b, int c) {
        sort_triplet(a, b, c);
        return is_losing[index3(a, b, c)] != 0U;
    };

    for (int n = 1; n <= x; ++n) {
        if (can_reach_losing(x - n, y, z)) {
            return true;
        }
    }
    for (int n = 1; n <= y; ++n) {
        if (can_reach_losing(x, y - n, z)) {
            return true;
        }
    }
    for (int n = 1; n <= z; ++n) {
        if (can_reach_losing(x, y, z - n)) {
            return true;
        }
    }

    for (int n = 1; n <= x; ++n) {
        if (can_reach_losing(x - n, y - n, z)) {
            return true;
        }
    }
    for (int n = 1; n <= x; ++n) {
        if (can_reach_losing(x - n, y, z - n)) {
            return true;
        }
    }
    for (int n = 1; n <= y; ++n) {
        if (can_reach_losing(x, y - n, z - n)) {
            return true;
        }
    }

    for (int n = 1; n <= x; ++n) {
        if (can_reach_losing(x - n, y - n, z - n)) {
            return true;
        }
    }

    return false;
}

BruteResult solve_bruteforce_parallel(const int limit, unsigned thread_count) {
    if (limit < 0) {
        throw std::invalid_argument("Brute-force limit must be non-negative");
    }

    const int side = limit + 1;
    BruteResult result;
    result.limit = limit;
    result.side = side;
    result.is_losing.assign(checked_cube_size(side), 0U);

    const auto index3 = [side](const int a, const int b, const int c) -> std::size_t {
        return (static_cast<std::size_t>(a) * static_cast<std::size_t>(side) +
                static_cast<std::size_t>(b)) *
                   static_cast<std::size_t>(side) +
               static_cast<std::size_t>(c);
    };

    thread_count = std::max(1U, thread_count);

    for (int total = 0; total <= 3 * limit; ++total) {
        std::vector<std::array<int, 3>> layer_states;

        const int x_hi = std::min(limit, total);
        for (int x = 0; x <= x_hi; ++x) {
            const int rem = total - x;
            const int y_hi = std::min(limit, rem / 2);
            for (int y = x; y <= y_hi; ++y) {
                const int z = rem - y;
                if (z < y || z > limit) {
                    continue;
                }
                layer_states.push_back({x, y, z});
            }
        }

        if (layer_states.empty()) {
            continue;
        }

        const unsigned workers =
            std::min<unsigned>(thread_count, static_cast<unsigned>(layer_states.size()));

        struct ChunkResult {
            u64 partial_sum = 0ULL;
            std::vector<std::size_t> losing_indices;
        };

        std::vector<ChunkResult> chunks(workers);

        const std::size_t per_worker =
            (layer_states.size() + static_cast<std::size_t>(workers) - 1U) /
            static_cast<std::size_t>(workers);

        auto worker_fn = [&](const unsigned worker_id) {
            const std::size_t begin = static_cast<std::size_t>(worker_id) * per_worker;
            const std::size_t end =
                std::min(layer_states.size(), begin + per_worker);

            ChunkResult& chunk = chunks[worker_id];
            chunk.losing_indices.reserve(end > begin ? (end - begin) / 4 + 1 : 1);

            for (std::size_t i = begin; i < end; ++i) {
                const int x = layer_states[i][0];
                const int y = layer_states[i][1];
                const int z = layer_states[i][2];

                if (has_move_to_losing_state(result.is_losing, side, x, y, z)) {
                    continue;
                }

                chunk.partial_sum += static_cast<u64>(x + y + z);
                chunk.losing_indices.push_back(index3(x, y, z));
            }
        };

        std::vector<std::thread> pool;
        pool.reserve(workers > 0U ? workers - 1U : 0U);

        for (unsigned worker = 1U; worker < workers; ++worker) {
            pool.emplace_back(worker_fn, worker);
        }

        worker_fn(0U);

        for (auto& t : pool) {
            t.join();
        }

        for (const ChunkResult& chunk : chunks) {
            result.losing_sum += chunk.partial_sum;
            for (const std::size_t idx : chunk.losing_indices) {
                result.is_losing[idx] = 1U;
            }
        }
    }

    return result;
}

bool check_equal_u64(const u64 actual, const u64 expected, const std::string& label) {
    if (actual == expected) {
        return true;
    }

    std::cerr << "Checkpoint failed for " << label << ": expected " << expected
              << ", got " << actual << '\n';
    return false;
}

bool check_bool(const bool actual, const bool expected, const std::string& label) {
    if (actual == expected) {
        return true;
    }

    std::cerr << "Checkpoint failed for " << label << ": expected "
              << (expected ? "true" : "false") << ", got "
              << (actual ? "true" : "false") << '\n';
    return false;
}

bool run_checkpoints(const unsigned thread_count, const int brute_limit) {
    if (brute_limit < 13) {
        std::cerr << "Checkpoint brute-force limit must be at least 13.\n";
        return false;
    }

    const BruteResult brute = solve_bruteforce_parallel(brute_limit, thread_count);
    const u64 fast_at_brute_limit = solve_fast(brute_limit);
    if (!check_equal_u64(fast_at_brute_limit, brute.losing_sum,
                         "fast vs brute at n=" + std::to_string(brute_limit))) {
        return false;
    }

    if (!check_bool(brute.query_is_losing(0, 1, 2), true,
                    "losing example (0,1,2)")) {
        return false;
    }
    if (!check_bool(brute.query_is_losing(1, 3, 3), true,
                    "losing example (1,3,3)")) {
        return false;
    }
    if (!check_bool(brute.query_is_losing(0, 0, 13), false,
                    "winning example (0,0,13)")) {
        return false;
    }
    if (!check_bool(brute.query_is_losing(0, 11, 11), false,
                    "winning example (0,11,11)")) {
        return false;
    }
    if (!check_bool(brute.query_is_losing(5, 5, 5), false,
                    "winning example (5,5,5)")) {
        return false;
    }

    const u64 sample = solve_fast(100);
    if (!check_equal_u64(sample, 173895ULL, "sample sum at n=100")) {
        return false;
    }

    return true;
}

} // namespace

int main(int argc, char** argv) {
    try {
        Options options;
        if (!parse_arguments(argc, argv, options)) {
            return 1;
        }

        if (options.limit < 0) {
            std::cerr << "--limit must be non-negative.\n";
            return 1;
        }
        if (options.checkpoint_bruteforce_limit < 0) {
            std::cerr << "--check-limit must be non-negative.\n";
            return 1;
        }

        const unsigned thread_count = pick_thread_count(options.requested_threads);

        if (options.run_checkpoints &&
            !run_checkpoints(thread_count, options.checkpoint_bruteforce_limit)) {
            return 1;
        }

        const u64 answer = solve_fast(options.limit);
        std::cout << answer << '\n';
    } catch (const std::exception& ex) {
        std::cerr << "Error: " << ex.what() << '\n';
        return 1;
    }

    return 0;
}

Python

def solve():
    limit = 1000
    side = limit + 1
    one = bytearray(side * side)
    two = bytearray(side * side)
    all_tab = bytearray(side * side)

    total = 0
    for x in range(limit + 1):
        for y in range(x, limit + 1):
            for z in range(y, limit + 1):
                yz = y * side + z
                xz = x * side + z
                xy = x * side + y
                yx_z = (y - x) * side + z
                zy_x = (z - y) * side + x
                zx_y = (z - x) * side + y
                yx_zx = (y - x) * side + (z - x)
                if (one[yz] or one[xz] or one[xy] or
                    two[yx_z] or two[zy_x] or two[zx_y] or
                    all_tab[yx_zx]):
                    continue
                total += x + y + z
                one[yz] = 1
                one[xz] = 1
                one[xy] = 1
                two[yx_z] = 1
                two[zy_x] = 1
                two[zx_y] = 1
                all_tab[yx_zx] = 1
    return str(total)

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

Java

public class Euler260 {
    public static String solve() {
        int limit = 1000;
        int side = limit + 1;
        int tableSize = side * side;

        byte[] one = new byte[tableSize];
        byte[] two = new byte[tableSize];
        byte[] all = new byte[tableSize];

        long sum = 0;

        for (int x = 0; x <= limit; ++x) {
            for (int y = x; y <= limit; ++y) {
                int idXy = x * side + y;
                if (one[idXy] != 0)
                    continue;

                int dx = y - x;
                int ix = x * side;
                int iy = y * side;
                int idx = dx * side;

                for (int z = y; z <= limit; ++z) {
                    int idYz = iy + z;
                    int idXz = ix + z;
                    int idDxZ = idx + z;
                    int idDzX = (z - y) * side + x;
                    int idDzY = (z - x) * side + y;
                    int idDxDz = idx + (z - x);

                    if (one[idYz] != 0 || one[idXz] != 0 ||
                            two[idDxZ] != 0 || two[idDzX] != 0 || two[idDzY] != 0 ||
                            all[idDxDz] != 0) {
                        continue;
                    }

                    sum += (x + y + z);

                    one[idYz] = 1;
                    one[idXz] = 1;
                    one[idXy] = 1;

                    two[idDxZ] = 1;
                    two[idDzX] = 1;
                    two[idDzY] = 1;

                    all[idDxDz] = 1;
                    break;
                }
            }
        }

        return String.valueOf(sum);
    }

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