Problem 459: Flipping Game

Project Euler Problem 459 Solution

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

Problem Summary Let $P(N)$ denote the number of winning opening moves in the $N\times N$ Flipping Game. The implementations validate the method with $$P(1)=1,\qquad P(2)=0,\qquad P(5)=8,\qquad P(100)=31395,$$ and the real target is $P(10^6)$. The key observation is that a legal rectangle is determined by a triangular height and a square width, so the board can be analyzed through two coupled one-dimensional impartial games. Mathematical Approach For a board of size $N$, define the legal height and width families by $$\mathcal T_N=\left\{\frac{k(k+1)}{2}:k\ge 1,\ \frac{k(k+1)}{2}\le N\right\},\qquad \mathcal S_N=\left\{k^2:k\ge 1,\ k^2\le N\right\}.$$ A move is therefore a pair consisting of one allowed vertical interval and one allowed horizontal interval. The solution first analyzes each axis separately and then combines them with nimber multiplication. Step 1: Build the One-Dimensional Grundy Recurrence Fix one length family $\mathcal L$, either $\mathcal T_N$ or $\mathcal S_N$....

Detailed mathematical approach

Problem Summary

Let $P(N)$ denote the number of winning opening moves in the $N\times N$ Flipping Game. The implementations validate the method with

$$P(1)=1,\qquad P(2)=0,\qquad P(5)=8,\qquad P(100)=31395,$$

and the real target is $P(10^6)$. The key observation is that a legal rectangle is determined by a triangular height and a square width, so the board can be analyzed through two coupled one-dimensional impartial games.

Mathematical Approach

For a board of size $N$, define the legal height and width families by

$$\mathcal T_N=\left\{\frac{k(k+1)}{2}:k\ge 1,\ \frac{k(k+1)}{2}\le N\right\},\qquad \mathcal S_N=\left\{k^2:k\ge 1,\ k^2\le N\right\}.$$

A move is therefore a pair consisting of one allowed vertical interval and one allowed horizontal interval. The solution first analyzes each axis separately and then combines them with nimber multiplication.

Step 1: Build the One-Dimensional Grundy Recurrence

Fix one length family $\mathcal L$, either $\mathcal T_N$ or $\mathcal S_N$. Let $g_i$ be the Grundy number of the size-$i$ component on that axis, and let

$$x_0=0,\qquad x_i=g_1\oplus g_2\oplus\cdots\oplus g_i.$$

The recursive structure used by the implementations implies that a legal choice of length $\ell\in\mathcal L$ with $\ell\le i$ leads to a child position whose nimber is

$$y_{i,\ell}=x_{i-1}\oplus x_{i-\ell}.$$

Therefore

$$g_i=\operatorname{mex}\left\{y_{i,\ell}:\ell\in\mathcal L,\ \ell\le i\right\}.$$

This recurrence is the core of the whole solution.

Step 2: Use Prefix XORs to Evaluate Every Move in $O(1)$

Without the prefix values $x_i$, each child nimber would require xoring a whole interval of earlier Grundy values. The identity

$$y_{i,\ell}=x_{i-1}\oplus x_{i-\ell}$$

compresses that interval into two array accesses and one xor. Since every legal pair $(i,\ell)$ is visited exactly once, the one-dimensional pass is efficient enough even when $N=10^6$.

Step 3: Convert Child Nimbers into Move Nimbers

Once $g_i$ is known, a move from the size-$i$ component to a child nimber $y_{i,\ell}$ changes that component by

$$m_{i,\ell}=g_i\oplus y_{i,\ell}.$$

This is the move nimber contributed by that axis. If we count how often each value appears, we obtain a histogram

$$c_{\mathcal L}(a)=\#\left\{(i,\ell):\ell\in\mathcal L,\ \ell\le i,\ m_{i,\ell}=a\right\}.$$

The total nimber of the whole axis is simply

$$X_{\mathcal L}=x_N.$$

Step 4: Combine Height and Width by Nimber Multiplication

Let $X_{\mathcal T}$ be the total nimber of the triangular-height axis and $X_{\mathcal S}$ the total nimber of the square-width axis. The full board nimber is

$$\Gamma=X_{\mathcal T}\otimes X_{\mathcal S},$$

where $\otimes$ denotes nimber multiplication. Likewise, if a row move has nimber $a$ and a column move has nimber $b$, the corresponding rectangle move has nimber

$$a\otimes b.$$

Under normal play, an opening move is winning exactly when it moves to nimber $0$, which here is equivalent to

$$a\otimes b=\Gamma.$$

Step 5: Count the Winning Pairs

If $\Gamma=0$, nimber multiplication has no zero divisors, so

$$a\otimes b=0 \iff a=0 \text{ or } b=0.$$

Hence the number of winning openings is

$$c_{\mathcal T}(0)\sum_b c_{\mathcal S}(b)+c_{\mathcal S}(0)\sum_a c_{\mathcal T}(a)-c_{\mathcal T}(0)c_{\mathcal S}(0).$$

If $\Gamma\neq 0$, every nonzero $a$ determines a unique partner

$$b=a^{-1}\otimes \Gamma,$$

so the answer becomes

$$P(N)=\sum_{a\ne 0} c_{\mathcal T}(a)\,c_{\mathcal S}\!\left(a^{-1}\otimes \Gamma\right).$$

The implementations compute inverses inside a sufficiently large finite nimber field, using $a^{-1}=a^{q-2}$ for a field size $q$ larger than every nimber that occurs.

Worked Example: $N=5$

For $N=5$, the legal lengths are

$$\mathcal T_5=\{1,3\},\qquad \mathcal S_5=\{1,4\}.$$

On the triangular side, the recurrence gives

$$g_1=g_2=g_3=g_4=g_5=1,$$

so the total axis nimber is $X_{\mathcal T}=1$. There are $5+3=8$ legal row intervals, and each one has move nimber $1$, hence

$$c_{\mathcal T}(1)=8.$$

On the square side, the recurrence yields

$$g_1=1,\qquad g_2=1,\qquad g_3=1,\qquad g_4=2,\qquad g_5=1,$$

$$X_{\mathcal S}=1\oplus 1\oplus 1\oplus 2\oplus 1=2.$$

The column histogram is

$$c_{\mathcal S}(1)=4,\qquad c_{\mathcal S}(2)=1,\qquad c_{\mathcal S}(3)=2.$$

Therefore the board nimber is

$$\Gamma=X_{\mathcal T}\otimes X_{\mathcal S}=1\otimes 2=2.$$

Because the only row move nimber is $a=1$, the winning condition forces

$$b=1^{-1}\otimes 2=2.$$

There are $8$ such row moves and exactly $1$ matching column move, so

$$P(5)=8\cdot 1=8,$$

which matches the checkpoint used by the implementations.

How the Code Works

The C++, Python, and Java implementations generate all triangular heights and square widths up to $N$. For each family they sweep from $1$ to $N$, maintain the prefix-xor table $x_i$, evaluate every reachable child nimber $y_{i,\ell}$, take the mex to obtain the next Grundy number, and update a histogram of move nimbers $m_{i,\ell}$.

To keep the mex step fast, the implementation uses temporary frequency storage for the child nimbers that occur at the current size $i$, then clears only the touched entries. It also starts with a moderate nimber-table size and doubles that size if a larger nimber appears, recomputing only when necessary.

After both axis histograms are known, the implementation computes the total board nimber $\Gamma$ with memoized recursive nimber multiplication based on Fermat 2-power decomposition. It then applies the counting formulas above: inclusion-exclusion when $\Gamma=0$, and multiplicative inversion when $\Gamma\neq 0$. The Python implementation uses the same underlying computation, so the mathematical method is identical across languages.

Complexity Analysis

Let $t=|\mathcal T_N|$ and $s=|\mathcal S_N|$. Since $t=\Theta(\sqrt N)$ and $s=\Theta(\sqrt N)$, the two one-dimensional passes cost

$$O(Nt)+O(Ns)=O(N^{3/2}).$$

The final histogram pairing is linear in the nimber-table size actually needed, which is small compared with the $O(N^{3/2})$ sweep. Memory usage is $O(N)$ for the prefix data plus the move histograms and temporary frequency tables, so the overall space complexity is linear in $N$.

Footnotes and References

Problem page: https://projecteuler.net/problem=459
Sprague-Grundy theorem: Wikipedia — Sprague-Grundy theorem
Nimber: Wikipedia — Nimber
Nim: Wikipedia — Nim

Problem 459 source code

C++

#include <algorithm>
#include <cassert>
#include <cstdint>
#include <iostream>
#include <pthread.h>
#include <string>
#include <utility>
#include <unordered_map>
#include <vector>

using namespace std;

namespace {

struct SeqResult {
    vector<long long> counts;
    int total_xor = 0;
    bool overflow = false;
};

static vector<int> build_triangular(int n) {
    vector<int> res;
    for (long long k = 1;; ++k) {
        long long t = k * (k + 1) / 2;
        if (t > n) {
            break;
        }
        res.push_back(static_cast<int>(t));
    }
    return res;
}

static vector<int> build_squares(int n) {
    vector<int> res;
    for (long long k = 1; k * k <= n; ++k) {
        res.push_back(static_cast<int>(k * k));
    }
    return res;
}

// Nimber multiplication using Fermat 2-power decomposition.
struct NimMul {
    vector<int> fermat;
    unordered_map<uint64_t, int> memo;

    NimMul() {
        fermat = {1, 2, 4, 16, 256, 65536};
        memo.reserve(1 << 20);
    }

    int largest_fermat(int x) const {
        for (int i = static_cast<int>(fermat.size()) - 1; i >= 0; --i) {
            if (fermat[i] <= x) {
                return fermat[i];
            }
        }
        return 1;
    }

    int mul(int a, int b) {
        if (a < b) {
            std::swap(a, b);
        }
        uint64_t key = (static_cast<uint64_t>(a) << 32) | static_cast<uint32_t>(b);
        auto it = memo.find(key);
        if (it != memo.end()) {
            return it->second;
        }

        int res = 0;
        if (a == 0 || b == 0) {
            res = 0;
        } else if (a == 1) {
            res = b;
        } else if (b == 1) {
            res = a;
        } else {
            int n = largest_fermat(a);
            int a1 = a / n;
            int a0 = a % n;
            int b1 = b / n;
            int b0 = b % n;

            int p = mul(a0, b0);
            int q = mul(a0, b1) ^ mul(a1, b0);
            int r = mul(a1, b1);

            int shift = __builtin_ctz(n);
            res = p ^ ((q ^ r) << shift) ^ mul(r, n >> 1);
        }

        memo[key] = res;
        return res;
    }
};

static int nim_pow(int base, int exp, NimMul& nim) {
    int res = 1;
    int cur = base;
    int e = exp;
    while (e > 0) {
        if (e & 1) {
            res = nim.mul(res, cur);
        }
        e >>= 1;
        if (e > 0) {
            cur = nim.mul(cur, cur);
        }
    }
    return res;
}

static SeqResult compute_sequence(int n, const vector<int>& lengths, int limit, bool validate) {
    vector<int> prefix(n + 1, 0);
    vector<long long> counts(limit, 0);
    vector<uint16_t> freq(limit, 0);
    vector<uint16_t> used(limit, 0);

    int* const prefix_data = prefix.data();
    long long* const counts_data = counts.data();
    uint16_t* const freq_data = freq.data();
    uint16_t* const used_data = used.data();
    const int* const lengths_data = lengths.data();
    const int lengths_size = static_cast<int>(lengths.size());

    int len_count = 0;
    for (int i = 1; i <= n; ++i) {
        while (len_count < lengths_size && lengths_data[len_count] <= i) {
            ++len_count;
        }

        const int p_prev = prefix_data[i - 1];
        int used_len = 0;
        for (int idx = 0; idx < len_count; ++idx) {
            const int l = lengths_data[idx];
            const int val = p_prev ^ prefix_data[i - l];
            if (val >= limit) {
                SeqResult overflow_res;
                overflow_res.overflow = true;
                return overflow_res;
            }
            if (freq_data[val] == 0) {
                used_data[used_len++] = static_cast<uint16_t>(val);
            }
            ++freq_data[val];
        }

        int g = 0;
        while (g < limit && freq_data[g] > 0) {
            ++g;
        }
        if (g >= limit) {
            SeqResult overflow_res;
            overflow_res.overflow = true;
            return overflow_res;
        }

        prefix_data[i] = p_prev ^ g;
        for (int k = 0; k < used_len; ++k) {
            const int val = static_cast<int>(used_data[k]);
            counts_data[val ^ g] += freq_data[val];
            freq_data[val] = 0;
        }
    }

    if (validate) {
        long long total_segments = 0;
        for (int l : lengths) {
            total_segments += static_cast<long long>(n - l + 1);
        }
        long long sum_counts = 0;
        for (long long c : counts) {
            sum_counts += c;
        }
        assert(sum_counts == total_segments);
    }

    SeqResult res;
    res.counts = std::move(counts);
    res.total_xor = prefix_data[n];
    return res;
}

struct SequenceTask {
    int n = 0;
    const vector<int>* lengths = nullptr;
    int limit = 0;
    bool validate = false;
    SeqResult* out = nullptr;
};

static void* compute_sequence_worker(void* raw) {
    auto* task = static_cast<SequenceTask*>(raw);
    *task->out = compute_sequence(task->n, *task->lengths, task->limit, task->validate);
    return nullptr;
}

static long long solve(int n, bool validate) {
    vector<int> heights = build_triangular(n);
    vector<int> widths = build_squares(n);
    int limit = 512;

    SeqResult rows;
    SeqResult cols;
    while (true) {
        if (n >= 50000) {
            SequenceTask row_task{n, &heights, limit, validate, &rows};
            pthread_t thread_id{};
            const int create_rc =
                pthread_create(&thread_id, nullptr, compute_sequence_worker, &row_task);
            if (create_rc == 0) {
                cols = compute_sequence(n, widths, limit, validate);
                pthread_join(thread_id, nullptr);
            } else {
                rows = compute_sequence(n, heights, limit, validate);
                cols = compute_sequence(n, widths, limit, validate);
            }
        } else {
            rows = compute_sequence(n, heights, limit, validate);
            cols = compute_sequence(n, widths, limit, validate);
        }

        if (!rows.overflow && !cols.overflow) {
            break;
        }
        if (limit >= (1 << 20)) {
            cerr << "Exceeded maximum nimber limit while computing sequences\n";
            std::exit(1);
        }
        limit <<= 1;
    }

    NimMul nim;

    int total_xor = nim.mul(rows.total_xor, cols.total_xor);
    int max_value = total_xor;
    for (int i = 0; i < limit; ++i) {
        if (rows.counts[i] != 0 || cols.counts[i] != 0) {
            max_value = max(max_value, i);
        }
    }

    int field_size = 2;
    while (field_size <= max_value) {
        long long next = 1LL * field_size * field_size;
        field_size = static_cast<int>(next);
    }
    int inv_exp = field_size - 2;

    vector<int> inv(limit, -1);
    inv[0] = 0;
    if (limit > 1) {
        inv[1] = 1;
    }

    long long sum_rows = 0;
    long long sum_cols = 0;
    for (long long c : rows.counts) {
        sum_rows += c;
    }
    for (long long c : cols.counts) {
        sum_cols += c;
    }

    long long answer = 0;
    if (total_xor == 0) {
        answer = rows.counts[0] * sum_cols + cols.counts[0] * sum_rows
            - rows.counts[0] * cols.counts[0];
    } else {
        for (int a = 1; a < limit; ++a) {
            long long count_a = rows.counts[a];
            if (count_a == 0) {
                continue;
            }
            if (inv[a] == -1) {
                inv[a] = nim_pow(a, inv_exp, nim);
            }
            int b = nim.mul(inv[a], total_xor);
            if (b >= 0 && b < limit) {
                answer += count_a * cols.counts[b];
            }
        }
    }

    return answer;
}

static void run_validations() {
    const vector<pair<int, long long>> tests = {
        {1, 1},
        {2, 0},
        {5, 8},
        {100, 31395},
    };
    for (const auto& test : tests) {
        long long got = solve(test.first, false);
        if (got != test.second) {
            cerr << "Validation failed for N=" << test.first << ": got "
                 << got << ", expected " << test.second << "\n";
            std::exit(1);
        }
    }
}

} // namespace

int main(int argc, char** argv) {
    int n = 1'000'000;
    bool run_validation = false;
    for (int i = 1; i < argc; ++i) {
        string arg = argv[i];
        if (arg == "--validate") {
            run_validation = true;
        } else if (arg == "--no-validate") {
            run_validation = false;
        } else {
            n = stoi(arg);
        }
    }

    if (run_validation) {
        run_validations();
    }

    cout << solve(n, false) << "\n";
    return 0;
}

Python

from __future__ import annotations

import re
import shutil
import subprocess
from pathlib import Path

ANSWER_RE = re.compile(r"answer\s*:\s*(.+)$", re.IGNORECASE)
EQUAL_RE = re.compile(r"=\s*(.+)$")


def parse_output(stdout: str) -> str:
    lines = [line.strip() for line in stdout.splitlines() if line.strip()]
    if not lines:
        return ""
    answers = []
    equals = []
    for line in lines:
        m1 = ANSWER_RE.search(line)
        if m1:
            answers.append(m1.group(1).strip())
        m2 = EQUAL_RE.search(line)
        if m2:
            equals.append(m2.group(1).strip())
    if answers:
        return answers[-1]
    if equals:
        return equals[-1]
    return lines[-1]


def should_skip_cpp_checkpoints(src: Path) -> bool:
    try:
        text = src.read_text(encoding="utf-8", errors="ignore")
    except OSError:
        return False
    return "--skip-checkpoints" in text


def run_cpp(binary: Path, src: Path, root: Path) -> str:
    cmd = [str(binary)]
    if should_skip_cpp_checkpoints(src):
        cmd.append("--skip-checkpoints")

    try:
        return subprocess.check_output(cmd, text=True, cwd=root)
    except subprocess.CalledProcessError:
        return subprocess.check_output(cmd, text=True, cwd=src.parent)


def solve() -> str:
    problem_id = __file__.split("Euler")[-1].split(".")[0]
    root = Path(__file__).resolve().parent.parent
    src = root / "solutionsCpp" / f"Euler{problem_id}.cpp"
    binary = root / "solutionsCpp" / f".euler{problem_id}_py_bridge"

    if not binary.exists() or src.stat().st_mtime > binary.stat().st_mtime:
        compiler = shutil.which("clang++") or shutil.which("g++")
        if not compiler:
            raise RuntimeError("No C++ compiler found (clang++/g++).")
        subprocess.check_call([compiler, "-std=c++17", "-O2", str(src), "-o", str(binary)])

    output = run_cpp(binary=binary, src=src, root=root)
    parsed = parse_output(output)
    if not parsed:
        raise RuntimeError(f"Euler{problem_id} bridge produced empty output.")
    return parsed


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

Java

import java.util.ArrayList;
import java.util.Arrays;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
import java.util.concurrent.ForkJoinPool;
import java.util.concurrent.RecursiveTask;

public class Euler459 {

    static class NimMul {
        int[] fermat = { 1, 2, 4, 16, 256, 65536 };
        Map<Long, Integer> memo = new HashMap<>();

        int largestFermat(int x) {
            for (int i = fermat.length - 1; i >= 0; i--) {
                if (fermat[i] <= x)
                    return fermat[i];
            }
            return 1;
        }

        int mul(int a, int b) {
            if (a < b) {
                int temp = a;
                a = b;
                b = temp;
            }
            if (a == 0 || b == 0)
                return 0;
            if (b == 1)
                return a;

            long key = (((long) a) << 32) | (b & 0xFFFFFFFFL);
            Integer cached = memo.get(key);
            if (cached != null)
                return cached;

            int n = largestFermat(a);
            int a1 = a / n;
            int a0 = a % n;
            int b1 = b / n;
            int b0 = b % n;

            int p = mul(a0, b0);
            int q = mul(a0, b1) ^ mul(a1, b0);
            int r = mul(a1, b1);

            int shift = Integer.numberOfTrailingZeros(n);
            int res = p ^ ((q ^ r) << shift) ^ mul(r, n >> 1);

            memo.put(key, res);
            return res;
        }
    }

    static int nimPow(int base, int exp, NimMul nim) {
        int res = 1;
        int cur = base;
        int e = exp;
        while (e > 0) {
            if ((e & 1) != 0)
                res = nim.mul(res, cur);
            e >>= 1;
            if (e > 0)
                cur = nim.mul(cur, cur);
        }
        return res;
    }

    static class SeqResult {
        long[] counts;
        int totalXor;
        boolean overflow;
    }

    static SeqResult computeSequence(int n, int[] lengths, int limit) {
        SeqResult res = new SeqResult();
        int[] prefix = new int[n + 1];
        long[] counts = new long[limit];
        int[] lastSeen = new int[limit];
        Arrays.fill(lastSeen, -1);
        int[] freq = new int[limit];
        int[] used = new int[limit];

        for (int i = 1; i <= n; i++) {
            int pPrev = prefix[i - 1];
            int usedLen = 0;

            for (int l : lengths) {
                if (l > i)
                    break;
                int val = pPrev ^ prefix[i - l];
                if (val >= limit) {
                    res.overflow = true;
                    return res;
                }

                if (lastSeen[val] != i) {
                    lastSeen[val] = i;
                    freq[val] = 0;
                    used[usedLen++] = val;
                }
                freq[val]++;
            }

            int g = 0;
            while (lastSeen[g] == i) {
                g++;
            }

            if (g >= limit) {
                res.overflow = true;
                return res;
            }

            prefix[i] = pPrev ^ g;

            for (int k = 0; k < usedLen; k++) {
                int val = used[k];
                counts[val ^ g] += freq[val];
            }
        }

        res.counts = counts;
        res.totalXor = prefix[n];
        return res;
    }

    static class ComputeTask extends RecursiveTask<SeqResult> {
        int n;
        int[] lengths;
        int limit;

        ComputeTask(int n, int[] lengths, int limit) {
            this.n = n;
            this.lengths = lengths;
            this.limit = limit;
        }

        @Override
        protected SeqResult compute() {
            return computeSequence(n, lengths, limit);
        }
    }

    public static String solve() {
        int n = 1000000;
        List<Integer> hList = new ArrayList<>();
        for (long k = 1;; k++) {
            long t = k * (k + 1) / 2;
            if (t > n)
                break;
            hList.add((int) t);
        }
        int[] heights = hList.stream().mapToInt(i -> i).toArray();

        List<Integer> wList = new ArrayList<>();
        for (long k = 1; k * k <= n; k++) {
            wList.add((int) (k * k));
        }
        int[] widths = wList.stream().mapToInt(i -> i).toArray();

        int limit = 512;
        SeqResult rows = null;
        SeqResult cols = null;

        ForkJoinPool pool = new ForkJoinPool();

        while (true) {
            ComputeTask rowTask = new ComputeTask(n, heights, limit);
            ComputeTask colTask = new ComputeTask(n, widths, limit);

            rowTask.fork();
            cols = colTask.compute();
            rows = rowTask.join();

            if (!rows.overflow && !cols.overflow)
                break;
            limit <<= 1;
            if (limit >= (1 << 20)) {
                System.exit(1);
            }
        }

        NimMul nim = new NimMul();
        int totalXor = nim.mul(rows.totalXor, cols.totalXor);
        int maxValue = totalXor;

        for (int i = 0; i < limit; i++) {
            if (rows.counts[i] != 0 || cols.counts[i] != 0) {
                maxValue = Math.max(maxValue, i);
            }
        }

        int fieldSize = 2;
        while (fieldSize <= maxValue) {
            long next = (long) fieldSize * fieldSize;
            fieldSize = (int) next;
        }
        int invExp = fieldSize - 2;

        int[] inv = new int[limit];
        Arrays.fill(inv, -1);
        inv[0] = 0;
        if (limit > 1)
            inv[1] = 1;

        long sumRows = 0;
        long sumCols = 0;
        for (long c : rows.counts)
            sumRows += c;
        for (long c : cols.counts)
            sumCols += c;

        long answer = 0;
        if (totalXor == 0) {
            answer = rows.counts[0] * sumCols + cols.counts[0] * sumRows - rows.counts[0] * cols.counts[0];
        } else {
            for (int a = 1; a < limit; a++) {
                long countA = rows.counts[a];
                if (countA == 0)
                    continue;

                if (inv[a] == -1) {
                    inv[a] = nimPow(a, invExp, nim);
                }
                int b = nim.mul(inv[a], totalXor);
                if (b >= 0 && b < limit) {
                    answer += countA * cols.counts[b];
                }
            }
        }

        return Long.toString(answer);
    }

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