Problem 550: Divisor Game
View on Project EulerProject Euler Problem 550 Solution
EulerSolve provides an optimized solution for Project Euler Problem 550, Divisor Game, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary We consider a game with \(k\) independent piles. Each pile starts at an integer \(x\) chosen from the interval \([2,n]\). For one pile, a legal move selects two proper divisors of the current number, both strictly between \(1\) and the number itself, and replaces that pile by the two resulting subgames. By the Sprague-Grundy theorem, the whole position is losing exactly when the xor of the pile Grundy values is \(0\). The task is therefore to count how many \(k\)-tuples \((x_1,\dots,x_k)\in [2,n]^k\) are winning, with the final answer taken modulo $$987654321.$$ Mathematical Approach The program splits the work into two layers: first compute the Grundy value of one pile, then combine those one-pile values across \(k\) piles by xor convolution. Step 1: Reduce a Number to Its Prime-Exponent Pattern Write $$x=\prod_{i=1}^{m} p_i^{a_i},\qquad a_1\ge a_2\ge \cdots \ge a_m\ge 1.$$ The ordered list \(\lambda(x)=(a_1,\dots,a_m)\) is the canonical state used by the implementation. The actual primes do not matter; only the multiset of exponents matters. Why is this valid? Every divisor of \(x\) is obtained by choosing exponents $$0\le b_i\le a_i$$ and forming \(\prod p_i^{b_i}\). After removing zero exponents and sorting the remaining positive ones in non-increasing order, we get the same canonical pattern no matter which prime labels were used....
Detailed mathematical approach
Problem Summary
We consider a game with \(k\) independent piles. Each pile starts at an integer \(x\) chosen from the interval \([2,n]\). For one pile, a legal move selects two proper divisors of the current number, both strictly between \(1\) and the number itself, and replaces that pile by the two resulting subgames. By the Sprague-Grundy theorem, the whole position is losing exactly when the xor of the pile Grundy values is \(0\).
The task is therefore to count how many \(k\)-tuples \((x_1,\dots,x_k)\in [2,n]^k\) are winning, with the final answer taken modulo
$$987654321.$$
Mathematical Approach
The program splits the work into two layers: first compute the Grundy value of one pile, then combine those one-pile values across \(k\) piles by xor convolution.
Step 1: Reduce a Number to Its Prime-Exponent Pattern
Write
$$x=\prod_{i=1}^{m} p_i^{a_i},\qquad a_1\ge a_2\ge \cdots \ge a_m\ge 1.$$
The ordered list \(\lambda(x)=(a_1,\dots,a_m)\) is the canonical state used by the implementation. The actual primes do not matter; only the multiset of exponents matters.
Why is this valid? Every divisor of \(x\) is obtained by choosing exponents
$$0\le b_i\le a_i$$
and forming \(\prod p_i^{b_i}\). After removing zero exponents and sorting the remaining positive ones in non-increasing order, we get the same canonical pattern no matter which prime labels were used. Therefore any two integers with the same exponent pattern have exactly the same family of reachable divisor patterns. By induction on the total exponent sum, they also have the same Grundy value.
For example,
$$12=2^2\cdot 3,\qquad 18=2\cdot 3^2$$
both have pattern \((2,1)\), so they belong to the same one-pile game state.
Step 2: Derive the One-Pile Grundy Recurrence
Fix a pattern \(\lambda\). Let \(\mathcal{D}(\lambda)\) be the set of canonical patterns obtained from all proper divisors of a number with pattern \(\lambda\), excluding \(1\) and excluding the number itself.
If \(\mathcal{G}(\mu)\) denotes the Grundy value of pattern \(\mu\), define
$$\Sigma(\lambda)=\{\mathcal{G}(\mu):\mu\in\mathcal{D}(\lambda)\}.$$
A legal move chooses two proper divisors, so the resulting position is the disjoint sum of two smaller games. Sprague-Grundy theory therefore gives the option nimbers
$$u\oplus v\qquad (u,v\in \Sigma(\lambda)),$$
where \(\oplus\) denotes bitwise xor and the two chosen divisors may be equal. Hence
$$\boxed{\mathcal{G}(\lambda)=\operatorname{mex}\{u\oplus v:\ u,v\in\Sigma(\lambda)\}.}$$
If \(\mathcal{D}(\lambda)\) is empty, there is no legal move and the Grundy value is \(0\).
This formula is exactly what the implementations evaluate recursively and memoize by pattern.
Step 3: Enumerate Proper Divisors in Exponent Space
Suppose \(\lambda=(a_1,\dots,a_m)\). Then every divisor pattern comes from a vector \((b_1,\dots,b_m)\) with
$$0\le b_i\le a_i.$$
The all-zero vector gives the forbidden divisor \(1\), while the vector \((a_1,\dots,a_m)\) gives the forbidden divisor \(x\) itself. All other vectors correspond to proper divisors.
After deleting zeros and sorting, several different exponent choices may collapse to the same canonical pattern, but that is harmless because only the resulting Grundy values matter. The recursion is finite because every proper divisor has strictly smaller total exponent sum than the original state.
Small examples:
$$\lambda=(1)\implies \mathcal{D}(\lambda)=\varnothing,\qquad \mathcal{G}(1)=0,$$
$$\lambda=(2)\implies \mathcal{D}(\lambda)=\{(1)\},\qquad \mathcal{G}(2)=\operatorname{mex}\{0\}=1,$$
$$\lambda=(3)\implies \Sigma(\lambda)=\{0,1\},\qquad \mathcal{G}(3)=\operatorname{mex}\{0,1\}=2.$$
Likewise the squarefree pattern \((1,1)\) has only prime proper divisors, so its option set again has nimber set \(\{0\}\), giving Grundy value \(1\).
Step 4: Convert One-Pile Values into a Frequency Distribution
Once the one-pile Grundy value is known for each \(x\in[2,n]\), define the histogram
$$C_r=\#\{x\in[2,n]:\mathcal{G}(x)=r\}.$$
For \(k\) independent piles, the number of positions whose total xor equals \(t\) is the \(k\)-fold xor convolution of the sequence \(C\):
$$F=\underbrace{C *_\oplus C *_\oplus \cdots *_\oplus C}_{k\text{ factors}},$$
where
$$ (A *_\oplus B)_t=\sum_{i\oplus j=t} A_iB_j. $$
The losing positions are exactly those with xor \(0\), so the desired answer is
$$\boxed{(n-1)^k-F_0 \pmod{987654321}.}$$
Step 5: Diagonalize XOR Convolution with the Walsh-Hadamard Transform
Xor convolution is handled efficiently by the Walsh-Hadamard transform. If \(\widehat{C}\) denotes the transform of \(C\), then
$$\widehat{F}_j=\bigl(\widehat{C}_j\bigr)^k.$$
So the procedure is:
$$C\ \longrightarrow\ \widehat{C}\ \longrightarrow\ \widehat{F}\ \longrightarrow\ F.$$
The inverse transform divides by the transform length. Because the modulus \(987654321\) is odd, every power of two has a modular inverse, so the inverse transform is valid modulo the required modulus.
Worked Example: The Checkpoint \((n,k)=(10,5)\)
The single-pile values for \(x\in[2,10]\) are easy to derive from the recurrence:
$$\begin{aligned} &2,3,5,7 &&\text{are prime} &&\Rightarrow \mathcal{G}=0,\\ &4,6,9,10 &&\text{have divisor nimber set }\{0\} &&\Rightarrow \mathcal{G}=1,\\ &8 &&\text{has divisor nimber set }\{0,1\} &&\Rightarrow \mathcal{G}=2. \end{aligned}$$
Therefore
$$C=(4,4,1,0).$$
Using transform length \(4\), the Walsh-Hadamard transform is
$$\widehat{C}=(9,1,7,-1).$$
Raise componentwise to the fifth power:
$$\widehat{F}=(9^5,1^5,7^5,(-1)^5)=(59049,1,16807,-1).$$
The zero-xor coefficient after the inverse transform is
$$F_0=\frac{59049+1+16807-1}{4}=18964.$$
Total positions are
$$9^5=59049,$$
so the number of winning positions is
$$59049-18964=40085,$$
which is the checkpoint used by the implementation.
How the Code Works
The C++, Python, and Java implementations follow the same pipeline. First they build a smallest-prime-factor sieve up to \(n\). Each integer from \(2\) to \(n\) is factored with that sieve, its prime exponents are sorted into canonical non-increasing order, and the corresponding one-pile Grundy value is obtained by memoized recursion on patterns.
During that recursion, the implementation enumerates every proper divisor in exponent space, converts it to its canonical pattern, gathers the Grundy values of all reachable proper-divisor states, and computes the mex of all pairwise xor combinations. The resulting one-pile Grundy histogram is then padded to the next power of two, transformed with the xor Walsh-Hadamard transform, raised pointwise to the \(k\)-th power, and inverse transformed. The coefficient of xor \(0\) is subtracted from \((n-1)^k\), always modulo \(987654321\).
Complexity Analysis
The smallest-prime-factor sieve up to \(n\) is linear or near-linear in practice and uses \(O(n)\) memory. Factoring all integers from \(2\) to \(n\) is near-linear on average. The recursive Grundy phase depends on the number of distinct exponent patterns actually encountered and on the number of divisor exponent vectors generated for each pattern; for \(n=10^7\), those patterns are heavily reused and each exponent list is short.
If the largest observed Grundy value is below a power of two \(L\), then the xor transform works on length \(L\) and costs \(O(L\log L)\) time, plus \(O(L\log k)\) for the pointwise modular exponentiations, with \(O(L)\) additional memory. This is what makes the \(k\)-pile aggregation feasible even when \(k\) is extremely large.
Footnotes and References
- Problem page: https://projecteuler.net/problem=550
- Sprague-Grundy theorem: Wikipedia — Sprague-Grundy theorem
- Prime factorization: Wikipedia — Prime factorization
- Divisor function and divisor structure: Wikipedia — Divisor
- Walsh-Hadamard transform and xor convolution: cp-algorithms — Walsh-Hadamard transform
Problem 550 source code
C++
#include <algorithm>
#include <cassert>
#include <cstdint>
#include <functional>
#include <iostream>
#include <numeric>
#include <unordered_map>
#include <vector>
namespace {
using u32 = std::uint32_t;
using u64 = std::uint64_t;
using i64 = std::int64_t;
static constexpr u64 MOD = 987654321ULL;
static u64 mod_pow(u64 a, u64 e) {
a %= MOD;
u64 r = 1 % MOD;
while (e) {
if (e & 1ULL) r = (r * a) % MOD;
a = (a * a) % MOD;
e >>= 1ULL;
}
return r;
}
static i64 ext_gcd(i64 a, i64 b, i64& x, i64& y) {
if (b == 0) {
x = 1;
y = 0;
return a;
}
i64 x1 = 0, y1 = 0;
i64 g = ext_gcd(b, a % b, x1, y1);
x = y1;
y = x1 - (a / b) * y1;
return g;
}
static u64 mod_inv(u64 a) {
i64 x = 0, y = 0;
i64 g = ext_gcd(static_cast<i64>(a % MOD), static_cast<i64>(MOD), x, y);
assert(g == 1);
i64 res = x % static_cast<i64>(MOD);
if (res < 0) res += static_cast<i64>(MOD);
return static_cast<u64>(res);
}
static void fwht_xor(std::vector<u64>& a, bool inverse) {
const std::size_t n = a.size();
for (std::size_t len = 1; len < n; len <<= 1U) {
for (std::size_t i = 0; i < n; i += (len << 1U)) {
for (std::size_t j = 0; j < len; ++j) {
u64 u = a[i + j];
u64 v = a[i + j + len];
u64 x = u + v;
if (x >= MOD) x -= MOD;
u64 y = (u >= v) ? (u - v) : (u + MOD - v);
a[i + j] = x;
a[i + j + len] = y;
}
}
}
if (inverse) {
const u64 inv_n = mod_inv(static_cast<u64>(n));
for (u64& x : a) x = (x * inv_n) % MOD;
}
}
static u64 encode_pattern(const int* exps, int len) {
// Encode non-increasing exponent multiset into 64-bit:
// top 4 bits: length (<=15), then 6 bits per exponent (<=63).
u64 key = (static_cast<u64>(len) << 60U);
for (int i = 0; i < len; ++i) key |= (static_cast<u64>(exps[i]) << (6U * i));
return key;
}
static void decode_pattern(u64 key, int* exps, int& len) {
len = static_cast<int>(key >> 60U);
for (int i = 0; i < len; ++i) exps[i] = static_cast<int>((key >> (6U * i)) & 63ULL);
}
struct GrundyComputer {
std::unordered_map<u64, u32> memo;
GrundyComputer() { memo.reserve(1024); }
u32 grundy(u64 key) {
auto it = memo.find(key);
if (it != memo.end()) return it->second;
int a[12];
int len = 0;
decode_pattern(key, a, len);
// Enumerate all proper divisors d of n (excluding 1 and n itself) in exponent space.
int cur[12];
std::vector<u32> div_nims;
div_nims.reserve(256);
std::function<void(int, bool, bool)> dfs = [&](int idx, bool any_pos, bool all_eq) {
if (idx == len) {
if (!any_pos) return; // d=1
if (all_eq) return; // d=n
int b[12];
int blen = 0;
for (int i = 0; i < len; ++i) {
if (cur[i] > 0) b[blen++] = cur[i];
}
// insertion sort descending, blen is small
for (int i = 1; i < blen; ++i) {
int v = b[i];
int j = i;
while (j > 0 && b[j - 1] < v) {
b[j] = b[j - 1];
--j;
}
b[j] = v;
}
const u64 dkey = encode_pattern(b, blen);
div_nims.push_back(grundy(dkey));
return;
}
for (int e = 0; e <= a[idx]; ++e) {
cur[idx] = e;
dfs(idx + 1, any_pos || (e > 0), all_eq && (e == a[idx]));
}
};
dfs(0, false, true);
if (div_nims.empty()) {
memo.emplace(key, 0U);
return 0U;
}
std::sort(div_nims.begin(), div_nims.end());
div_nims.erase(std::unique(div_nims.begin(), div_nims.end()), div_nims.end());
const u32 mx = div_nims.back();
std::vector<std::uint8_t> present(static_cast<std::size_t>(mx + 1), 0);
for (u32 v : div_nims) present[static_cast<std::size_t>(v)] = 1;
u32 mex = 0;
while (true) {
bool ok = false;
for (u32 v : div_nims) {
const u32 t = v ^ mex;
if (t <= mx && present[static_cast<std::size_t>(t)]) {
ok = true;
break;
}
}
if (!ok) break;
++mex;
}
memo.emplace(key, mex);
return mex;
}
};
static std::vector<u32> build_spf(u32 n) {
std::vector<u32> spf(static_cast<std::size_t>(n + 1), 0U);
std::vector<u32> primes;
primes.reserve(static_cast<std::size_t>(n / 10));
for (u32 i = 2; i <= n; ++i) {
if (spf[i] == 0) {
spf[i] = i;
primes.push_back(i);
}
for (u32 p : primes) {
const u64 v = static_cast<u64>(p) * static_cast<u64>(i);
if (v > n) break;
spf[static_cast<std::size_t>(v)] = p;
if (p == spf[i]) break;
}
}
return spf;
}
static u64 winning_positions(u32 n, u64 k, GrundyComputer& gc) {
const std::vector<u32> spf = build_spf(n);
std::vector<u64> counts(1, 0ULL);
u32 max_g = 0;
for (u32 x = 2; x <= n; ++x) {
u32 t = x;
int exps[12];
int len = 0;
while (t > 1) {
const u32 p = spf[t];
int c = 0;
while (t > 1 && spf[t] == p) {
t /= p;
++c;
}
exps[len++] = c;
}
// insertion sort descending
for (int i = 1; i < len; ++i) {
int v = exps[i];
int j = i;
while (j > 0 && exps[j - 1] < v) {
exps[j] = exps[j - 1];
--j;
}
exps[j] = v;
}
const u64 key = encode_pattern(exps, len);
const u32 g = gc.grundy(key);
if (g >= counts.size()) counts.resize(static_cast<std::size_t>(g + 1), 0ULL);
++counts[static_cast<std::size_t>(g)];
if (g > max_g) max_g = g;
}
std::size_t sz = 1;
while (sz <= static_cast<std::size_t>(max_g)) sz <<= 1U;
std::vector<u64> a(sz, 0ULL);
for (std::size_t i = 0; i < counts.size(); ++i) a[i] = counts[i] % MOD;
fwht_xor(a, false);
for (u64& x : a) x = mod_pow(x, k);
fwht_xor(a, true);
const u64 losing = a[0] % MOD;
const u64 total = mod_pow(static_cast<u64>(n - 1U), k);
return (total + MOD - losing) % MOD;
}
} // namespace
int main() {
std::ios::sync_with_stdio(false);
std::cin.tie(nullptr);
GrundyComputer gc;
// Validation point from the problem statement.
assert(winning_positions(10, 5, gc) == 40085ULL);
std::cout << winning_positions(10'000'000U, 1'000'000'000'000ULL, gc) % MOD << '\n';
return 0;
}
Python
def solve():
MOD = 987654321
N = 10_000_000
K = 1_000_000_000_000
def mod_pow(a, e, m=MOD):
a %= m
r = 1 % m
while e:
if e & 1:
r = r * a % m
a = a * a % m
e >>= 1
return r
def mod_inv(a, m=MOD):
b = m
u, v = 1, 0
while b:
t = a // b
a, b = b, a - t * b
u, v = v, u - t * v
return u % m
def fwht_xor(a, inverse=False):
n = len(a)
ln = 1
while ln < n:
for i in range(0, n, ln << 1):
for j in range(ln):
u = a[i + j]
v = a[i + j + ln]
x = u + v
if x >= MOD:
x -= MOD
y = u - v
if y < 0:
y += MOD
a[i + j] = x
a[i + j + ln] = y
ln <<= 1
if inverse:
iv = mod_inv(n)
for i in range(n):
a[i] = a[i] * iv % MOD
def encode_pattern(exps):
key = len(exps) << 60
for i, e in enumerate(exps):
key |= e << (6 * i)
return key
def decode_pattern(key):
ln = key >> 60
exps = []
for i in range(ln):
exps.append((key >> (6 * i)) & 63)
return exps
memo = {}
def grundy(key):
if key in memo:
return memo[key]
a = decode_pattern(key)
ln = len(a)
div_nims = []
def dfs(idx, any_pos, all_eq, cur):
if idx == ln:
if not any_pos:
return
if all_eq:
return
b = sorted([c for c in cur if c > 0], reverse=True)
dkey = encode_pattern(b)
div_nims.append(grundy(dkey))
return
for e in range(a[idx] + 1):
cur.append(e)
dfs(idx + 1, any_pos or e > 0, all_eq and e == a[idx], cur)
cur.pop()
dfs(0, False, True, [])
if not div_nims:
memo[key] = 0
return 0
div_nims = sorted(set(div_nims))
present = set(div_nims)
mex = 0
while True:
ok = False
for v in div_nims:
t = v ^ mex
if t in present:
ok = True
break
if not ok:
break
mex += 1
memo[key] = mex
return mex
# Build SPF
spf = list(range(N + 1))
primes = []
for i in range(2, N + 1):
if spf[i] == i:
primes.append(i)
for p in primes:
if i * p > N:
break
spf[i * p] = p
if p == spf[i]:
break
counts = {}
max_g = 0
for x in range(2, N + 1):
t = x
exps = []
while t > 1:
p = spf[t]
c = 0
while t > 1 and spf[t] == p:
t //= p
c += 1
exps.append(c)
exps.sort(reverse=True)
key = encode_pattern(exps)
g = grundy(key)
counts[g] = counts.get(g, 0) + 1
if g > max_g:
max_g = g
sz = 1
while sz <= max_g:
sz <<= 1
a = [0] * sz
for g, cnt in counts.items():
a[g] = cnt % MOD
fwht_xor(a)
for i in range(sz):
a[i] = mod_pow(a[i], K)
fwht_xor(a, True)
losing = a[0] % MOD
total = mod_pow(N - 1, K)
return str((total + MOD - losing) % MOD)
if __name__ == '__main__':
print(solve())
Java
import java.util.ArrayList;
import java.util.Arrays;
import java.util.Collections;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
public class Euler550 {
static final long MOD = 987654321L;
static long modPow(long a, long e) {
a %= MOD;
long r = 1 % MOD;
while (e > 0) {
if ((e & 1) != 0)
r = (r * a) % MOD;
a = (a * a) % MOD;
e >>= 1;
}
return r;
}
static long[] extGcd(long a, long b) {
if (b == 0)
return new long[] { a, 1, 0 };
long[] res = extGcd(b, a % b);
long g = res[0];
long x1 = res[1];
long y1 = res[2];
return new long[] { g, y1, x1 - (a / b) * y1 };
}
static long modInv(long a) {
long[] res = extGcd(a % MOD, MOD);
long x = res[1] % MOD;
if (x < 0)
x += MOD;
return x;
}
static void fwhtXor(long[] a, boolean inverse) {
int n = a.length;
for (int len = 1; len < n; len <<= 1) {
for (int i = 0; i < n; i += (len << 1)) {
for (int j = 0; j < len; j++) {
long u = a[i + j];
long v = a[i + j + len];
long x = u + v;
if (x >= MOD)
x -= MOD;
long y = u - v;
if (y < 0)
y += MOD;
a[i + j] = x;
a[i + j + len] = y;
}
}
}
if (inverse) {
long invN = modInv(n);
for (int i = 0; i < n; i++) {
a[i] = (a[i] * invN) % MOD;
}
}
}
static class GrundyComputer {
Map<Long, Integer> memo = new HashMap<>();
int grundy(long key) {
if (memo.containsKey(key))
return memo.get(key);
List<Integer> aList = new ArrayList<>();
long tempKey = key;
int aLen = (int) (tempKey >>> 60);
for (int i = 0; i < aLen; i++) {
aList.add((int) ((tempKey >> (6 * i)) & 63L));
}
int[] a = new int[aLen];
for (int i = 0; i < aLen; i++)
a[i] = aList.get(i);
int[] cur = new int[aLen];
List<Integer> divNims = new ArrayList<>();
dfs(0, false, true, a, cur, divNims);
if (divNims.isEmpty()) {
memo.put(key, 0);
return 0;
}
Collections.sort(divNims);
List<Integer> uniqueNims = new ArrayList<>();
for (int v : divNims) {
if (uniqueNims.isEmpty() || uniqueNims.get(uniqueNims.size() - 1) != v) {
uniqueNims.add(v);
}
}
int mx = uniqueNims.get(uniqueNims.size() - 1);
boolean[] present = new boolean[mx + 1];
for (int v : uniqueNims)
present[v] = true;
int mex = 0;
while (true) {
boolean ok = false;
for (int v : uniqueNims) {
int t = v ^ mex;
if (t <= mx && present[t]) {
ok = true;
break;
}
}
if (!ok)
break;
mex++;
}
memo.put(key, mex);
return mex;
}
void dfs(int idx, boolean anyPos, boolean allEq, int[] a, int[] cur, List<Integer> divNims) {
if (idx == a.length) {
if (!anyPos)
return;
if (allEq)
return;
int[] b = new int[a.length];
int blen = 0;
for (int i = 0; i < a.length; i++) {
if (cur[i] > 0)
b[blen++] = cur[i];
}
for (int i = 1; i < blen; i++) {
int v = b[i];
int j = i;
while (j > 0 && b[j - 1] < v) {
b[j] = b[j - 1];
j--;
}
b[j] = v;
}
long dkey = ((long) blen) << 60;
for (int i = 0; i < blen; i++) {
dkey |= ((long) b[i]) << (6 * i);
}
divNims.add(grundy(dkey));
return;
}
for (int e = 0; e <= a[idx]; e++) {
cur[idx] = e;
dfs(idx + 1, anyPos || (e > 0), allEq && (e == a[idx]), a, cur, divNims);
}
}
}
static int[] buildSpf(int n) {
int[] spf = new int[n + 1];
int numPrimes = 0;
int maxPrimes = Math.max(1000, (int) (n / Math.max(1, Math.log(n)) * 1.2));
int[] primes = new int[maxPrimes];
for (int i = 2; i <= n; i++) {
if (spf[i] == 0) {
spf[i] = i;
if (numPrimes < primes.length) {
primes[numPrimes++] = i;
} else {
int[] newPrimes = new int[primes.length * 2];
System.arraycopy(primes, 0, newPrimes, 0, primes.length);
primes = newPrimes;
primes[numPrimes++] = i;
}
}
for (int j = 0; j < numPrimes; j++) {
int p = primes[j];
long v = (long) p * i;
if (v > n)
break;
spf[(int) v] = p;
if (p == spf[i])
break;
}
}
return spf;
}
static long winningPositions(int n, long k, GrundyComputer gc) {
int[] spf = buildSpf(n);
long[] counts = new long[1];
int maxG = 0;
int[] exps = new int[32];
for (int x = 2; x <= n; x++) {
int t = x;
int len = 0;
while (t > 1) {
int p = spf[t];
int c = 0;
while (t > 1 && spf[t] == p) {
t /= p;
c++;
}
exps[len++] = c;
}
for (int i = 1; i < len; i++) {
int v = exps[i];
int j = i;
while (j > 0 && exps[j - 1] < v) {
exps[j] = exps[j - 1];
j--;
}
exps[j] = v;
}
long key = ((long) len) << 60;
for (int i = 0; i < len; i++) {
key |= ((long) exps[i]) << (6 * i);
}
int g = gc.grundy(key);
if (g >= counts.length) {
long[] newCounts = new long[g + 1];
System.arraycopy(counts, 0, newCounts, 0, counts.length);
counts = newCounts;
}
counts[g]++;
if (g > maxG)
maxG = g;
}
int sz = 1;
while (sz <= maxG)
sz <<= 1;
long[] a = new long[sz];
for (int i = 0; i < counts.length; i++) {
a[i] = counts[i] % MOD;
}
fwhtXor(a, false);
for (int i = 0; i < a.length; i++) {
a[i] = modPow(a[i], k);
}
fwhtXor(a, true);
long losing = a[0] % MOD;
long total = modPow(n - 1, k);
long ans = (total + MOD - losing) % MOD;
return ans;
}
public static String solve() {
GrundyComputer gc = new GrundyComputer();
return Long.toString(winningPositions(10000000, 1000000000000L, gc));
}
public static void main(String[] args) {
System.out.println(solve());
}
}