Problem 888: 1249 Nim

Project Euler Problem 888 Solution

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

Problem Summary We study an impartial heap game in which a heap of size $n$ may be changed in two ways: remove $1$, $2$, $4$, or $9$ counters, or split the heap into two positive heaps whose sizes add up to $n$. For the large parameters $N=12{,}491{,}249$ and $m=1249$, the required quantity $S(N,m)$ is the number of size-$m$ multisets of heap sizes chosen from $\{1,2,\dots,N\}$ whose combined position is losing, taken modulo $912491249$. By the Sprague-Grundy theorem, a disjoint sum of heaps is losing exactly when the xor of the single-heap Grundy values is $0$. The whole problem therefore becomes: compute the Grundy value $g(n)$ for one heap, understand how often each Grundy value appears among $1\le n\le N$, and count size-$m$ multisets whose xor is zero. Mathematical Approach The implementations combine game theory, eventual periodicity, and a compact xor dynamic program. The key point is that the huge bound $N$ is handled through periodic structure, not by computing every $g(n)$ all the way to $N$. Step 1: Write the Single-Heap Grundy Recurrence Let $g(n)$ be the Grundy value of one heap of size $n$. Then $g(0)=0$, and for $n\ge 1$ every legal move contributes one reachable Grundy value. The subtraction moves contribute $g(n-1)$, $g(n-2)$, $g(n-4)$, and $g(n-9)$ whenever the index is nonnegative....

Detailed mathematical approach

Problem Summary

We study an impartial heap game in which a heap of size $n$ may be changed in two ways: remove $1$, $2$, $4$, or $9$ counters, or split the heap into two positive heaps whose sizes add up to $n$. For the large parameters $N=12{,}491{,}249$ and $m=1249$, the required quantity $S(N,m)$ is the number of size-$m$ multisets of heap sizes chosen from $\{1,2,\dots,N\}$ whose combined position is losing, taken modulo $912491249$.

By the Sprague-Grundy theorem, a disjoint sum of heaps is losing exactly when the xor of the single-heap Grundy values is $0$. The whole problem therefore becomes: compute the Grundy value $g(n)$ for one heap, understand how often each Grundy value appears among $1\le n\le N$, and count size-$m$ multisets whose xor is zero.

Mathematical Approach

The implementations combine game theory, eventual periodicity, and a compact xor dynamic program. The key point is that the huge bound $N$ is handled through periodic structure, not by computing every $g(n)$ all the way to $N$.

Step 1: Write the Single-Heap Grundy Recurrence

Let $g(n)$ be the Grundy value of one heap of size $n$. Then $g(0)=0$, and for $n\ge 1$ every legal move contributes one reachable Grundy value. The subtraction moves contribute $g(n-1)$, $g(n-2)$, $g(n-4)$, and $g(n-9)$ whenever the index is nonnegative. A split $n=a+(n-a)$ contributes the xor $g(a)\oplus g(n-a)$.

Therefore

$$g(n)=\operatorname{mex}\left(\{g(n-s): s\in\{1,2,4,9\},\ s\le n\}\cup\{g(a)\oplus g(n-a):1\le a\le \lfloor n/2\rfloor\}\right).$$

The upper limit $\lfloor n/2\rfloor$ is enough because the split $a,n-a$ is the same position as $n-a,a$. Computing the first values gives

$$g(1..10)=1,2,0,3,4,6,1,2,5,3.$$

These initial values already show that the sequence is nontrivial: subtraction moves alone would be much simpler, but the split move injects xor structure into the recurrence.

Step 2: Separate Even and Odd Indices and Exploit Periodicity

The implementations do not try to prove a closed form for $g(n)$. Instead they compute a long prefix and then look for eventual periodicity separately in the even and odd subsequences

$$e_k=g(2k),\qquad o_k=g(2k-1).$$

This parity split is exactly what works in practice. The computed data show that the even subsequence becomes periodic with period $79$ after the first $22$ even terms, while the odd subsequence becomes periodic with period $70$ after the first $161$ odd terms.

If $c_r(N)$ denotes the number of heap sizes $n\in[1,N]$ with Grundy value $r$, then

$$c_r(N)=c_r^{\text{odd}}\left(\left\lceil\frac{N}{2}\right\rceil\right)+c_r^{\text{even}}\left(\left\lfloor\frac{N}{2}\right\rfloor\right).$$

For each parity subsequence, once a preperiod and a period are known, the count of any value $r$ among the first $T$ terms is obtained by prefix plus whole cycles plus a short tail. That turns the enormous value $N=12{,}491{,}249$ into a straightforward counting exercise.

Step 3: Collapse Heap Sizes into Grundy Classes

After periodic counting, all heap sizes are grouped only by their Grundy value. Let

$$c_r=\#\{n\in\{1,\dots,N\}:g(n)=r\}.$$

Now consider one fixed class $r$. There are $c_r$ distinct heap sizes inside this class, and the game position allows repeated heap sizes, so choosing exactly $k$ heaps from that class is a multiset-counting problem. The number of possibilities is

$$\binom{c_r+k-1}{k}.$$

This is the usual combinations-with-repetition formula, equivalently the coefficient of $z^k$ in

$$\frac{1}{(1-z)^{c_r}}=\sum_{k\ge 0}\binom{c_r+k-1}{k}z^k.$$

The xor contribution of those $k$ heaps is especially simple: since $r\oplus r=0$, the total xor from class $r$ is $0$ when $k$ is even and $r$ when $k$ is odd.

Step 4: Run a Dynamic Program on Xor and Heap Count

Process the Grundy classes one by one. Let $DP[\xi][t]$ be the number of ways to choose $t$ heaps from the classes processed so far so that the current xor is $\xi$. Initially, before any class is used,

$$DP[0][0]=1.$$

When the class $r$ is added, choosing $k$ heaps from it changes the count by $k$ and changes the xor only according to the parity of $k$. Thus the transition is

$$DP'[\xi][t+k]\mathrel{+}=DP[\xi][t]\binom{c_r+k-1}{k}\qquad\text{for even }k,$$

$$DP'[\xi\oplus r][t+k]\mathrel{+}=DP[\xi][t]\binom{c_r+k-1}{k}\qquad\text{for odd }k.$$

After every Grundy class has been processed, the answer is

$$S(N,m)=DP[0][m]\pmod{912491249}.$$

The computed prefix shows that only Grundy values $0$ through $10$ occur, so the xor state space has size $16$, which keeps the dynamic program compact.

Step 5: Worked Example $(N,m)=(10,2)$

Using the first ten Grundy values

$$g(1..10)=1,2,0,3,4,6,1,2,5,3,$$

the class counts are

$$c_0=1,\quad c_1=2,\quad c_2=2,\quad c_3=2,\quad c_4=1,\quad c_5=1,\quad c_6=1.$$

With exactly two heaps, the total xor is $0$ precisely when both chosen heaps come from the same Grundy class. So

$$S(10,2)=\sum_r \binom{c_r+1}{2}=\binom{2}{2}+\binom{3}{2}+\binom{3}{2}+\binom{3}{2}+\binom{2}{2}+\binom{2}{2}+\binom{2}{2}=13.$$

This matches the dynamic-program interpretation perfectly: every class contributes its even-sized selections to xor $0$, and for $m=2$ that means selecting two heaps from one class.

How the Code Works

The implementation first computes a prefix of the Grundy sequence up to $20000$. For each heap size it gathers all reachable Grundy values into a bit mask and then takes the mex. That gives both the sequence itself and the maximum Grundy value appearing in the prefix.

Next it extracts the even-indexed and odd-indexed subsequences and searches for a short period in each one. Candidate periods up to $300$ are tested, the tail is required to contain at least six consecutive repetitions, and then the start is moved left to the earliest index where the same period still works. In the final data this yields period $79$ starting after $22$ even terms and period $70$ starting after $161$ odd terms.

Using those periodic descriptions, the implementation counts how many times each Grundy value occurs in $\{1,\dots,N\}$. For each class size $c_r$, it then builds the coefficients $\binom{c_r+k-1}{k}\bmod 912491249$ for $0\le k\le m$. The coefficients are split into even-$k$ and odd-$k$ lists so that the xor update is just “stay in the same mask” or “toggle by $r$”.

Finally, a rolling two-layer dynamic program over xor masks and selected-heap count is executed. Since the observed Grundy values lie in $0,\dots,10$, only $16$ xor masks are needed, and the answer is the entry with xor $0$ and total size $1249$.

Complexity Analysis

Let $M=20000$ be the precomputation cutoff, let $R$ be the number of Grundy classes actually used, and let $X$ be the number of xor masks. The Grundy prefix computation costs

$$O\left(\sum_{n=1}^{M} n\right)=O(M^2)$$

time, because every $n$ checks all splits up to $n/2$, and it uses $O(M)$ memory for the stored sequence. The period search over the even and odd subsequences is much smaller and is negligible compared with the quadratic prefix build.

The counting stage after periodic compression is linear in the period lengths. The main combinatorial stage uses a dynamic program of size $X\times (m+1)$, updated for each class over all $k\le m$, so its running time is $O(RXm^2)$ and its memory use is $O(Xm)$. In this problem $R=11$ and $X=16$, so the practical bottleneck is the one-time Grundy precomputation, not the final xor DP.

Footnotes and References

Problem page: Project Euler 888
Sprague-Grundy theorem: Wikipedia - Sprague-Grundy theorem
Mex function in combinatorial game theory: Wikipedia - mex
Nim and xor evaluation: Wikipedia - Nim
Combinations with repetition: Wikipedia - combinations with repetition

Problem 888 source code

C++

#include <algorithm>
#include <cstdint>
#include <iostream>
#include <thread>
#include <vector>
#include <boost/multiprecision/cpp_int.hpp>

using namespace std;
using boost::multiprecision::cpp_int;

namespace {

constexpr int64_t MOD = 912491249LL;

inline int64_t mod_add(int64_t a, int64_t b) {
    int64_t v = a + b;
    if (v >= MOD) v -= MOD;
    return v;
}

inline int64_t mod_mul(int64_t a, int64_t b) {
    return static_cast<int64_t>((__int128)a * b % MOD);
}

struct PeriodInfo {
    int period = 0;
    int start = 0;
    vector<int> seq;
};

struct GrundyData {
    vector<int> g;
    int max_g = 0;
    PeriodInfo even;
    PeriodInfo odd;
};

vector<int> compute_grundy(int max_n, int& max_g) {
    vector<int> g(static_cast<size_t>(max_n + 1), 0);
    max_g = 0;
    const int moves[] = {1, 2, 4, 9};

    for (int n = 1; n <= max_n; ++n) {
        uint64_t mask = 0;
        for (int mv : moves) {
            if (n >= mv) {
                int val = g[n - mv];
                if (val >= 63) {
                    cerr << "Grundy value too large for mask\n";
                    exit(1);
                }
                mask |= (1ULL << val);
            }
        }
        int half = n / 2;
        for (int a = 1; a <= half; ++a) {
            int val = g[a] ^ g[n - a];
            if (val >= 63) {
                cerr << "Grundy value too large for mask\n";
                exit(1);
            }
            mask |= (1ULL << val);
        }
        int mex = 0;
        while (mask & (1ULL << mex)) ++mex;
        if (mex >= 63) {
            cerr << "Grundy value too large for mask\n";
            exit(1);
        }
        g[n] = mex;
        if (mex > max_g) max_g = mex;
    }
    return g;
}

PeriodInfo detect_period(const vector<int>& seq, int max_period = 300, int min_reps = 6) {
    const int n = static_cast<int>(seq.size());
    for (int p = 1; p <= max_period; ++p) {
        int L = p * min_reps;
        if (L + p > n) continue;
        int start_tail = n - L;
        bool ok = true;
        for (int i = start_tail; i + p < n; ++i) {
            if (seq[i] != seq[i + p]) {
                ok = false;
                break;
            }
        }
        if (!ok) continue;
        for (int s = 0; s + p < n; ++s) {
            bool good = true;
            for (int i = s; i + p < n; ++i) {
                if (seq[i] != seq[i + p]) {
                    good = false;
                    break;
                }
            }
            if (good) {
                PeriodInfo info;
                info.period = p;
                info.start = s;
                info.seq = seq;
                return info;
            }
        }
    }
    return PeriodInfo{};
}

GrundyData build_grundy_data(int max_n) {
    GrundyData data;
    data.g = compute_grundy(max_n, data.max_g);

    vector<int> even;
    vector<int> odd;
    even.reserve(static_cast<size_t>(max_n / 2));
    odd.reserve(static_cast<size_t>((max_n + 1) / 2));

    for (int k = 1; 2 * k <= max_n; ++k) {
        even.push_back(data.g[2 * k]);
    }
    for (int k = 1; 2 * k - 1 <= max_n; ++k) {
        odd.push_back(data.g[2 * k - 1]);
    }

    data.even = detect_period(even);
    data.odd = detect_period(odd);
    if (data.even.period == 0 || data.odd.period == 0) {
        cerr << "Failed to detect periods\n";
        exit(1);
    }
    return data;
}

void add_counts(const PeriodInfo& info, long long total, vector<long long>& counts) {
    if (total <= 0) return;
    const vector<int>& seq = info.seq;
    const int start = info.start;
    const int period = info.period;

    if (total <= start) {
        for (long long i = 0; i < total; ++i) {
            counts[seq[static_cast<size_t>(i)]]++;
        }
        return;
    }

    for (int i = 0; i < start; ++i) {
        counts[seq[static_cast<size_t>(i)]]++;
    }

    long long rem = total - start;
    vector<long long> freq(counts.size(), 0);
    for (int i = 0; i < period; ++i) {
        freq[seq[static_cast<size_t>(start + i)]]++;
    }
    long long cycles = rem / period;
    long long tail = rem % period;
    for (size_t g = 0; g < counts.size(); ++g) {
        counts[g] += freq[g] * cycles;
    }
    for (int i = 0; i < tail; ++i) {
        counts[seq[static_cast<size_t>(start + i)]]++;
    }
}

vector<long long> grundy_counts(long long N, const GrundyData& data) {
    vector<long long> counts(static_cast<size_t>(data.max_g + 1), 0);
    long long total_even = N / 2;
    long long total_odd = (N + 1) / 2;
    add_counts(data.even, total_even, counts);
    add_counts(data.odd, total_odd, counts);
    return counts;
}

vector<int64_t> build_comb(long long c, int m) {
    vector<int64_t> comb(static_cast<size_t>(m + 1), 0);
    cpp_int val = 1;
    comb[0] = 1;
    for (int k = 1; k <= m; ++k) {
        val *= static_cast<long long>(c + k - 1);
        val /= k;
        comb[static_cast<size_t>(k)] = (val % MOD).convert_to<long long>();
    }
    return comb;
}

vector<vector<int64_t>> update_dp(const vector<vector<int64_t>>& dp,
                                  int gval,
                                  const vector<pair<int, int64_t>>& even_terms,
                                  const vector<pair<int, int64_t>>& odd_terms,
                                  int m) {
    const int states = static_cast<int>(dp.size());
    int thread_count = static_cast<int>(thread::hardware_concurrency());
    if (thread_count <= 0) thread_count = 1;
    if (thread_count > states) thread_count = states;

    vector<vector<int64_t>> next(states, vector<int64_t>(static_cast<size_t>(m + 1), 0));
    if (thread_count == 1) {
        for (int mask = 0; mask < states; ++mask) {
            const auto& cur = dp[mask];
            auto& dest_same = next[mask];
            auto& dest_xor = next[mask ^ gval];
            for (int i = 0; i <= m; ++i) {
                int64_t curv = cur[static_cast<size_t>(i)];
                if (curv == 0) continue;
                for (const auto& term : even_terms) {
                    int j = term.first;
                    if (i + j > m) break;
                    int64_t add = mod_mul(curv, term.second);
                    dest_same[static_cast<size_t>(i + j)] =
                        mod_add(dest_same[static_cast<size_t>(i + j)], add);
                }
                for (const auto& term : odd_terms) {
                    int j = term.first;
                    if (i + j > m) break;
                    int64_t add = mod_mul(curv, term.second);
                    dest_xor[static_cast<size_t>(i + j)] =
                        mod_add(dest_xor[static_cast<size_t>(i + j)], add);
                }
            }
        }
        return next;
    }

    vector<vector<vector<int64_t>>> partial(
        static_cast<size_t>(thread_count),
        vector<vector<int64_t>>(states, vector<int64_t>(static_cast<size_t>(m + 1), 0)));

    vector<thread> workers;
    workers.reserve(static_cast<size_t>(thread_count));
    for (int t = 0; t < thread_count; ++t) {
        workers.emplace_back([&, t]() {
            for (int mask = t; mask < states; mask += thread_count) {
                const auto& cur = dp[mask];
                auto& dest_same = partial[static_cast<size_t>(t)][mask];
                auto& dest_xor = partial[static_cast<size_t>(t)][mask ^ gval];
                for (int i = 0; i <= m; ++i) {
                    int64_t curv = cur[static_cast<size_t>(i)];
                    if (curv == 0) continue;
                    for (const auto& term : even_terms) {
                        int j = term.first;
                        if (i + j > m) break;
                        int64_t add = mod_mul(curv, term.second);
                        dest_same[static_cast<size_t>(i + j)] =
                            mod_add(dest_same[static_cast<size_t>(i + j)], add);
                    }
                    for (const auto& term : odd_terms) {
                        int j = term.first;
                        if (i + j > m) break;
                        int64_t add = mod_mul(curv, term.second);
                        dest_xor[static_cast<size_t>(i + j)] =
                            mod_add(dest_xor[static_cast<size_t>(i + j)], add);
                    }
                }
            }
        });
    }
    for (auto& th : workers) th.join();

    for (int t = 0; t < thread_count; ++t) {
        for (int mask = 0; mask < states; ++mask) {
            auto& dest = next[mask];
            const auto& src = partial[static_cast<size_t>(t)][mask];
            for (int i = 0; i <= m; ++i) {
                int64_t val = src[static_cast<size_t>(i)];
                if (val == 0) continue;
                dest[static_cast<size_t>(i)] = mod_add(dest[static_cast<size_t>(i)], val);
            }
        }
    }

    return next;
}

int64_t compute_S(long long N, int m, const GrundyData& data) {
    vector<long long> counts = grundy_counts(N, data);

    int max_g = data.max_g;
    int bits = 0;
    while ((1 << bits) <= max_g) ++bits;
    int states = 1 << bits;

    vector<vector<int64_t>> dp(states, vector<int64_t>(static_cast<size_t>(m + 1), 0));
    dp[0][0] = 1;

    for (int gval = 0; gval <= max_g; ++gval) {
        long long c = counts[static_cast<size_t>(gval)];
        if (c == 0) continue;

        vector<int64_t> comb = build_comb(c, m);
        vector<pair<int, int64_t>> even_terms;
        vector<pair<int, int64_t>> odd_terms;
        even_terms.reserve(static_cast<size_t>(m / 2 + 1));
        odd_terms.reserve(static_cast<size_t>(m / 2 + 1));
        for (int k = 0; k <= m; ++k) {
            int64_t v = comb[static_cast<size_t>(k)];
            if (v == 0) continue;
            if (k & 1) {
                odd_terms.push_back({k, v});
            } else {
                even_terms.push_back({k, v});
            }
        }

        dp = update_dp(dp, gval, even_terms, odd_terms, m);
    }

    return dp[0][static_cast<size_t>(m)];
}

}  // namespace

int main() {
    const int max_n = 20000;
    GrundyData data = build_grundy_data(max_n);

    const int64_t check1 = compute_S(12, 4, data);
    if (check1 != 204) {
        cerr << "Validation failure: S(12,4)\n";
        return 1;
    }

    const int64_t check2 = compute_S(124, 9, data);
    const int64_t expected2 = 2259208528408LL % MOD;
    if (check2 != expected2) {
        cerr << "Validation failure: S(124,9)\n";
        return 1;
    }

    const long long N = 12'491'249LL;
    const int m = 1249;
    cout << compute_S(N, m, data) << '\n';
    return 0;
}

Python

import sys

MOD = 912491249

def mod_add(a, b):
    return (a + b) % MOD

def compute_grundy(max_n):
    g = [0] * (max_n + 1)
    max_g = 0
    moves = [1, 2, 4, 9]

    for n in range(1, max_n + 1):
        mask = 0
        for mv in moves:
            if n >= mv:
                mask |= (1 << g[n - mv])
        
        half = n // 2
        for a in range(1, half + 1):
            mask |= (1 << (g[a] ^ g[n - a]))
            
        mex = 0
        while mask & (1 << mex):
            mex += 1
            
        g[n] = mex
        if mex > max_g:
            max_g = mex
            
    return g, max_g

def detect_period(seq, max_period=300, min_reps=6):
    n = len(seq)
    for p in range(1, max_period + 1):
        L = p * min_reps
        if L + p > n:
            continue
        start_tail = n - L
        
        ok = True
        for i in range(start_tail, n - p):
            if seq[i] != seq[i + p]:
                ok = False
                break
        if not ok:
            continue
            
        for s in range(0, n - p):
            good = True
            for i in range(s, n - p):
                if seq[i] != seq[i + p]:
                    good = False
                    break
            if good:
                return {'period': p, 'start': s, 'seq': seq}
    return None

def build_grundy_data(max_n):
    g, max_g = compute_grundy(max_n)
    even = [g[2 * k] for k in range(1, max_n // 2 + 1)]
    odd = [g[2 * k - 1] for k in range(1, (max_n + 1) // 2 + 1)]

    even_info = detect_period(even)
    odd_info = detect_period(odd)
    
    return {'g': g, 'max_g': max_g, 'even': even_info, 'odd': odd_info}

def add_counts(info, total, counts):
    if total <= 0:
        return
    seq = info['seq']
    start = info['start']
    period = info['period']

    if total <= start:
        for i in range(total):
            counts[seq[i]] += 1
        return

    for i in range(start):
        counts[seq[i]] += 1

    rem = total - start
    freq = [0] * len(counts)
    for i in range(period):
        freq[seq[start + i]] += 1
        
    cycles = rem // period
    tail = rem % period
    
    for gidx in range(len(counts)):
        counts[gidx] += freq[gidx] * cycles
        
    for i in range(tail):
        counts[seq[start + i]] += 1

def grundy_counts(N, data):
    counts = [0] * (data['max_g'] + 1)
    total_even = N // 2
    total_odd = (N + 1) // 2
    
    add_counts(data['even'], total_even, counts)
    add_counts(data['odd'], total_odd, counts)
    return counts

def build_comb(c, m):
    comb = [0] * (m + 1)
    val = 1
    comb[0] = 1
    for k in range(1, m + 1):
        val = (val * (c + k - 1)) // k
        comb[k] = val % MOD
    return comb

def compute_S(N, m, data):
    counts = grundy_counts(N, data)
    
    max_g = data['max_g']
    bits = 0
    while (1 << bits) <= max_g:
        bits += 1
    states = 1 << bits

    dp = [[0] * (m + 1) for _ in range(states)]
    dp[0][0] = 1

    for gval in range(max_g + 1):
        c = counts[gval]
        if c == 0:
            continue

        comb = build_comb(c, m)
        even_terms = []
        odd_terms = []
        
        for k in range(m + 1):
            v = comb[k]
            if v == 0:
                continue
            if k & 1:
                odd_terms.append((k, v))
            else:
                even_terms.append((k, v))

        next_dp = [[0] * (m + 1) for _ in range(states)]
        
        for mask in range(states):
            cur = dp[mask]
            dest_same = next_dp[mask]
            dest_xor = next_dp[mask ^ gval]
            
            for i in range(m + 1):
                curv = cur[i]
                if curv == 0:
                    continue
                    
                for j, v in even_terms:
                    if i + j > m:
                        break
                    dest_same[i + j] = (dest_same[i + j] + curv * v) % MOD
                    
                for j, v in odd_terms:
                    if i + j > m:
                        break
                    dest_xor[i + j] = (dest_xor[i + j] + curv * v) % MOD
                    
        dp = next_dp

    return dp[0][m]

def solve():
    max_n = 20000
    data = build_grundy_data(max_n)
    
    N = 12491249
    m = 1249
    return str(compute_S(N, m, data))

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

Java

import java.math.BigInteger;

public class Euler888 {

    static final long MOD = 912491249L;

    static class PeriodInfo {
        int period = 0;
        int start = 0;
        int[] seq;
    }

    static class GrundyData {
        int[] g;
        int maxG = 0;
        PeriodInfo even;
        PeriodInfo odd;
    }

    static int[] computeGrundy(int maxN, int[] maxGRef) {
        int[] g = new int[maxN + 1];
        int maxG = 0;
        int[] moves = { 1, 2, 4, 9 };

        for (int n = 1; n <= maxN; ++n) {
            long mask = 0;
            for (int mv : moves) {
                if (n >= mv) {
                    mask |= (1L << g[n - mv]);
                }
            }
            int half = n / 2;
            for (int a = 1; a <= half; ++a) {
                mask |= (1L << (g[a] ^ g[n - a]));
            }
            int mex = 0;
            while ((mask & (1L << mex)) != 0)
                ++mex;
            g[n] = mex;
            if (mex > maxG)
                maxG = mex;
        }
        maxGRef[0] = maxG;
        return g;
    }

    static PeriodInfo detectPeriod(int[] seq, int maxPeriod, int minReps) {
        int n = seq.length;
        for (int p = 1; p <= maxPeriod; ++p) {
            int L = p * minReps;
            if (L + p > n)
                continue;
            int startTail = n - L;
            boolean ok = true;
            for (int i = startTail; i + p < n; ++i) {
                if (seq[i] != seq[i + p]) {
                    ok = false;
                    break;
                }
            }
            if (!ok)
                continue;

            for (int s = 0; s + p < n; ++s) {
                boolean good = true;
                for (int i = s; i + p < n; ++i) {
                    if (seq[i] != seq[i + p]) {
                        good = false;
                        break;
                    }
                }
                if (good) {
                    PeriodInfo info = new PeriodInfo();
                    info.period = p;
                    info.start = s;
                    info.seq = seq;
                    return info;
                }
            }
        }
        return new PeriodInfo();
    }

    static GrundyData buildGrundyData(int maxN) {
        GrundyData data = new GrundyData();
        int[] maxGRef = new int[1];
        data.g = computeGrundy(maxN, maxGRef);
        data.maxG = maxGRef[0];

        int[] even = new int[maxN / 2];
        int[] odd = new int[(maxN + 1) / 2];

        int evenIdx = 0, oddIdx = 0;
        for (int k = 1; 2 * k <= maxN; ++k) {
            even[evenIdx++] = data.g[2 * k];
        }
        for (int k = 1; 2 * k - 1 <= maxN; ++k) {
            odd[oddIdx++] = data.g[2 * k - 1];
        }

        data.even = detectPeriod(even, 300, 6);
        data.odd = detectPeriod(odd, 300, 6);
        return data;
    }

    static void addCounts(PeriodInfo info, long total, long[] counts) {
        if (total <= 0)
            return;
        int[] seq = info.seq;
        int start = info.start;
        int period = info.period;

        if (total <= start) {
            for (int i = 0; i < total; ++i) {
                counts[seq[i]]++;
            }
            return;
        }

        for (int i = 0; i < start; ++i) {
            counts[seq[i]]++;
        }

        long rem = total - start;
        long[] freq = new long[counts.length];
        for (int i = 0; i < period; ++i) {
            freq[seq[start + i]]++;
        }

        long cycles = rem / period;
        long tail = rem % period;

        for (int g = 0; g < counts.length; ++g) {
            counts[g] += freq[g] * cycles;
        }
        for (int i = 0; i < tail; ++i) {
            counts[seq[start + i]]++;
        }
    }

    static long[] grundyCounts(long N, GrundyData data) {
        long[] counts = new long[data.maxG + 1];
        long totalEven = N / 2;
        long totalOdd = (N + 1) / 2;
        addCounts(data.even, totalEven, counts);
        addCounts(data.odd, totalOdd, counts);
        return counts;
    }

    static long[] buildComb(long c, int m) {
        long[] comb = new long[m + 1];
        BigInteger val = BigInteger.ONE;
        comb[0] = 1;
        BigInteger bMod = BigInteger.valueOf(MOD);
        for (int k = 1; k <= m; ++k) {
            val = val.multiply(BigInteger.valueOf(c + k - 1));
            val = val.divide(BigInteger.valueOf(k));
            comb[k] = val.mod(bMod).longValue();
        }
        return comb;
    }

    static class Term {
        int k;
        long v;

        Term(int k, long v) {
            this.k = k;
            this.v = v;
        }
    }

    static long computeS(long N, int m, GrundyData data) {
        long[] counts = grundyCounts(N, data);

        int bits = 0;
        while ((1 << bits) <= data.maxG)
            ++bits;
        int states = 1 << bits;

        long[][] dp = new long[states][m + 1];
        dp[0][0] = 1;

        for (int gval = 0; gval <= data.maxG; ++gval) {
            long c = counts[gval];
            if (c == 0)
                continue;

            long[] comb = buildComb(c, m);
            java.util.List<Term> evenTerms = new java.util.ArrayList<>();
            java.util.List<Term> oddTerms = new java.util.ArrayList<>();

            for (int k = 0; k <= m; ++k) {
                long v = comb[k];
                if (v == 0)
                    continue;
                if ((k & 1) == 1)
                    oddTerms.add(new Term(k, v));
                else
                    evenTerms.add(new Term(k, v));
            }

            long[][] nextDp = new long[states][m + 1];

            for (int mask = 0; mask < states; ++mask) {
                long[] cur = dp[mask];
                long[] destSame = nextDp[mask];
                long[] destXor = nextDp[mask ^ gval];

                for (int i = 0; i <= m; ++i) {
                    long curv = cur[i];
                    if (curv == 0)
                        continue;

                    for (Term term : evenTerms) {
                        int j = term.k;
                        if (i + j > m)
                            break;
                        long add = (curv * term.v) % MOD;
                        destSame[i + j] = (destSame[i + j] + add) % MOD;
                    }
                    for (Term term : oddTerms) {
                        int j = term.k;
                        if (i + j > m)
                            break;
                        long add = (curv * term.v) % MOD;
                        destXor[i + j] = (destXor[i + j] + add) % MOD;
                    }
                }
            }
            dp = nextDp;
        }

        return dp[0][m];
    }

    public static String solve() {
        GrundyData data = buildGrundyData(20000);
        return Long.toString(computeS(12491249L, 1249, data));
    }

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