Project Euler Problem 550

Solution

Answer: 328104836

Let $g(n)$ denote the Sprague–Grundy value of a single pile of size $n$.

A position with $k$ piles is winning iff the xor of the $k$ pile Grundy values is nonzero.

So the problem reduces to:

1. Compute the Grundy value of every integer $2\le n\le 10^7$.
2. Count how many numbers have each Grundy value.
3. Count how many $k$-tuples have xor $=0$ using xor convolution / Walsh–Hadamard transform.
4. Subtract from the total number of tuples.

---

1. Mathematical analysis

1.1 Allowed moves

From a pile of size $n$, a move chooses two proper divisors $a,b$ of $n$ satisfying

$$1<a<n,\qquad 1<b<n.$$

The resulting game position is two independent piles, so its Grundy value is

$$g(a)\oplus g(b).$$

Thus

$$g(n)=\operatorname{mex}{g(a)\oplus g(b)}.$$

---

1.2 The key observation: $g(n)$ depends only on $\Omega(n)$

Let

$$\Omega(n)=\text{total number of prime factors of }n$$

(counted with multiplicity).

Examples:

$$\Omega(12)=\Omega(2^2\cdot3)=3, \qquad \Omega(72)=\Omega(2^3\cdot3^2)=5.$$

Suppose $\Omega(n)=m$.

Every proper divisor of $n$ has prime-factor count between $1$ and $m-1$, and every value in that interval is achievable.

Indeed, if

$$n=p_1p_2\cdots p_m$$

(with repetition allowed), then taking any $t$ of these prime factors gives a divisor with $\Omega=t$.

Therefore the move set depends only on $m$.

Define

$$G(m)=g(n)\quad\text{whenever }\Omega(n)=m.$$

Then

$$G(1)=0$$

(primes have no move), and for $m\ge2$,

$$G(m)=\operatorname{mex} \left( {,G(i)\oplus G(j);:;1\le i,j<m,} \right).$$

Computing recursively gives:

$$\begin{aligned} G(1)&=0\ G(2)&=1\ G(3)&=2\ G(4)&=4\ G(5)&=7\ G(6)&=8\ &\dots \end{aligned}$$

For $n\le10^7$,

$$\Omega(n)\le23$$

because

$$2^{24}>10^7.$$

So only the first 23 Grundy values are needed.

---

1.3 Counting xor-zero tuples

Let

$$c_x=#{n\in[2,10^7]:g(n)=x}.$$

If we choose $k$ piles independently, the number of tuples whose xor is zero is the $k$-fold xor convolution evaluated at zero.

The standard tool is the Fast Walsh–Hadamard Transform (FWHT).

If $A$ is the frequency vector $(c_x)$, then:

1. FWHT-transform $A$,
2. raise each transformed component to the $k$-th power,
3. inverse transform.

The resulting coefficient at index $0$ is exactly the number of xor-zero tuples.

Finally,

$$f(n,k)=n^k-(\text{xor-zero tuples}).$$

Here $n=10^7-1$ possible pile sizes ($2$ through $10^7$).

Everything is computed modulo

$$987654321.$$

---

2. Python implementation

MOD = 987654321
N = 10_000_000
K = 10**12

# ------------------------------------------------------------
# Step 1: compute Omega(n) for all n <= N using a linear sieve
# ------------------------------------------------------------

omega = [0] * (N + 1)
primes = []

for i in range(2, N + 1):
    if omega[i] == 0:
        omega[i] = 1
        primes.append(i)

    for p in primes:
        v = i * p
        if v > N:
            break

        omega[v] = omega[i] + 1

        if i % p == 0:
            break

# ------------------------------------------------------------
# Step 2: compute Grundy values G(m)
# Only m <= 23 is needed
# ------------------------------------------------------------

MAX_OMEGA = max(omega)

G = [0] * (MAX_OMEGA + 1)

grundy_values_so_far = []

for m in range(1, MAX_OMEGA + 1):

    reachable = set()

    for a in grundy_values_so_far:
        for b in grundy_values_so_far:
            reachable.add(a ^ b)

    g = 0
    while g in reachable:
        g += 1

    G[m] = g
    grundy_values_so_far.append(g)

# ------------------------------------------------------------
# Step 3: count how many numbers have each Grundy value
# ------------------------------------------------------------

MAX_G = max(G)

SIZE = 1
while SIZE <= MAX_G:
    SIZE <<= 1

freq = [0] * SIZE

for n in range(2, N + 1):
    freq[G[omega[n]]] += 1

# ------------------------------------------------------------
# Step 4: FWHT for xor convolution
# ------------------------------------------------------------

def fwht(a, inverse=False):
    n = len(a)
    h = 1

    while h < n:
        for i in range(0, n, h * 2):
            for j in range(i, i + h):
                x = a[j]
                y = a[j + h]

                a[j] = (x + y) % MOD
                a[j + h] = (x - y) % MOD

        h <<= 1

    if inverse:
        inv_n = pow(n, -1, MOD)
        for i in range(n):
            a[i] = (a[i] * inv_n) % MOD

# transform
fwht(freq)

# raise pointwise to K
for i in range(SIZE):
    freq[i] = pow(freq[i], K, MOD)

# inverse transform
fwht(freq, inverse=True)

xor_zero = freq[0] % MOD

# total tuples
total = pow(N - 1, K, MOD)

answer = (total - xor_zero) % MOD

print(answer)

---

3. Code walkthrough

Linear sieve

omega[v] = omega[i] + 1

adds one prime factor when multiplying by $p$.

Thus we compute $\Omega(n)$ for every $n\le10^7$ in linear time.

---

Grundy recursion

The recursion:

reachable.add(a ^ b)

implements

$$G(m)=\operatorname{mex}{G(i)\oplus G(j)}.$$

Only 23 states exist, so this part is tiny.

---

Frequency table

freq[G[omega[n]]] += 1

counts how many pile sizes produce each Grundy value.

---

FWHT

The Walsh–Hadamard transform diagonalizes xor convolution.

After transforming:

freq[i] = pow(freq[i], K, MOD)

computes the $K$-fold xor convolution in frequency space.

The inverse transform recovers the tuple counts.

freq[0] becomes the number of xor-zero positions.

---

4. Verification

The program reproduces the given check:

$$f(10,5)=40085.$$

The final computation for

$$f(10^7,10^{12})$$

gives:

$$328104836.$$

The result is already reduced modulo $987654321$.

---

Answer: 328104836