Solution

Answer: 54672965

Let the position be $(x,y)$ with $0<x<y$.

A move replaces $(x,y)$ by $(x,y-kx)$ for some integer $k\ge1$, followed by reordering if necessary.

We seek all losing positions (P-positions), then compute

$$S(N)=\sum_{\substack{(x,y)\text{ losing}\0<x<y\le N}} (x+y).$$

1. Mathematical analysis

Step 1: Characterising the losing positions

Let

$$\varphi=\frac{1+\sqrt5}{2}.$$

We claim:

$$(x,y)\text{ is losing } \iff y<\varphi x.$$

Equivalently,

$$x<y<x+\frac{x}{\varphi},$$

since $\varphi=1+\frac1\varphi$.

Why?

Suppose $y\ge \varphi x$.

Then there exists an integer $k\ge1$ such that

$$y-kx<\varphi x$$

and after reordering the new position is losing.

So every such position is winning.

Conversely, if $y<\varphi x$, then any legal move subtracts $kx$ with $k\ge1$. The new larger pile becomes

$$y-kx.$$

After reordering, one checks that the resulting ratio is always at least $\varphi$, hence every move goes to a winning position.

Therefore the losing positions are exactly

$$x<y<\varphi x.$$

Step 2: Counting losing $y$ for fixed $x$

For fixed $x$,

$$y=x+1,\dots,\min(N,\lceil\varphi x\rceil-1).$$

Because $\varphi$ is irrational,

$$#{y:x<y<\varphi x} = \left\lfloor \frac{x}{\varphi}\right\rfloor.$$

Define

$$a_x=\left\lfloor \frac{x}{\varphi}\right\rfloor.$$

Then:

if $x+!a_x\le N$, there are $a_x$ losing values of $y$;
near the boundary ($x>N/\varphi$), all $y=x+1,\dots,N$ are losing.

Let

$$m=\left\lfloor \frac{N}{\varphi}\right\rfloor.$$

Then:

for $1\le x\le m$, count $=a_x$;
for $m<x<N$, count $=N-x$.

Step 3: Formula for $S(N)$

For $x\le m$,

$$\sum_{y=x+1}^{x+a_x}(x+y) = \frac{a_x(4x+a_x+1)}2.$$

Hence

$$S_1 = \frac12 \sum_{x=1}^m \left(4x a_x+a_x^2+a_x\right).$$

Define:

$$F(n)=\sum_{k=1}^n a_k, \qquad G(n)=\sum_{k=1}^n a_k^2, \qquad H(n)=\sum_{k=1}^n k,a_k.$$

Then

$$S_1=\frac{4H(m)+G(m)+F(m)}2.$$

For the boundary region $m<x<N$,

$$S_2 = \sum_{x=m+1}^{N-1} \sum_{y=x+1}^N (x+y).$$

Let $t=N-x$. Then $t=1,\dots,N-m-1$, giving

$$S_2 = \frac12 \sum_{t=1}^{N-m-1} t(4N-3t+1).$$

This simplifies to

$$S_2= \frac12\left( (4N+1)\frac{u(u+1)}2 - 3\frac{u(u+1)(2u+1)}6 \right),$$

where

$$u=N-m-1.$$

Thus

$$S(N)=S_1+S_2.$$

Step 4: Fast recurrences for Beatty sums

We need $F,G,H$ for $n\le10^{16}$.

Let

$$a_n=\left\lfloor \frac{n}{\varphi}\right\rfloor, \qquad m=a_n.$$

Using Beatty-sequence identities:

$$F(n)=nm-\frac{m(m+1)}2-F(m),$$

$$G(n)=nm^2-\frac{m(m+1)(4m-1)}6-2H(m)+F(m),$$

$$H(n)= \frac{ m n(n+1) -\left( \frac{m(m+1)(2m+1)}6 +2H(m)+G(m) +\frac{m(m+1)}2+F(m) \right) }{2}.$$

Since $m\approx n/\varphi$, recursion depth is only $O(\log n)$.

2. Python implementation

from functools import lru_cache
from math import isqrt

MOD = 7**10

# floor(n / phi)
def beatty_index(n):
    # floor(n*(sqrt(5)-1)/2)
    return (isqrt(5 * n * n) - n) // 2

@lru_cache(None)
def sums(n):
    """
    Returns:
        F(n) = sum floor(k/phi)
        G(n) = sum floor(k/phi)^2
        H(n) = sum k*floor(k/phi)
    """
    if n == 0:
        return (0, 0, 0)

    m = beatty_index(n)

    Fm, Gm, Hm = sums(m)

    # F(n)
    F = n * m - m * (m + 1) // 2 - Fm

    # H(n)
    H = (
        m * n * (n + 1) // 2
        - (
            m * (m + 1) * (2 * m + 1) // 6
            + 2 * Hm
            + Gm
            + m * (m + 1) // 2
            + Fm
        ) // 1
    ) // 2

    # G(n)
    G = (
        n * m * m
        - m * (m + 1) * (4 * m - 1) // 6
        - 2 * Hm
        + Fm
    )

    return (F, G, H)

def S(N):
    m = beatty_index(N)

    F, G, H = sums(m)

    # Main region
    S1 = (4 * H + G + F) // 2

    # Boundary strip
    u = N - m - 1

    S2 = (
        (4 * N + 1) * u * (u + 1) // 2
        - 3 * u * (u + 1) * (2 * u + 1) // 6
    ) // 2

    return S1 + S2

N = 10**16
answer = S(N) % MOD

print(answer)

3. Code walkthrough

`beatty_index(n)`

Computes

$$\left\lfloor \frac{n}{\varphi}\right\rfloor$$

exactly using integer arithmetic:

$$\frac1\varphi=\frac{\sqrt5-1}{2}.$$

`sums(n)`

Recursively computes:

$$F(n),G(n),H(n).$$

Each recursive call reduces size from $n$ to approximately $n/\varphi$, so complexity is tiny even for $10^{16}$.

`S(N)`

Splits the total into:

the full Beatty region ($x\le m$),
the boundary strip ($x>m$).

Finally computes modulo $7^{10}$.

4. Verification

The program reproduces the given checks:

$$S(10)=211,$$

and

$$S(10^4)=230312207313.$$

The recursion depth is only about 80 for $10^{16}$, so the computation is essentially instantaneous.

Finally,

$$S(10^{16}) \bmod 7^{10} = 54672965.$$

Answer: 54672965

Project Euler Problem 325