WARNING: The text below provides a guidance for a Project Euler problem.

  • If you pride yourself in tackling problems entirely on your own, you are encouraged to exit this post right now.

  • If you are simply looking for a hint for the problem, please visit the »Project Euler’s official forum«. You may find some hints here, but if you don’t want your problem spoiled, scroll cautiously.

  • If you are looking for a readily available copy/paste answer, you will not find it here. All of the code snippets in the post have been stripped of some crucial lines (these are clearly marked). It is my aim that anyone able to replicate the answer based on such snippets will have to understand the solution conceptually. At that point, by virtue of learning something new, I believe he deserves it anyway.

  • Guidance for future problems will not be published before 100 people have solved the problem and there are at least two more recent problems.

Scatterstone Nim

The problem tasks us with finding the number of winning positions in a certain variation of »Nim« defined by parameters . We will, in the beginning, focus on one particular given pair . To reiterate, the game goes like this:

  • Two players start the game with some number of piles of stones in total.
  • They take turns, where they pick one pile and split it into nonzero piles.
  • Whoever doesn’t have any valid move anymore, loses.

This being a combinatorial impartial game, the »Sprague-Grundy theorem« guarantees that we will be able to analyse it with »Grundy numbers (nimbers)«. Consider the set of (sorted) tuples with ; these will represent all possible positions in our game ( is the size of -th pile). Each position will be associated with its Grundy number, which is standardly defined recursively as:

The in the above expression stands for the minimal excludant: . The standard combinatorial game theory tells us that position will be losing if and only if , and moreover, we can understand playing the position as playing the games in positions simultaneously. The nimber addidion gives us the relation ( denotes the logical xor):

The calculation of Grundy numbers using is quite expensive, but the above tells us that we may find any readily just by precalculating all for . For the precalculation of these singleton positions, we shall use . But before dwelling unnecessarily into brute force, there are helpful simplifications worth noting:

Theorem: If , then .

Proof. Comes from the fact that we can split an even pile only into two piles of equal parity, and that we can split an odd pile only into two piles of distinct parity.

Theorem: If , then .

Proof. The claim holds for , since . We will now show by induction that if we have a pile of size , positions with Grundy numbers are reachable. In particular, we may perform splits (for any valid ):

Thus the only case for where we actually need to compute the Grundy numbers is . In that case, letting denote the set of integer partitions of of cardinality , we use:

def P(n, k, l=1):
    if k < 1:
        raise StopIteration
    if k == 1:
        if n >= l:
            yield (n,)
        raise StopIteration
    for i in range(l,n+1):
        for result in P(n-i,k-1,i):
            # <REMOVED>

grundy3 = [0, 0]
for n in range(2, maxn):
    grundy3.append(mex([xor([grundy3[i] for i in t]) for t in P(n, 3)] + [xor([grundy3[i] for i in t]) # <REMOVED>

def grundy(n, k):
    if k == 2:
        return 0 if n&1 else 1
    if k == 3:
        return grundy3[n]
        return # <REMOVED>

Once we have computed all Grundy numbers, it will be quite easy to write a recursive function to count the number of winning positions; these are just the ones with Grundy number . Specifically, we are looking for the number of integer partitions of such that the xor sum of the Grundy numbers of components is nonzero.

def F(n, k, M, xor):
    if n == 0:
        return int(xor != 0)
    res = 0
    for i in range(1, min(M, n) + 1):
        res += # <REMOVED>
        res %= MOD
    return res

def f(n, k):
    return F(n, k, n, 0)
def g(n):
    res = 0
    for k in range(2, n+1):
        res += f(n, min(k, 4))
        res %= MOD
    return res

This piece of Python code returns the answer in about 9 seconds when run using PyPy.