Problem 895: Gold & Silver Coin Game II

Project Euler Problem 895 Solution

EulerSolve provides an optimized solution for Project Euler Problem 895, Gold & Silver Coin Game II, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Problem 895 asks for the value of $G(9898)$ in the Gold & Silver Coin Game II, reduced modulo $M=989898989$. The C++, Python, and Java implementations do not simulate game positions directly. Instead they evaluate an exact closed decomposition of $G(n)$ consisting of a parity-sensitive kernel, a recurrence-defined coefficient sequence, and a short correction term. This reformulation is what makes the computation practical. Every ingredient can be generated in a single forward pass, and the only divisions are modular inverses, so the whole calculation stays inside $\mathbb{Z}/M\mathbb{Z}$. Mathematical Approach Let $u$ denote the modular inverse of $3$, so $3u\equiv 1\pmod M$. For $n\ge 2$, the implementations compute $G(n)$ from three pieces: a kernel sequence $K(m)$, a coefficient sequence $c_i$, and a parity-dependent prefactor $B(n)$. Step 1: Isolate the kernel sequence The inner term depends only on the argument $m$ and on its parity. Written in mathematical form, the implemented kernel is $$K(m)=\begin{cases} 0, & m\le 0,\\ 1, & m=1,\\ 14\cdot 2^{m/2}u+6m-15 \pmod M, & m\ge 2,\ m\text{ even},\\ 20\cdot 2^{(m-1)/2}u+6m-15 \pmod M, & m\ge 3,\ m\text{ odd}. \end{cases}$$ The power of $2$ grows only with $\lfloor m/2\rfloor$, while the coefficient in front of it switches from $14$ to $20$ when parity changes....

Detailed mathematical approach

Problem Summary

Problem 895 asks for the value of $G(9898)$ in the Gold & Silver Coin Game II, reduced modulo $M=989898989$. The C++, Python, and Java implementations do not simulate game positions directly. Instead they evaluate an exact closed decomposition of $G(n)$ consisting of a parity-sensitive kernel, a recurrence-defined coefficient sequence, and a short correction term.

This reformulation is what makes the computation practical. Every ingredient can be generated in a single forward pass, and the only divisions are modular inverses, so the whole calculation stays inside $\mathbb{Z}/M\mathbb{Z}$.

Mathematical Approach

Let $u$ denote the modular inverse of $3$, so $3u\equiv 1\pmod M$. For $n\ge 2$, the implementations compute $G(n)$ from three pieces: a kernel sequence $K(m)$, a coefficient sequence $c_i$, and a parity-dependent prefactor $B(n)$.

Step 1: Isolate the kernel sequence

The inner term depends only on the argument $m$ and on its parity. Written in mathematical form, the implemented kernel is

$$K(m)=\begin{cases} 0, & m\le 0,\\ 1, & m=1,\\ 14\cdot 2^{m/2}u+6m-15 \pmod M, & m\ge 2,\ m\text{ even},\\ 20\cdot 2^{(m-1)/2}u+6m-15 \pmod M, & m\ge 3,\ m\text{ odd}. \end{cases}$$

The power of $2$ grows only with $\lfloor m/2\rfloor$, while the coefficient in front of it switches from $14$ to $20$ when parity changes. That even/odd split is one of the main structural features of the final formula.

Step 2: Build the coefficient sequence

The convolution weights form a second-order recurrence:

$$c_0=1,\qquad c_1=10,$$

$$c_i\equiv \left((20i-10)c_{i-1}-(64i-64)c_{i-2}\right)i^{-1}\pmod M \qquad (i\ge 2).$$

The first few values are

$$1,\ 10,\ 118,\ 1540,\ 21286,\ \dots$$

Because $M$ is prime and every index $i$ used in the computation satisfies $1\le i<M$, the inverse $i^{-1}$ always exists modulo $M$.

Step 3: Package the recurrence with a generating function

If we define

$$A(x)=\sum_{i\ge 0} c_i x^i,$$

then the recurrence is equivalent to the differential equation

$$\left(1-20x+64x^2\right)A'(x)=\left(10-64x\right)A(x),\qquad A(0)=1.$$

Solving it gives

$$A(x)=\frac{1}{\sqrt{1-20x+64x^2}}=\frac{1}{\sqrt{(1-4x)(1-16x)}}.$$

This viewpoint explains why the recurrence weights naturally appear in a one-dimensional convolution: they are exactly the coefficients of a fixed algebraic generating function.

Step 4: Assemble the final formula for $G(n)$

The remaining term is a small parity-dependent correction:

$$B(n)\equiv \begin{cases} 35\cdot 2^{n/2}u-6n-35u \pmod M, & n\text{ even},\\ 50\cdot 2^{(n-1)/2}u-6n-35u \pmod M, & n\text{ odd}. \end{cases}$$

With that notation, the implemented identity is

$$\boxed{G(n)\equiv B(n)+\sum_{i=0}^{\lfloor (n-1)/2\rfloor} c_i\,K(n-2i)\pmod M.}$$

The shift by $2$ means the parity of $n$ is preserved throughout the sum. Even inputs only touch the even branch of the kernel, while odd inputs walk through the odd branch and finally terminate at $K(1)=1$.

Step 5: Why the modular divisions are legitimate

The formula contains divisions by $3$ and by the recurrence index $i$. In modular arithmetic these are replaced by multiplicative inverses. Since $M=989898989$ is prime, every nonzero residue has an inverse, so the computation is exact in the finite field $\mathbb{Z}/M\mathbb{Z}$; nothing is being approximated numerically.

Worked Example: $n=5$

The recurrence gives

$$c_2=\frac{(20\cdot 2-10)c_1-(64\cdot 2-64)c_0}{2}=\frac{30\cdot 10-64}{2}=118.$$

Because $n=5$ is odd, the prefactor is

$$B(5)=50\cdot 2^2u-30-35u=165u-30\equiv 25\pmod M.$$

The kernel values needed by the sum are

$$K(5)=20\cdot 2^2u+15=80u+15,\qquad K(3)=20\cdot 2u+3=40u+3,\qquad K(1)=1.$$

Therefore

$$\begin{aligned} G(5)&\equiv B(5)+c_0K(5)+c_1K(3)+c_2K(1) \\ &\equiv 25+(80u+15)+10(40u+3)+118 \\ &\equiv 188+480u \\ &\equiv 188+160 \\ &\equiv 348 \pmod M. \end{aligned}$$

This matches the known checkpoint used by the implementations.

How the Code Works

The C++, Python, and Java implementations all follow the same numerical plan. They first compute the inverse of $3$, evaluate the parity-dependent prefactor $B(n)$, and determine the limit $\lfloor (n-1)/2\rfloor$. They then build the coefficient list $c_0,c_1,\dots,c_{\lfloor (n-1)/2\rfloor}$ iteratively from the recurrence above.

After that, the implementation scans through the same index range once more, evaluates the kernel at $n-2i$, multiplies by the stored coefficient $c_i$, and accumulates the answer modulo $M$. The C++ version obtains inverses with the extended Euclidean algorithm, while the Python and Java versions use fast modular exponentiation under the prime modulus. In every language, the algorithmic structure is identical.

Complexity Analysis

Let $L=\lfloor (n-1)/2\rfloor$. Building the recurrence coefficients costs $O(L)$ arithmetic steps, and the final convolution sum costs another $O(L)$ steps. Because the modulus is fixed for the problem, this is linear time in $n$. If one counts modular inverse computation explicitly, the running time is $O(L\log M)$, which is still effectively linear here. The memory usage is $O(L)$ because the coefficient list is stored for the final pass.

Footnotes and References

Problem page: https://projecteuler.net/problem=895
Recurrence relation: Wikipedia — Recurrence relation
Generating function: Wikipedia — Generating function
Cauchy product: Wikipedia — Cauchy product
Modular multiplicative inverse: Wikipedia — Modular multiplicative inverse

Problem 895 source code

C++

#include <cstdint>
#include <iostream>
#include <stdexcept>
#include <vector>

using namespace std;

static constexpr long long MOD = 989898989LL;

static long long norm_mod(long long x) {
    x %= MOD;
    if (x < 0) x += MOD;
    return x;
}

static long long add_mod(long long a, long long b) {
    long long v = a + b;
    if (v >= MOD) v -= MOD;
    return v;
}

static long long sub_mod(long long a, long long b) {
    long long v = a - b;
    if (v < 0) v += MOD;
    return v;
}

static long long mul_mod(long long a, long long b) {
    return static_cast<long long>(
        (static_cast<unsigned long long>(a) * static_cast<unsigned long long>(b)) %
        static_cast<unsigned long long>(MOD));
}

static long long mod_pow(long long base, long long exp) {
    long long result = 1;
    long long cur = norm_mod(base);
    while (exp > 0) {
        if (exp & 1LL) result = mul_mod(result, cur);
        cur = mul_mod(cur, cur);
        exp >>= 1LL;
    }
    return result;
}

static long long mod_inv(long long a) {
    long long t = 0, new_t = 1;
    long long r = MOD, new_r = norm_mod(a);
    while (new_r != 0) {
        long long q = r / new_r;
        long long tmp_t = t - q * new_t;
        t = new_t;
        new_t = tmp_t;
        long long tmp_r = r - q * new_r;
        r = new_r;
        new_r = tmp_r;
    }
    if (r != 1) {
        throw runtime_error("mod inverse does not exist");
    }
    return norm_mod(t);
}

static long long g8(long long n, long long inv3) {
    if (n <= 0) return 0;
    if (n == 1) return 1;
    long long coeff = (n % 2 == 0) ? 14 : 20;
    long long term1 = mul_mod(coeff, mul_mod(mod_pow(2, n / 2), inv3));
    long long term2 = norm_mod(6 * n - 15);
    return add_mod(term1, term2);
}

static long long compute_G(int n) {
    if (n <= 1) return 0;

    long long inv3 = mod_inv(3);
    long long pow2 = mod_pow(2, n / 2);
    long long pref = (n % 2 == 0) ? 35 : 50;
    long long counter = mul_mod(pref, mul_mod(pow2, inv3));
    counter = sub_mod(counter, norm_mod(6LL * n));
    counter = sub_mod(counter, mul_mod(35, inv3));

    int lim = (n - 1) / 2;
    vector<long long> nums(lim + 1, 0);
    nums[0] = 1;
    if (lim >= 1) nums[1] = 10;

    for (int i = 2; i <= lim; ++i) {
        long long a = mul_mod(norm_mod(20LL * i - 10), nums[i - 1]);
        long long b = mul_mod(norm_mod(64LL * i - 64), nums[i - 2]);
        long long val = sub_mod(a, b);
        nums[i] = mul_mod(val, mod_inv(i));
    }

    for (int i = 0; i <= lim; ++i) {
        counter = add_mod(counter, mul_mod(nums[i], g8(n - 2LL * i, inv3)));
    }
    return counter;
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    struct Check {
        int n;
        long long expected;
    };
    const Check checks[] = {
        {2, 6},
        {5, 348},
        {20, 125825982708LL},
    };

    for (const auto &check : checks) {
        long long got = compute_G(check.n);
        long long exp = check.expected % MOD;
        if (got != exp) {
            cerr << "Validation failed for G(" << check.n << "): got " << got
                 << ", expected " << exp << ".\n";
            return 1;
        }
    }
    cerr << "Validation checkpoints passed.\n";

    cout << compute_G(9898) << "\n";
    return 0;
}

Python

MOD = 989898989

def norm_mod(x):
    return x % MOD

def add_mod(a, b):
    return (a + b) % MOD

def sub_mod(a, b):
    return (a - b) % MOD

def mul_mod(a, b):
    return (a * b) % MOD

def mod_pow(base, exp):
    return pow(base, exp, MOD)

def mod_inv(a):
    return pow(a, MOD - 2, MOD)

def g8(n, inv3):
    if n <= 0: return 0
    if n == 1: return 1
    coeff = 14 if n % 2 == 0 else 20
    term1 = mul_mod(coeff, mul_mod(mod_pow(2, n // 2), inv3))
    term2 = norm_mod(6 * n - 15)
    return add_mod(term1, term2)

def compute_G(n):
    if n <= 1: return 0

    inv3 = mod_inv(3)
    pow2 = mod_pow(2, n // 2)
    pref = 35 if n % 2 == 0 else 50
    counter = mul_mod(pref, mul_mod(pow2, inv3))
    counter = sub_mod(counter, norm_mod(6 * n))
    counter = sub_mod(counter, mul_mod(35, inv3))

    lim = (n - 1) // 2
    nums = [0] * (lim + 1)
    nums[0] = 1
    if lim >= 1:
        nums[1] = 10

    for i in range(2, lim + 1):
        a = mul_mod(norm_mod(20 * i - 10), nums[i - 1])
        b = mul_mod(norm_mod(64 * i - 64), nums[i - 2])
        val = sub_mod(a, b)
        nums[i] = mul_mod(val, mod_inv(i))

    for i in range(lim + 1):
        counter = add_mod(counter, mul_mod(nums[i], g8(n - 2 * i, inv3)))
        
    return counter % MOD

def solve():
    return str(compute_G(9898))

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

Java

public class Euler895 {
    static final long MOD = 989898989L;

    static long normMod(long x) {
        x %= MOD;
        if (x < 0)
            x += MOD;
        return x;
    }

    static long addMod(long a, long b) {
        long v = a + b;
        if (v >= MOD)
            v -= MOD;
        return v;
    }

    static long subMod(long a, long b) {
        long v = a - b;
        if (v < 0)
            v += MOD;
        return v;
    }

    static long mulMod(long a, long b) {
        return (a * b) % MOD;
    }

    static long modPow(long base, long exp) {
        long result = 1;
        long cur = normMod(base);
        while (exp > 0) {
            if ((exp & 1) == 1)
                result = mulMod(result, cur);
            cur = mulMod(cur, cur);
            exp >>= 1;
        }
        return result;
    }

    static long modInv(long a) {
        return modPow(a, MOD - 2);
    }

    static long g8(long n, long inv3) {
        if (n <= 0)
            return 0;
        if (n == 1)
            return 1;
        long coeff = (n % 2 == 0) ? 14 : 20;
        long term1 = mulMod(coeff, mulMod(modPow(2, n / 2), inv3));
        long term2 = normMod(6 * n - 15);
        return addMod(term1, term2);
    }

    static long computeG(int n) {
        if (n <= 1)
            return 0;

        long inv3 = modInv(3);
        long pow2 = modPow(2, n / 2);
        long pref = (n % 2 == 0) ? 35 : 50;
        long counter = mulMod(pref, mulMod(pow2, inv3));
        counter = subMod(counter, normMod(6L * n));
        counter = subMod(counter, mulMod(35, inv3));

        int lim = (n - 1) / 2;
        long[] nums = new long[lim + 1];
        nums[0] = 1;
        if (lim >= 1)
            nums[1] = 10;

        for (int i = 2; i <= lim; ++i) {
            long a = mulMod(normMod(20L * i - 10), nums[i - 1]);
            long b = mulMod(normMod(64L * i - 64), nums[i - 2]);
            long val = subMod(a, b);
            nums[i] = mulMod(val, modInv(i));
        }

        for (int i = 0; i <= lim; ++i) {
            counter = addMod(counter, mulMod(nums[i], g8(n - 2L * i, inv3)));
        }
        return counter;
    }

    public static String solve() {
        return Long.toString(computeG(9898));
    }

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