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

1. Problem page: https://projecteuler.net/problem=550
2. Sprague-Grundy theorem: Wikipedia — Sprague-Grundy theorem
3. Prime factorization: Wikipedia — Prime factorization
4. Divisor function and divisor structure: Wikipedia — Divisor
5. Walsh-Hadamard transform and xor convolution: cp-algorithms — Walsh-Hadamard transform