Project Euler

Stone Game II

A game is played with two piles of stones. On each turn, a player removes k (smaller pile) stones from the larger pile, where k >= 1. The player who empties a pile wins. A position (a, b) with a <...

Source sync Apr 19, 2026

Problem #0325

Level Level 18

Solved By 720

Languages C++, Python

Answer 54672965

Length 464 words

modular_arithmeticgame_theorynumber_theory

A game is played with two piles of stones and two players.

On each player's turn, the player may remove a number of stones from the larger pile.

The number of stones removed must be a positive multiple of the number of stones in the smaller pile.

E.g. Let the ordered pair $(6,14)$ describe a configuration with $6$ stones in the smaller pile and $14$ stones in the larger pile, then the first player can remove $6$ or $12$ stones from the larger pile.

The player taking all the stones from a pile wins the game.

A winning configuration is one where the first player can force a win. For example, $(1,5)$, $(2,6)$, and $(3,12)$ are winning configurations because the first player can immediately remove all stones in the second pile.

A losing configuration is one where the second player can force a win, no matter what the first player does. For example, $(2,3)$ and $(3,4)$ are losing configurations: any legal move leaves a winning configuration for the second player.

Define $S(N)$ as the sum of $(x_i + y_i)$ for all losing configurations $(x_i, y_i), 0 < x_i < y_i \le N$.

We can verify that $S(10) = 211$ and $S(10^4) = 230312207313$.

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

Problem 325: Stone Game II

Mathematical Foundation

Theorem 1 (Losing positions have quotient 1). All losing positions $(a, b)$ with $a < b$ satisfy $a < b < 2a$ (i.e., $\lfloor b/a \rfloor = 1$ ).

Proof. Let $q = \lfloor b/a \rfloor$ and $r = b \bmod a$ . From position $(a, b)$ with $q \geq 2$ , the current player can move to $(r, a)$ (by choosing $k = q$ ) or to $(a, a + r)$ (by choosing $k = q - 1$ ). Now $(a, a + r)$ is losing iff $(r, a)$ is winning (since $(a, a+r)$ with quotient $\geq 2$ would recurse, but after one move from $(a, a+r)$ we reach $(r, a)$ ). Hence exactly one of $(r, a)$ and $(a, a+r)$ is losing, and the current player can always move to the losing one. Therefore $(a, b)$ is winning whenever $q \geq 2$ . $\square$

Theorem 2 (Beatty threshold). Position $(a, a + r)$ with $0 < r < a$ is losing if and only if $r \leq L(a)$ , where

L(a) = \left\lfloor \frac{a}{\varphi} \right\rfloor = \left\lfloor \frac{a(\sqrt{5} - 1)}{2} \right\rfloor.

Proof. Since all losing positions have quotient 1, the only move from $(a, a+r)$ is to $(r, a)$ , so $(a, a+r)$ is losing iff $(r, a)$ is winning. Define the set of losing ratios $\mathcal{L} = \{r/a : (a, a+r) \text{ is losing}\}$ . The game tree mirrors the Euclidean algorithm, and the parity of the first index at which the continued fraction $[q_1; q_2, \ldots]$ of $b/a$ has $q_i \geq 2$ determines the winner. Specifically, $(a, a+r)$ is losing iff the continued fraction of $a/r$ has its first entry $\geq 2$ at an odd index.

The losing ratios form the union of intervals:

\mathcal{L} = \bigcup_{k=0}^{\infty} \left(\frac{F_{2k}}{F_{2k+1}}, \frac{F_{2k+2}}{F_{2k+3}}\right]

where $F_i$ are the Fibonacci numbers. These intervals are contiguous and their union converges to $(0, 1/\varphi)$ where $\varphi = (1+\sqrt{5})/2$ . Therefore $L(a) = \lfloor a/\varphi \rfloor = \lfloor a(\sqrt{5}-1)/2 \rfloor$ . $\square$

Lemma 1 (Sum decomposition). We have

S(N) = \sum_{a=1}^{N-1} \sum_{r=1}^{R(a)} (2a + r)

where $R(a) = \min(L(a), N - a)$ , and this splits at $a^* = \max\{a : L(a) \leq N - a\} \approx N/\varphi$ into:

Range 1 ( $a \leq a^*$ ): $R(a) = L(a)$ , contributing Beatty-type sums.
Range 2 ( $a > a^*$ ): $R(a) = N - a$ , contributing a closed-form polynomial sum.

Proof. Each losing position is $(a, a+r)$ with $1 \leq r \leq L(a)$ and $a + r \leq N$ . The split at $a^*$ separates the regime where $L(a)$ is the binding constraint from the regime where $N - a$ is binding. $\square$

Theorem 3 (Quadratic floor-sum algorithm). The sums

f_0 = \sum_{x=1}^{A} \lfloor x\alpha \rfloor, \quad f_1 = \sum_{x=1}^{A} x\lfloor x\alpha \rfloor, \quad f_2 = \sum_{x=1}^{A} \lfloor x\alpha \rfloor^2

with $\alpha = (\sqrt{5}-1)/2$ approximated by a sufficiently good rational $p/q$ (e.g., $F_{79}/F_{80}$ ), can be computed simultaneously in $O(\log A)$ arithmetic operations via a Euclidean reciprocal reduction.

Proof. The standard floor-sum identity $\sum_{x=0}^{N-1} \lfloor (ax+b)/m \rfloor = NM - 1 + \text{reciprocal sum}$ where $M = \lfloor (a(N-1)+b)/m \rfloor$ reduces $(a, m)$ in the manner of the Euclidean algorithm. Extending this to track the quadratic and cross terms $f_1, f_2$ yields a system of six simultaneous recurrences that undergo the same Euclidean reduction. Since $\alpha \approx F_{79}/F_{80}$ (a ratio of consecutive Fibonacci numbers), the algorithm terminates in at most 80 steps. $\square$

Editorial

S(N) = sum of (a+b) for all losing positions (a,b) with 0 < a < b <= N. A position (a, a+r) is losing iff r <= floor(a*(sqrt(5)-1)/2). Uses quadratic floor sum algorithm with Fibonacci rational approximation. We find a* = max a such that floor(a * alpha) <= N - a. We then range 1: a = 1..a_star, R(a) = floor(a * alpha). Finally, range 2: a = a_star+1..N-1, R(a) = N - a.

Pseudocode

Find a* = max a such that floor(a * alpha) <= N - a
Range 1: a = 1..a_star, R(a) = floor(a * alpha)
Range 2: a = a_star+1..N-1, R(a) = N - a

Complexity Analysis

Time: $O(\log N)$ for the quadratic floor-sum algorithm. The Euclidean recursion depth for $F_{79}/F_{80}$ is at most 80 steps, each involving $O(1)$ arithmetic operations modulo $7^{10}$ .
Space: $O(\log N)$ for the recursion stack.

Answer

\boxed{54672965}

C++ project_euler/problem_325/solution.cpp

#include <bits/stdc++.h>
using namespace std;

typedef long long ll;
typedef __int128 lll;

const ll MOD = 282475249LL; // 7^10 (NOT prime!)

ll md(lll x) {
    ll r = (ll)(x % (lll)MOD);
    return r < 0 ? r + MOD : r;
}

// Extended GCD for modular inverse (MOD is not prime, can't use Fermat)
ll modinv(ll a, ll m) {
    ll g = m, x = 0, y = 1;
    ll aa = a % m;
    if (aa < 0) aa += m;
    while (aa != 0) {
        ll q = g / aa;
        g -= q * aa; swap(g, aa);
        x -= q * y; swap(x, y);
    }
    return (x % m + m) % m;
}

const ll inv2 = modinv(2, MOD);
const ll inv6 = modinv(6, MOD);

ll mul(ll a, ll b) { return md((lll)a * b); }
ll add(ll a, ll b) { return (a + b) % MOD; }
ll sub(ll a, ll b) { return (a - b % MOD + MOD) % MOD; }

ll S1n(lll n) { return mul(mul(md(n), md(n - 1)), inv2); }
ll S2n(lll n) { return mul(mul(mul(md(n), md(n - 1)), md(2 * n - 1)), inv6); }

struct T3 { ll f0, f1, f2; };

T3 fsq(lll N, lll a, lll b, lll m) {
    if (N <= 0) return {0, 0, 0};

    if (a >= m) {
        lll qa = a / m;
        auto g = fsq(N, a % m, b, m);
        ll s1 = S1n(N), s2 = S2n(N);
        ll qm = md(qa);
        ll f0 = add(mul(qm, s1), g.f0);
        ll f1 = add(mul(qm, s2), g.f1);
        ll f2 = add(add(mul(mul(qm, qm), s2), mul(mul(2, qm), g.f1)), g.f2);
        return {f0, f1, f2};
    }

    if (b >= m) {
        lll qb = b / m;
        auto g = fsq(N, a, b % m, m);
        ll qm = md(qb), nm = md(N);
        ll f0 = add(mul(qm, nm), g.f0);
        ll f1 = add(mul(qm, S1n(N)), g.f1);
        ll f2 = add(add(mul(mul(qm, qm), nm), mul(mul(2, qm), g.f0)), g.f2);
        return {f0, f1, f2};
    }

    if (a == 0) return {0, 0, 0};

    lll Mv = (a * (N - 1) + b) / m;
    if (Mv == 0) return {0, 0, 0};

    auto g = fsq(Mv, m, m - b - 1, a);

    ll nm = md(N - 1), mm = md(Mv), T = S1n(N);

    ll f0 = sub(mul(nm, mm), g.f0);
    ll f1 = sub(mul(mm, T), mul(inv2, add(g.f2, g.f0)));
    ll f2 = sub(sub(mul(mul(nm, mm), mm), mul(2, g.f1)), g.f0);

    return {f0, f1, f2};
}

ll isqrt_safe(lll n) {
    if (n <= 0) return 0;
    lll s = (lll)sqrtl((long double)n);
    while (s > 0 && s * s > n) s--;
    while ((s + 1) * (s + 1) <= n) s++;
    return (ll)s;
}

ll Lf(ll a) {
    lll n = (lll)a * a * 5;
    ll s = isqrt_safe(n);
    return (s - a) / 2;
}

int main() {
    ll N = 10000000000000000LL; // 10^16
    lll P = 14472334024676221LL; // F_79
    lll Q = 23416728348467685LL; // F_80

    // Binary search for a_star
    ll lo = 1, hi = N - 1;
    while (lo < hi) {
        ll mid = lo + (hi - lo + 1) / 2;
        if (Lf(mid) <= N - mid)
            lo = mid;
        else
            hi = mid - 1;
    }
    ll a_star = lo;

    // Range 1: a = 1..a_star, R(a) = L(a)
    auto r = fsq((lll)a_star + 1, P, 0, Q);
    ll range1 = add(mul(2, r.f1), mul(inv2, add(r.f2, r.f0)));

    // Range 2: a = a_star+1..N-1, R(a) = N-a
    ll J = N - a_star - 1;
    ll range2 = 0;
    if (J > 0) {
        ll jm = md(J), j1 = md(J+1), j2 = md(2*J+1), nm = md(N);
        ll sj = mul(mul(jm, j1), inv2);
        ll sj2 = mul(mul(mul(jm, j1), j2), inv6);
        range2 = sub(mul(add(mul(2, nm), inv2), sj), mul(mul(3, inv2), sj2));
    }

    cout << add(range1, range2) << endl;
    return 0;
}

Python project_euler/problem_325/solution.py

"""
Problem 325: Stone Game II

S(N) = sum of (a+b) for all losing positions (a,b) with 0 < a < b <= N.
A position (a, a+r) is losing iff r <= floor(a*(sqrt(5)-1)/2).

Uses quadratic floor sum algorithm with Fibonacci rational approximation.
"""
from math import isqrt

MOD = 282475249  # 7^10
inv2 = pow(2, -1, MOD)
inv6 = pow(6, -1, MOD)

# Large Fibonacci numbers for rational approximation of (sqrt(5)-1)/2
fib = [0, 1]
for i in range(100):
    fib.append(fib[-1] + fib[-2])
P, Q = fib[79], fib[80]  # alpha ~ P/Q, exact for a <= 10^16

def S1(n):
    """sum_{x=0}^{n-1} x mod MOD"""
    return n % MOD * ((n - 1) % MOD) % MOD * inv2 % MOD

def S2(n):
    """sum_{x=0}^{n-1} x^2 mod MOD"""
    return n % MOD * ((n - 1) % MOD) % MOD * ((2 * n - 1) % MOD) % MOD * inv6 % MOD

def floor_sum_quad(N, a, b, m):
    """
    Compute mod MOD:
    f0 = sum_{x=0}^{N-1} floor((ax+b)/m)
    f1 = sum_{x=0}^{N-1} x * floor((ax+b)/m)
    f2 = sum_{x=0}^{N-1} floor((ax+b)/m)^2
    """
    if N <= 0:
        return (0, 0, 0)

    if a >= m:
        qa = a // m
        g0, g1, g2 = floor_sum_quad(N, a % m, b, m)
        s1, s2 = S1(N), S2(N)
        qa_m = qa % MOD
        f0 = (qa_m * s1 + g0) % MOD
        f1 = (qa_m * s2 + g1) % MOD
        f2 = (qa_m * qa_m % MOD * s2 + 2 * qa_m * g1 + g2) % MOD
        return (f0 % MOD, f1 % MOD, f2 % MOD)

    if b >= m:
        qb = b // m
        g0, g1, g2 = floor_sum_quad(N, a, b % m, m)
        qb_m, Nm = qb % MOD, N % MOD
        f0 = (qb_m * Nm + g0) % MOD
        f1 = (qb_m * S1(N) + g1) % MOD
        f2 = (qb_m * qb_m % MOD * Nm + 2 * qb_m * g0 + g2) % MOD
        return (f0 % MOD, f1 % MOD, f2 % MOD)

    if a == 0:
        return (0, 0, 0)

    M_val = (a * (N - 1) + b) // m
    if M_val == 0:
        return (0, 0, 0)

    g0, g1, g2 = floor_sum_quad(M_val, m, m - b - 1, a)

    Nm = (N - 1) % MOD
    Mm = M_val % MOD
    T = S1(N)

    f0 = (Nm * Mm - g0) % MOD
    f1 = (Mm * T - inv2 * (g2 + g0)) % MOD
    f2 = (Nm * Mm % MOD * Mm - 2 * g1 - g0) % MOD

    return (f0 % MOD, f1 % MOD, f2 % MOD)

def L(a):
    """L(a) = floor(a * (sqrt(5)-1)/2)"""
    s = isqrt(5 * a * a)
    while (s + 1) * (s + 1) <= 5 * a * a:
        s += 1
    return (s - a) // 2

def compute_S(N):
    """Compute S(N) mod MOD."""
    # Binary search for threshold a*: largest a with L(a) <= N - a
    lo, hi = 1, N - 1
    while lo < hi:
        mid = (lo + hi + 1) // 2
        if L(mid) <= N - mid:
            lo = mid
        else:
            hi = mid - 1
    a_star = lo

    # Range 1: a = 1..a*, R(a) = L(a)
    # Contribution = 2*f1 + (f2 + f0)/2
    f0_1, f1_1, f2_1 = floor_sum_quad(a_star + 1, P, 0, Q)
    range1 = (2 * f1_1 + inv2 * (f2_1 + f0_1)) % MOD

    # Range 2: a = a*+1..N-1, R(a) = N-a
    # Let j = N - a, j = 1..J where J = N - a* - 1
    # sum_{j=1}^{J} [(2N + 1/2)*j - 3/2*j^2]
    J = N - a_star - 1
    if J > 0:
        Jm = J % MOD
        J1m = (J + 1) % MOD
        J2m = (2 * J + 1) % MOD
        Nm = N % MOD
        sum_j = Jm * J1m % MOD * inv2 % MOD
        sum_j2 = Jm * J1m % MOD * J2m % MOD * inv6 % MOD
        range2 = ((2 * Nm + inv2) % MOD * sum_j - 3 * inv2 % MOD * sum_j2) % MOD
    else:
        range2 = 0

    return (range1 + range2) % MOD

def solve():
    N = 10**16
    result = compute_S(N)
    print(result)

if __name__ == "__main__":
    solve()