Problem 366: Stone Game III

View on Project Euler

Project Euler Problem 366 Solution

EulerSolve provides an optimized solution for Project Euler Problem 366, Stone Game III, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Let \(M(n)\) denote the largest winning first move in the Stone Game III position with \(n\) stones, and define $$S(N)=\sum_{n=1}^{N} M(n).$$ The Project Euler task asks for \(S(10^{18}) \bmod 10^8\). A direct scan of all positions up to \(10^{18}\) is impossible, so the local solution replaces per-position game search by a recurrence on Fibonacci intervals. Mathematical Approach 1. Winning-Move Criterion Used by the Verifier The C++ file contains a slow checker based on Zeckendorf representations. For a positive integer \(m\), let \(z(m)\) be the smallest Fibonacci term appearing in its Zeckendorf decomposition. The checker tests every move \(x\) with \(1 \le x \lt n\) and uses the characterization $$2x \lt z(n-x).$$ Therefore the game-theoretic quantity implemented in all three language versions is $$M(n)=\max\left\{x \in \{1,\dots,n-1\}: 2x \lt z(n-x)\right\}.$$ The optimized solver never evaluates this inequality for every \(n\). Instead, it exploits the block pattern that this rule creates when \(n\) is grouped by Fibonacci ranges. 2. Fibonacci Interval Decomposition Use the same Fibonacci indexing as the code: $$F_0=1,\qquad F_1=1,\qquad F_{k+1}=F_k+F_{k-1}.$$ For each \(k\), define the interval $$I_k=[F_k,\,F_{k+1}-1].$$ Its length is $$|I_k|=F_{k+1}-F_k=F_{k-1}.$$ Now let \(A_k\) be the sequence of values \(M(n)\) on \(I_k\)....

Detailed mathematical approach

Problem Summary

Let \(M(n)\) denote the largest winning first move in the Stone Game III position with \(n\) stones, and define

$$S(N)=\sum_{n=1}^{N} M(n).$$

The Project Euler task asks for \(S(10^{18}) \bmod 10^8\). A direct scan of all positions up to \(10^{18}\) is impossible, so the local solution replaces per-position game search by a recurrence on Fibonacci intervals.

Mathematical Approach

1. Winning-Move Criterion Used by the Verifier

The C++ file contains a slow checker based on Zeckendorf representations. For a positive integer \(m\), let \(z(m)\) be the smallest Fibonacci term appearing in its Zeckendorf decomposition. The checker tests every move \(x\) with \(1 \le x \lt n\) and uses the characterization

$$2x \lt z(n-x).$$

Therefore the game-theoretic quantity implemented in all three language versions is

$$M(n)=\max\left\{x \in \{1,\dots,n-1\}: 2x \lt z(n-x)\right\}.$$

The optimized solver never evaluates this inequality for every \(n\). Instead, it exploits the block pattern that this rule creates when \(n\) is grouped by Fibonacci ranges.

2. Fibonacci Interval Decomposition

Use the same Fibonacci indexing as the code:

$$F_0=1,\qquad F_1=1,\qquad F_{k+1}=F_k+F_{k-1}.$$

For each \(k\), define the interval

$$I_k=[F_k,\,F_{k+1}-1].$$

Its length is

$$|I_k|=F_{k+1}-F_k=F_{k-1}.$$

Now let \(A_k\) be the sequence of values \(M(n)\) on \(I_k\). Evaluating the slow definition on small intervals gives

$$A_5=(0,1,2,3,1),\qquad A_6=(0,1,2,3,4,5,6,2),$$

$$A_7=(0,1,\dots,10,3,1),\qquad A_8=(0,1,\dots,16,4,5,6,2).$$

This already shows the structure used by the fast solver: each Fibonacci interval starts with a long linear prefix, and the remaining tail is copied recursively from two intervals earlier.

3. Prefix Lengths and the Tail Recurrence

Let \(P_k\) be the length of the initial linear prefix of \(A_k\). Then the first \(P_k\) values in interval \(I_k\) are exactly

$$0,1,2,\dots,P_k-1.$$

After that prefix comes a tail sequence \(B_k\). The base cases visible in the code are

$$B_5=(1),\qquad B_6=(2).$$

For \(k \ge 7\), the tail begins with the consecutive block

$$P_{k-3},\,P_{k-3}+1,\dots,P_{k-2}-1,$$

and after that block it continues with the old tail \(B_{k-2}\).

Comparing lengths gives a recurrence for \(P_k\). Since the full interval has length \(F_{k-1}\), the tail length is \(F_{k-1}-P_k\). On the other hand, the recursive description says that this tail length is

$$\left(P_{k-2}-P_{k-3}\right)+\left(F_{k-3}-P_{k-2}\right)=F_{k-3}-P_{k-3}.$$

Hence

$$F_{k-1}-P_k=F_{k-3}-P_{k-3},$$

so

$$P_k=F_{k-2}+P_{k-3}\qquad (k \ge 6).$$

The implementation uses the initial values

$$P_1=1,\qquad P_2=1,\qquad P_3=2,\qquad P_4=3,\qquad P_5=4.$$

4. Summing a Full Fibonacci Interval

Let \(T_k\) be the sum of the tail \(B_k\). From the recursive tail structure we get

$$T_5=1,\qquad T_6=2,$$

and for \(k \ge 7\),

$$T_k=T_{k-2}+\sum_{i=P_{k-3}}^{P_{k-2}-1} i.$$

Therefore the total contribution of the whole interval \(I_k\) is

$$U_k=\sum_{n=F_k}^{F_{k+1}-1} M(n)=\sum_{i=0}^{P_k-1} i + T_k=\frac{(P_k-1)P_k}{2}+T_k.$$

This is exactly the quantity stored in interval_sum[k] in the C++ code and in the corresponding arrays in Python and Java.

5. Prefix Queries Inside the Last Interval

For a general limit \(N\), choose the unique \(k\) such that

$$F_k \le N \lt F_{k+1},$$

and write

$$j=N-F_k.$$

Then

$$S(N)=\sum_{r=1}^{k-1} U_r+\operatorname{Pref}_k(j),$$

where \(\operatorname{Pref}_k(j)\) is the sum of the first \(j+1\) entries of \(A_k\).

If \(j \lt P_k\), we are still inside the linear prefix, so

$$\operatorname{Pref}_k(j)=\sum_{i=0}^{j} i=\frac{j(j+1)}{2}.$$

If \(j \ge P_k\), we first take the whole linear prefix and then continue through the tail. Because the tail itself is built from a consecutive block plus the tail from interval \(k-2\), the prefix calculation descends by \(k \mapsto k-2\). That recursive evaluation is implemented by prefix_interval and prefix_tail.

6. Worked Checkpoint: \(S(100)=728\)

The slow verifier in the C++ program checks this value explicitly. The complete intervals \(I_1\) through \(I_9=[55,88]\) contribute

$$U_1+\cdots+U_9=662.$$

Now \(100\) lies in \(I_{10}=[89,143]\), so \(j=100-89=11\). Because \(P_{10}=45\), this index is still inside the linear prefix of interval 10. Therefore

$$\operatorname{Pref}_{10}(11)=0+1+\cdots+11=\frac{11 \cdot 12}{2}=66.$$

Thus

$$S(100)=662+66=728,$$

matching the checkpoint used by the repository implementation.

How the Code Works

The three solution files implement the same recurrence. They first generate Fibonacci numbers up to \(N\), then precompute the arrays prefix_len, tail_sum, interval_sum, and cumulative interval_prefix_sum. A query \(S(N)\) is answered by locating the last full interval below \(N\), adding the precomputed total before it, and evaluating the remaining partial interval recursively.

The language-specific differences are only numeric types: C++ uses unsigned __int128 for exact sums, Python relies on native big integers, and Java uses BigInteger. The C++ version also keeps the slow Zeckendorf-based verifier, checking \(S(100)=728\) and matching the optimized method against brute force for \(N=1000\).

Complexity Analysis

There are only \(O(\log N)\) Fibonacci levels up to \(N\). Building all tables therefore costs \(O(\log N)\) time and \(O(\log N)\) memory. Answering one query \(S(N)\) also takes \(O(\log N)\) time, because interval lookup and the recursive tail descent both cross only logarithmically many Fibonacci blocks.

Footnotes and References

  1. Problem page: https://projecteuler.net/problem=366
  2. Zeckendorf theorem: Wikipedia — Zeckendorf's theorem
  3. Fibonacci numbers: Wikipedia — Fibonacci number
  4. Arithmetic series: Wikipedia — Arithmetic progression

Problem 366 source code

C++

#include <algorithm>
#include <cstdint>
#include <iostream>
#include <string>
#include <vector>

namespace {

using u64 = std::uint64_t;
using u128 = unsigned __int128;

constexpr u64 kLimit = 1000000000000000000ULL;
constexpr u64 kMod = 100000000ULL;

struct Precomputed {
    std::vector<u64> fib;        // fib[0]=1, fib[1]=1, fib[2]=2, ...
    std::vector<u64> prefix_len; // P_k
    std::vector<u128> tail_sum;  // tail sum in interval k
    std::vector<u128> interval_sum;
    std::vector<u128> interval_prefix_sum;
};

u128 tri(const u64 n) {
    return static_cast<u128>(n) * static_cast<u128>(n + 1ULL) / 2U;
}

u128 sum_range(const u64 a, const u64 b) {
    if (b < a) {
        return 0;
    }
    const u64 len = b - a + 1ULL;
    return static_cast<u128>(a + b) * static_cast<u128>(len) / 2U;
}

Precomputed build_precomputed(const u64 max_n) {
    Precomputed data;
    data.fib = {1ULL, 1ULL, 2ULL};  // index 0..2
    while (data.fib.back() <= max_n) {
        const std::size_t sz = data.fib.size();
        data.fib.push_back(data.fib[sz - 1] + data.fib[sz - 2]);
    }

    const std::size_t kmax = data.fib.size() - 2;  // largest interval index with fib[k] <= max_n

    data.prefix_len.assign(data.fib.size(), 0ULL);
    data.tail_sum.assign(data.fib.size(), 0U);
    data.interval_sum.assign(data.fib.size(), 0U);
    data.interval_prefix_sum.assign(data.fib.size(), 0U);

    // Base prefix lengths.
    data.prefix_len[1] = 1;
    data.prefix_len[2] = 1;
    data.prefix_len[3] = 2;
    data.prefix_len[4] = 3;
    data.prefix_len[5] = 4;

    for (std::size_t k = 6; k <= kmax; ++k) {
        data.prefix_len[k] = data.fib[k - 2] + data.prefix_len[k - 3];
    }

    // Tail sums for intervals.
    data.tail_sum[5] = 1;
    data.tail_sum[6] = 2;
    for (std::size_t k = 7; k <= kmax; ++k) {
        const u64 a = data.prefix_len[k - 3];
        const u64 b = data.prefix_len[k - 2] - 1ULL;
        data.tail_sum[k] = data.tail_sum[k - 2] + sum_range(a, b);
    }

    // Sum of M(n) over a full interval [F_k, F_{k+1}-1].
    data.interval_sum[1] = 0;
    data.interval_sum[2] = 0;
    for (std::size_t k = 3; k <= kmax; ++k) {
        const u64 p = data.prefix_len[k];
        data.interval_sum[k] = tri(p - 1ULL) + data.tail_sum[k];
    }

    for (std::size_t k = 1; k <= kmax; ++k) {
        data.interval_prefix_sum[k] = data.interval_prefix_sum[k - 1] + data.interval_sum[k];
    }

    return data;
}

u128 prefix_tail(const Precomputed& data, const std::size_t k, const u64 j) {
    if (j == static_cast<u64>(-1)) {
        return 0;
    }
    if (k <= 4) {
        return 0;
    }
    if (k == 5) {
        return 1;
    }
    if (k == 6) {
        return 2;
    }

    const u64 start = data.prefix_len[k - 3];
    const u64 end = data.prefix_len[k - 2] - 1ULL;
    const u64 block_len = end - start + 1ULL;

    if (j < block_len) {
        return sum_range(start, start + j);
    }
    return sum_range(start, end) + prefix_tail(data, k - 2, j - block_len);
}

u128 prefix_interval(const Precomputed& data, const std::size_t k, const u64 idx) {
    const u64 interval_len = data.fib[k - 1];
    if (idx >= interval_len - 1ULL) {
        return data.interval_sum[k];
    }
    if (k <= 2) {
        return 0;
    }

    const u64 p = data.prefix_len[k];
    if (idx < p) {
        return tri(idx);
    }

    const u64 tail_idx = idx - p;
    return tri(p - 1ULL) + prefix_tail(data, k, tail_idx);
}

u128 fast_sum_M(const u64 n, const Precomputed& data) {
    if (n == 0ULL) {
        return 0;
    }

    std::size_t k = 1;
    while (k + 1 < data.fib.size() && data.fib[k + 1] <= n) {
        ++k;
    }

    const u128 full_before = data.interval_prefix_sum[k - 1];
    const u64 idx = n - data.fib[k];
    return full_before + prefix_interval(data, k, idx);
}

u64 smallest_zeckendorf_term(u64 n, const std::vector<u64>& fib12) {
    std::size_t i = fib12.size() - 1;
    u64 smallest = 0;
    while (n > 0ULL) {
        while (fib12[i] > n) {
            --i;
        }
        n -= fib12[i];
        smallest = fib12[i];
        if (i >= 2) {
            i -= 2;
        } else {
            break;
        }
    }
    return smallest;
}

u64 slow_sum_M(const int limit) {
    std::vector<u64> fib12 = {1ULL, 2ULL};
    while (fib12.back() <= static_cast<u64>(limit)) {
        const std::size_t sz = fib12.size();
        fib12.push_back(fib12[sz - 1] + fib12[sz - 2]);
    }

    std::vector<u64> z(static_cast<std::size_t>(limit + 1), 0ULL);
    for (int i = 1; i <= limit; ++i) {
        z[static_cast<std::size_t>(i)] = smallest_zeckendorf_term(static_cast<u64>(i), fib12);
    }

    u64 total = 0;
    for (int n = 1; n <= limit; ++n) {
        u64 best = 0;
        for (int x = 1; x < n; ++x) {
            if (2ULL * static_cast<u64>(x) < z[static_cast<std::size_t>(n - x)]) {
                best = static_cast<u64>(x);
            }
        }
        total += best;
    }
    return total;
}

bool run_checkpoints() {
    const Precomputed data = build_precomputed(1000000ULL);

    const u64 slow_100 = slow_sum_M(100);
    if (slow_100 != 728ULL) {
        std::cerr << "Checkpoint failed: slow sum for n<=100\n";
        return false;
    }

    const u128 fast_100 = fast_sum_M(100ULL, data);
    if (fast_100 != 728ULL) {
        std::cerr << "Checkpoint failed: fast sum for n<=100\n";
        return false;
    }

    const u64 slow_1000 = slow_sum_M(1000);
    const u128 fast_1000 = fast_sum_M(1000ULL, data);
    if (fast_1000 != slow_1000) {
        std::cerr << "Checkpoint failed: slow/fast mismatch for n<=1000\n";
        return false;
    }

    return true;
}

}  // namespace

int main(int argc, char** argv) {
    bool skip_checkpoints = false;
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            skip_checkpoints = true;
        } else {
            std::cerr << "Unknown argument: " << arg << '\n';
            return 1;
        }
    }

    if (!skip_checkpoints && !run_checkpoints()) {
        return 2;
    }

    const Precomputed data = build_precomputed(kLimit);
    const u128 total = fast_sum_M(kLimit, data);
    const u64 answer = static_cast<u64>(total % static_cast<u128>(kMod));
    std::cout << answer << '\n';
    return 0;
}

Python

def tri(n):
    return n * (n + 1) // 2

def sum_range(a, b):
    if b < a: return 0
    return (a + b) * (b - a + 1) // 2

class Precomputed:
    def __init__(self, max_n):
        self.fib = [1, 1, 2]
        while self.fib[-1] <= max_n:
            self.fib.append(self.fib[-1] + self.fib[-2])
            
        kmax = len(self.fib) - 2
        
        self.prefix_len = [0] * len(self.fib)
        self.tail_sum = [0] * len(self.fib)
        self.interval_sum = [0] * len(self.fib)
        self.interval_prefix_sum = [0] * len(self.fib)
        
        self.prefix_len[1] = 1
        self.prefix_len[2] = 1
        self.prefix_len[3] = 2
        self.prefix_len[4] = 3
        if len(self.prefix_len) > 5:
            self.prefix_len[5] = 4
            
        for k in range(6, kmax + 1):
            self.prefix_len[k] = self.fib[k - 2] + self.prefix_len[k - 3]
            
        if len(self.tail_sum) > 5: self.tail_sum[5] = 1
        if len(self.tail_sum) > 6: self.tail_sum[6] = 2
        for k in range(7, kmax + 1):
            a = self.prefix_len[k - 3]
            b = self.prefix_len[k - 2] - 1
            self.tail_sum[k] = self.tail_sum[k - 2] + sum_range(a, b)
            
        self.interval_sum[1] = 0
        self.interval_sum[2] = 0
        for k in range(3, kmax + 1):
            p = self.prefix_len[k]
            self.interval_sum[k] = tri(p - 1) + self.tail_sum[k]
            
        for k in range(1, kmax + 1):
            self.interval_prefix_sum[k] = self.interval_prefix_sum[k - 1] + self.interval_sum[k]

def prefix_tail(data, k, j):
    if j < 0: return 0
    if k <= 4: return 0
    if k == 5: return 1
    if k == 6: return 2
    
    start = data.prefix_len[k - 3]
    end = data.prefix_len[k - 2] - 1
    block_len = end - start + 1
    
    if j < block_len:
        return sum_range(start, start + j)
    return sum_range(start, end) + prefix_tail(data, k - 2, j - block_len)

def prefix_interval(data, k, idx):
    interval_len = data.fib[k - 1]
    if idx >= interval_len - 1:
        return data.interval_sum[k]
    if k <= 2: return 0
    
    p = data.prefix_len[k]
    if idx < p:
        return tri(idx)
        
    tail_idx = idx - p
    return tri(p - 1) + prefix_tail(data, k, tail_idx)

def fast_sum_M(n, data):
    if n == 0: return 0
    k = 1
    while k + 1 < len(data.fib) and data.fib[k + 1] <= n:
        k += 1
        
    full_before = data.interval_prefix_sum[k - 1]
    idx = n - data.fib[k]
    return full_before + prefix_interval(data, k, idx)

def solve():
    kLimit = 1000000000000000000
    kMod = 100000000
    data = Precomputed(kLimit)
    total = fast_sum_M(kLimit, data)
    ans = total % kMod
    return str(ans)

if __name__ == '__main__':
    print(solve())

Java

import java.math.BigInteger;
import java.util.ArrayList;
import java.util.List;

public class Euler366 {
    static final long kLimit = 1000000000000000000L;
    static final long kMod = 100000000L;

    static BigInteger tri(long n) {
        return BigInteger.valueOf(n).multiply(BigInteger.valueOf(n + 1)).divide(BigInteger.valueOf(2));
    }

    static BigInteger sumRange(long a, long b) {
        if (b < a)
            return BigInteger.ZERO;
        long len = b - a + 1;
        return BigInteger.valueOf(a + b).multiply(BigInteger.valueOf(len)).divide(BigInteger.valueOf(2));
    }

    static class Precomputed {
        List<Long> fib;
        long[] prefixLen;
        BigInteger[] tailSum;
        BigInteger[] intervalSum;
        BigInteger[] intervalPrefixSum;

        Precomputed(long maxN) {
            fib = new ArrayList<>();
            fib.add(1L);
            fib.add(1L);
            fib.add(2L);
            while (fib.get(fib.size() - 1) <= maxN) {
                int sz = fib.size();
                fib.add(fib.get(sz - 1) + fib.get(sz - 2));
            }

            int kmax = fib.size() - 2;
            int size = fib.size();

            prefixLen = new long[size];
            tailSum = new BigInteger[size];
            intervalSum = new BigInteger[size];
            intervalPrefixSum = new BigInteger[size];

            for (int i = 0; i < size; i++) {
                tailSum[i] = BigInteger.ZERO;
                intervalSum[i] = BigInteger.ZERO;
                intervalPrefixSum[i] = BigInteger.ZERO;
            }

            prefixLen[1] = 1;
            if (size > 2)
                prefixLen[2] = 1;
            if (size > 3)
                prefixLen[3] = 2;
            if (size > 4)
                prefixLen[4] = 3;
            if (size > 5)
                prefixLen[5] = 4;

            for (int k = 6; k <= kmax; k++) {
                prefixLen[k] = fib.get(k - 2) + prefixLen[k - 3];
            }

            if (size > 5)
                tailSum[5] = BigInteger.valueOf(1);
            if (size > 6)
                tailSum[6] = BigInteger.valueOf(2);

            for (int k = 7; k <= kmax; k++) {
                long a = prefixLen[k - 3];
                long b = prefixLen[k - 2] - 1;
                tailSum[k] = tailSum[k - 2].add(sumRange(a, b));
            }

            intervalSum[1] = BigInteger.ZERO;
            if (size > 2)
                intervalSum[2] = BigInteger.ZERO;

            for (int k = 3; k <= kmax; k++) {
                long p = prefixLen[k];
                intervalSum[k] = tri(p - 1).add(tailSum[k]);
            }

            for (int k = 1; k <= kmax; k++) {
                intervalPrefixSum[k] = intervalPrefixSum[k - 1].add(intervalSum[k]);
            }
        }
    }

    static BigInteger prefixTail(Precomputed data, int k, long j) {
        if (j < 0)
            return BigInteger.ZERO;
        if (k <= 4)
            return BigInteger.ZERO;
        if (k == 5)
            return BigInteger.ONE;
        if (k == 6)
            return BigInteger.valueOf(2);

        long start = data.prefixLen[k - 3];
        long end = data.prefixLen[k - 2] - 1;
        long blockLen = end - start + 1;

        if (j < blockLen) {
            return sumRange(start, start + j);
        }
        return sumRange(start, end).add(prefixTail(data, k - 2, j - blockLen));
    }

    static BigInteger prefixInterval(Precomputed data, int k, long idx) {
        long intervalLen = data.fib.get(k - 1);
        if (idx >= intervalLen - 1) {
            return data.intervalSum[k];
        }
        if (k <= 2)
            return BigInteger.ZERO;

        long p = data.prefixLen[k];
        if (idx < p) {
            return tri(idx);
        }

        long tailIdx = idx - p;
        return tri(p - 1).add(prefixTail(data, k, tailIdx));
    }

    static BigInteger fastSumM(long n, Precomputed data) {
        if (n == 0)
            return BigInteger.ZERO;

        int k = 1;
        while (k + 1 < data.fib.size() && data.fib.get(k + 1) <= n) {
            k++;
        }

        BigInteger fullBefore = data.intervalPrefixSum[k - 1];
        long idx = n - data.fib.get(k);
        return fullBefore.add(prefixInterval(data, k, idx));
    }

    static String solve() {
        Precomputed data = new Precomputed(kLimit);
        BigInteger total = fastSumM(kLimit, data);
        long ans = total.remainder(BigInteger.valueOf(kMod)).longValue();
        return Long.toString(ans);
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}