Problem 984: Knights and Horses

Project Euler Problem 984 Solution

EulerSolve provides an optimized solution for Project Euler Problem 984, Knights and Horses, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Problem 984, Knights and Horses , ultimately asks for the value of a specific integer sequence at the enormous index $n=10^{18}$, reduced modulo $10^9+7$. The implementations make that sequence concrete through its first terms $1,4,9,92,903,4411,14959,41083,98200,212418,425756,803074,1441065,2479669$. The key point is that the computation is not done by any direct combinatorial search at size $10^{18}$. Instead, the sequence has already been compressed into two exact algebraic descriptions: a parity-sensitive closed form for the tail of the sequence, and a linear-recurrence evaluator that can reproduce the remaining cases modulo the prime modulus. Mathematical Approach The three implementations expose the same structure: after the first few exceptional indices, the Knights and Horses sequence is a quasi-polynomial in $n$ with period 2. The recurrence engine and the explicit even formula are two faces of the same underlying mathematics. The Sequence Encoded by the Implementations From the stored initial values and the signed recurrence used in the exact checkpoint, the tail of the sequence satisfies $$ f_n=7f_{n-1}-19f_{n-2}+21f_{n-3}+6f_{n-4}-42f_{n-5}+42f_{n-6}-6f_{n-7}-21f_{n-8}+19f_{n-9}-7f_{n-10}+f_{n-11}, \qquad n\ge 15. $$ This is the actual algebraic law extracted from the implementations....

Detailed mathematical approach

Problem Summary

Problem 984, Knights and Horses, ultimately asks for the value of a specific integer sequence at the enormous index $n=10^{18}$, reduced modulo $10^9+7$. The implementations make that sequence concrete through its first terms $1,4,9,92,903,4411,14959,41083,98200,212418,425756,803074,1441065,2479669$.

The key point is that the computation is not done by any direct combinatorial search at size $10^{18}$. Instead, the sequence has already been compressed into two exact algebraic descriptions: a parity-sensitive closed form for the tail of the sequence, and a linear-recurrence evaluator that can reproduce the remaining cases modulo the prime modulus.

Mathematical Approach

The three implementations expose the same structure: after the first few exceptional indices, the Knights and Horses sequence is a quasi-polynomial in $n$ with period 2. The recurrence engine and the explicit even formula are two faces of the same underlying mathematics.

The Sequence Encoded by the Implementations

From the stored initial values and the signed recurrence used in the exact checkpoint, the tail of the sequence satisfies

$$ f_n=7f_{n-1}-19f_{n-2}+21f_{n-3}+6f_{n-4}-42f_{n-5}+42f_{n-6}-6f_{n-7}-21f_{n-8}+19f_{n-9}-7f_{n-10}+f_{n-11}, \qquad n\ge 15. $$

This is the actual algebraic law extracted from the implementations. The final Project Euler answer is simply the value of this sequence at a huge even index.

The Transition Polynomial and Its Factorization

The generic evaluator stores fourteen basis values, so its transition polynomial is

$$ \chi(x)=x^{14}-7x^{13}+19x^{12}-21x^{11}-6x^{10}+42x^9-42x^8+6x^7+21x^6-19x^5+7x^4-x^3. $$

Factoring it gives

$$ \chi(x)=x^3(x-1)^9(x+1)^2. $$

The factor $x^3$ explains why the smallest indices are treated separately. After that transient prefix, the behavior is governed by the roots $1$ and $-1$: multiplicity 9 at $1$, and multiplicity 2 at $-1$.

Why the Tail Becomes a Parity-Dependent Polynomial

A root $1$ of multiplicity 9 contributes a polynomial of degree at most 8, while a root $-1$ of multiplicity 2 contributes $(-1)^n$ times a polynomial of degree at most 1. Therefore, for $n\ge 4$, the sequence must have the form

$$ f(n)=A_8(n)+(-1)^nB_1(n), $$

with $\deg A_8\le 8$ and $\deg B_1\le 1$. Matching this shape against the values encoded in the implementations yields

$$ f(n)=\frac{31}{40320}n^8+\frac{31}{3360}n^7+\frac{67}{1440}n^6+\frac{41}{320}n^5+\frac{313}{1440}n^4-\frac{5699}{240}n^3+\frac{16049}{420}n^2+\frac{941251}{4480}n-\frac{107261}{256}-(-1)^n\frac{2n+3}{256}, \qquad n\ge 4. $$

This single formula explains why the code can switch between a polynomial shortcut and a generic recurrence solver without changing the underlying sequence.

Recovering the Even Closed Form

The required query is $n=10^{18}$, which is even. Setting $(-1)^n=1$ collapses the quasi-polynomial to the degree-8 polynomial used directly by the implementations:

$$ P(n)=\frac{31}{40320}n^8+\frac{31}{3360}n^7+\frac{67}{1440}n^6+\frac{41}{320}n^5+\frac{313}{1440}n^4-\frac{5699}{240}n^3+\frac{16049}{420}n^2+\frac{29413}{140}n-419, \qquad n>3,\ n\equiv 0\pmod 2. $$

Because the modulus $10^9+7$ is prime and larger than every denominator that appears here, each division is implemented as multiplication by a modular inverse.

Worked Example: Evaluating $f(16)$

The parity reduction is already enough to show the method on a concrete value. Since $16$ is even and greater than 3, we use the closed polynomial:

$$ f(16)=P(16)=6623481. $$

The recurrence branch produces the same value from earlier terms, so the polynomial and the recurrence are not independent tricks; they are two exact descriptions of the same tail sequence.

How the Code Works

Closed-Form Evaluation for the Target Index

The C++, Python, and Java implementations first notice that the requested index $10^{18}$ is even. They compute $n,n^2,\dots,n^8$ modulo $10^9+7$, replace each rational denominator by its modular inverse, and evaluate the polynomial $P(n)$ term by term. That is enough to obtain the final answer for the problem instance.

Generic Linear-Recurrence Evaluator

The implementations also keep a general recurrence solver. It represents $x^{n-1}$ modulo the transition polynomial $\chi(x)$, using repeated squaring on coefficient vectors. After each coefficient-vector multiplication, higher powers are reduced with the recurrence coefficients, which is the standard quotient-ring or Kitamasa viewpoint for linear recurrences. Although fourteen basis slots are stored, the last three recurrence coefficients are zero, so the genuine tail relation is the 11-term recurrence written above.

Validation Checkpoints

One implementation performs explicit sanity checks before the final evaluation. It verifies small values such as $f(3)=9$ and $f(5)=903$, checks the exact integer value $f(100)=8658918531876$ using signed 128-bit recurrence arithmetic, and confirms the modular checkpoint $f(10000)\equiv 377956308\pmod{10^9+7}$. These checks protect against indexing mistakes, sign errors, and incorrect modular reductions.

Complexity Analysis

For the actual Project Euler query, the runtime is $O(1)$: the target index is even, so only a fixed number of modular multiplications and inverses are needed to evaluate the degree-8 polynomial. Memory usage is also $O(1)$.

The retained generic evaluator is more general. With a fixed state size $K=14$, repeated squaring of reduced coefficient vectors costs $O(K^2\log n)$ time and $O(K)$ memory. Since only 11 recurrence coefficients are nonzero, the constant factor is small.

Footnotes and References

Project Euler problem page: https://projecteuler.net/problem=984
Linear recurrences with constant coefficients: Wikipedia - Linear recurrence with constant coefficients
Characteristic polynomial: Wikipedia - Characteristic polynomial
Finite differences and polynomial sequences: Wikipedia - Finite difference
Modular multiplicative inverse: Wikipedia - Modular multiplicative inverse
Exponentiation by squaring: Wikipedia - Exponentiation by squaring

Problem 984 source code

C++

#include <array>
#include <cstdint>
#include <iostream>
#include <vector>

using i64 = std::int64_t;
using u64 = std::uint64_t;
using i128 = __int128_t;

static constexpr i64 MOD = 1000000007LL;
static constexpr int K = 14;

static constexpr std::array<i64, K> INIT = {
    1LL,
    4LL,
    9LL,
    92LL,
    903LL,
    4411LL,
    14959LL,
    41083LL,
    98200LL,
    212418LL,
    425756LL,
    803074LL,
    1441065LL,
    2479669LL
};

static constexpr std::array<i64, K> REC = {
    7LL,
    MOD - 19LL,
    21LL,
    6LL,
    MOD - 42LL,
    42LL,
    MOD - 6LL,
    MOD - 21LL,
    19LL,
    MOD - 7LL,
    1LL,
    0LL,
    0LL,
    0LL
};

static constexpr std::array<i64, 11> REC_SIGNED = {
    7LL,
    -19LL,
    21LL,
    6LL,
    -42LL,
    42LL,
    -6LL,
    -21LL,
    19LL,
    -7LL,
    1LL
};

static i64 mod_pow(i64 base, u64 exp) {
    i64 result = 1;
    base %= MOD;
    while (exp > 0ULL) {
        if (exp & 1ULL) {
            result = static_cast<i128>(result) * base % MOD;
        }
        base = static_cast<i128>(base) * base % MOD;
        exp >>= 1ULL;
    }
    return result;
}

static i64 mod_inv(i64 x) {
    return mod_pow((x % MOD + MOD) % MOD, MOD - 2);
}

static i64 f_even_closed_form_mod(const u64 n) {
    const i64 x = static_cast<i64>(n % static_cast<u64>(MOD));
    std::array<i64, 9> p{};
    p[0] = 1;
    for (int i = 1; i <= 8; ++i) {
        p[i] = static_cast<i128>(p[i - 1]) * x % MOD;
    }

    const i64 inv40320 = mod_inv(40320);
    const i64 inv3360 = mod_inv(3360);
    const i64 inv1440 = mod_inv(1440);
    const i64 inv320 = mod_inv(320);
    const i64 inv240 = mod_inv(240);
    const i64 inv420 = mod_inv(420);
    const i64 inv140 = mod_inv(140);

    i64 ans = 0;
    ans = (ans + static_cast<i128>(31) * p[8] % MOD * inv40320) % MOD;
    ans = (ans + static_cast<i128>(31) * p[7] % MOD * inv3360) % MOD;
    ans = (ans + static_cast<i128>(67) * p[6] % MOD * inv1440) % MOD;
    ans = (ans + static_cast<i128>(41) * p[5] % MOD * inv320) % MOD;
    ans = (ans + static_cast<i128>(313) * p[4] % MOD * inv1440) % MOD;
    ans = (ans + static_cast<i128>(MOD - 5699) * p[3] % MOD * inv240) % MOD;
    ans = (ans + static_cast<i128>(16049) * p[2] % MOD * inv420) % MOD;
    ans = (ans + static_cast<i128>(29413) * p[1] % MOD * inv140) % MOD;
    ans = (ans + MOD - 419) % MOD;
    return ans;
}

static std::array<i64, K> combine(const std::array<i64, K>& a,
                                   const std::array<i64, K>& b) {
    std::array<i64, 2 * K> prod{};

    for (int i = 0; i < K; ++i) {
        if (a[i] == 0) continue;
        for (int j = 0; j < K; ++j) {
            if (b[j] == 0) continue;
            prod[i + j] = (prod[i + j] + static_cast<i128>(a[i]) * b[j]) % MOD;
        }
    }

    for (int i = 2 * K - 2; i >= K; --i) {
        if (prod[i] == 0) continue;
        for (int j = 1; j <= K; ++j) {
            prod[i - j] = (prod[i - j] + static_cast<i128>(prod[i]) * REC[j - 1]) % MOD;
        }
    }

    std::array<i64, K> out{};
    for (int i = 0; i < K; ++i) out[i] = prod[i];
    return out;
}

static i64 f_mod(const u64 n) {
    if (n == 0ULL) {
        return 0;
    }
    if (n > 3ULL && (n & 1ULL) == 0ULL) {
        return f_even_closed_form_mod(n);
    }
    if (n <= static_cast<u64>(K)) {
        return INIT[static_cast<std::size_t>(n - 1)] % MOD;
    }

    u64 exp = n - 1ULL;
    std::array<i64, K> poly{};
    std::array<i64, K> step{};
    poly[0] = 1;
    step[1] = 1;

    while (exp > 0ULL) {
        if (exp & 1ULL) {
            poly = combine(poly, step);
        }
        step = combine(step, step);
        exp >>= 1ULL;
    }

    i64 ans = 0;
    for (int i = 0; i < K; ++i) {
        ans = (ans + static_cast<i128>(poly[i]) * INIT[i]) % MOD;
    }
    return ans;
}

static i64 f_100_exact() {
    std::vector<i128> f(101, 0);
    for (int n = 1; n <= K; ++n) {
        f[n] = INIT[static_cast<std::size_t>(n - 1)];
    }

    // This recurrence is valid from n >= 15.
    for (int n = 15; n <= 100; ++n) {
        i128 cur = 0;
        for (int j = 1; j <= 11; ++j) {
            cur += static_cast<i128>(REC_SIGNED[static_cast<std::size_t>(j - 1)]) * f[n - j];
        }
        f[n] = cur;
    }

    return static_cast<i64>(f[100]);
}

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

    if (f_mod(3ULL) != 9LL ||
        f_mod(5ULL) != 903LL ||
        f_100_exact() != 8658918531876LL ||
        f_mod(10000ULL) != 377956308LL) {
        std::cerr << "Validation failed\n";
        return 1;
    }

    const u64 target = 1000000000000000000ULL;
    std::cout << f_mod(target) << '\n';
    return 0;
}

Python

MOD = 1000000007
K = 14

INIT = [
    1, 4, 9, 92, 903, 4411, 14959, 41083, 98200,
    212418, 425756, 803074, 1441065, 2479669
]

REC = [
    7, MOD - 19, 21, 6, MOD - 42, 42, MOD - 6,
    MOD - 21, 19, MOD - 7, 1, 0, 0, 0
]

def combine(a, b):
    prod = [0] * (2 * K)
    for i in range(K):
        if a[i] == 0: continue
        for j in range(K):
            if b[j] == 0: continue
            prod[i + j] = (prod[i + j] + a[i] * b[j]) % MOD
            
    for i in range(2 * K - 2, K - 1, -1):
        if prod[i] == 0: continue
        for j in range(1, K + 1):
            prod[i - j] = (prod[i - j] + prod[i] * REC[j - 1]) % MOD
            
    return prod[:K]

def f_even_closed_form_mod(n):
    x = n % MOD
    p = [0] * 9
    p[0] = 1
    for i in range(1, 9):
        p[i] = (p[i - 1] * x) % MOD
        
    def mod_inv(val):
        return pow(val % MOD, MOD - 2, MOD)
        
    inv40320 = mod_inv(40320)
    inv3360 = mod_inv(3360)
    inv1440 = mod_inv(1440)
    inv320 = mod_inv(320)
    inv240 = mod_inv(240)
    inv420 = mod_inv(420)
    inv140 = mod_inv(140)
    
    ans = 0
    ans = (ans + 31 * p[8] % MOD * inv40320) % MOD
    ans = (ans + 31 * p[7] % MOD * inv3360) % MOD
    ans = (ans + 67 * p[6] % MOD * inv1440) % MOD
    ans = (ans + 41 * p[5] % MOD * inv320) % MOD
    ans = (ans + 313 * p[4] % MOD * inv1440) % MOD
    ans = (ans + (MOD - 5699) * p[3] % MOD * inv240) % MOD
    ans = (ans + 16049 * p[2] % MOD * inv420) % MOD
    ans = (ans + 29413 * p[1] % MOD * inv140) % MOD
    ans = (ans + MOD - 419) % MOD
    return ans

def f_mod(n):
    if n == 0:
        return 0
    if n > 3 and (n & 1) == 0:
        return f_even_closed_form_mod(n)
    if n <= K:
        return INIT[n - 1] % MOD
        
    exp = n - 1
    poly = [0] * K
    step = [0] * K
    poly[0] = 1
    step[1] = 1
    
    while exp > 0:
        if exp & 1:
            poly = combine(poly, step)
        step = combine(step, step)
        exp >>= 1
        
    ans = 0
    for i in range(K):
        ans = (ans + poly[i] * INIT[i]) % MOD
    return ans

def solve():
    target = 1000000000000000000
    return str(f_mod(target))

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

Java

import java.math.BigInteger;

public class Euler984 {
    static final long MOD = 1000000007L;
    static final int K = 14;

    static final long[] INIT = {
            1L, 4L, 9L, 92L, 903L, 4411L, 14959L, 41083L, 98200L,
            212418L, 425756L, 803074L, 1441065L, 2479669L
    };

    static final long[] REC = {
            7L, MOD - 19L, 21L, 6L, MOD - 42L, 42L, MOD - 6L,
            MOD - 21L, 19L, MOD - 7L, 1L, 0L, 0L, 0L
    };

    static long modPow(long base, long exp) {
        long result = 1;
        base %= MOD;
        while (exp > 0) {
            if ((exp & 1) != 0) {
                result = BigInteger.valueOf(result).multiply(BigInteger.valueOf(base)).mod(BigInteger.valueOf(MOD))
                        .longValue();
            }
            base = BigInteger.valueOf(base).multiply(BigInteger.valueOf(base)).mod(BigInteger.valueOf(MOD)).longValue();
            exp >>= 1;
        }
        return result;
    }

    static long modInv(long x) {
        return modPow((x % MOD + MOD) % MOD, MOD - 2);
    }

    static long fEvenClosedFormMod(long n) {
        long x = n % MOD;
        long[] p = new long[9];
        p[0] = 1;
        for (int i = 1; i <= 8; ++i) {
            p[i] = BigInteger.valueOf(p[i - 1]).multiply(BigInteger.valueOf(x)).mod(BigInteger.valueOf(MOD))
                    .longValue();
        }

        long inv40320 = modInv(40320);
        long inv3360 = modInv(3360);
        long inv1440 = modInv(1440);
        long inv320 = modInv(320);
        long inv240 = modInv(240);
        long inv420 = modInv(420);
        long inv140 = modInv(140);

        long ans = 0;
        ans = (ans + BigInteger.valueOf(31).multiply(BigInteger.valueOf(p[8])).mod(BigInteger.valueOf(MOD))
                .multiply(BigInteger.valueOf(inv40320)).mod(BigInteger.valueOf(MOD)).longValue()) % MOD;
        ans = (ans + BigInteger.valueOf(31).multiply(BigInteger.valueOf(p[7])).mod(BigInteger.valueOf(MOD))
                .multiply(BigInteger.valueOf(inv3360)).mod(BigInteger.valueOf(MOD)).longValue()) % MOD;
        ans = (ans + BigInteger.valueOf(67).multiply(BigInteger.valueOf(p[6])).mod(BigInteger.valueOf(MOD))
                .multiply(BigInteger.valueOf(inv1440)).mod(BigInteger.valueOf(MOD)).longValue()) % MOD;
        ans = (ans + BigInteger.valueOf(41).multiply(BigInteger.valueOf(p[5])).mod(BigInteger.valueOf(MOD))
                .multiply(BigInteger.valueOf(inv320)).mod(BigInteger.valueOf(MOD)).longValue()) % MOD;
        ans = (ans + BigInteger.valueOf(313).multiply(BigInteger.valueOf(p[4])).mod(BigInteger.valueOf(MOD))
                .multiply(BigInteger.valueOf(inv1440)).mod(BigInteger.valueOf(MOD)).longValue()) % MOD;
        ans = (ans + BigInteger.valueOf(MOD - 5699).multiply(BigInteger.valueOf(p[3])).mod(BigInteger.valueOf(MOD))
                .multiply(BigInteger.valueOf(inv240)).mod(BigInteger.valueOf(MOD)).longValue()) % MOD;
        ans = (ans + BigInteger.valueOf(16049).multiply(BigInteger.valueOf(p[2])).mod(BigInteger.valueOf(MOD))
                .multiply(BigInteger.valueOf(inv420)).mod(BigInteger.valueOf(MOD)).longValue()) % MOD;
        ans = (ans + BigInteger.valueOf(29413).multiply(BigInteger.valueOf(p[1])).mod(BigInteger.valueOf(MOD))
                .multiply(BigInteger.valueOf(inv140)).mod(BigInteger.valueOf(MOD)).longValue()) % MOD;
        ans = (ans + MOD - 419) % MOD;

        return ans;
    }

    static long[] combine(long[] a, long[] b) {
        long[] prod = new long[2 * K];
        for (int i = 0; i < K; ++i) {
            if (a[i] == 0)
                continue;
            for (int j = 0; j < K; ++j) {
                if (b[j] == 0)
                    continue;
                prod[i + j] = (prod[i + j] + BigInteger.valueOf(a[i]).multiply(BigInteger.valueOf(b[j]))
                        .mod(BigInteger.valueOf(MOD)).longValue()) % MOD;
            }
        }

        for (int i = 2 * K - 2; i >= K; --i) {
            if (prod[i] == 0)
                continue;
            for (int j = 1; j <= K; ++j) {
                prod[i - j] = (prod[i - j] + BigInteger.valueOf(prod[i]).multiply(BigInteger.valueOf(REC[j - 1]))
                        .mod(BigInteger.valueOf(MOD)).longValue()) % MOD;
            }
        }

        long[] out = new long[K];
        System.arraycopy(prod, 0, out, 0, K);
        return out;
    }

    static long fMod(long n) {
        if (n == 0)
            return 0;
        if (n > 3 && (n & 1) == 0)
            return fEvenClosedFormMod(n);
        if (n <= K)
            return INIT[(int) (n - 1)] % MOD;

        long exp = n - 1;
        long[] poly = new long[K];
        long[] step = new long[K];
        poly[0] = 1;
        step[1] = 1;

        while (exp > 0) {
            if ((exp & 1) != 0) {
                poly = combine(poly, step);
            }
            step = combine(step, step);
            exp >>= 1;
        }

        long ans = 0;
        for (int i = 0; i < K; ++i) {
            ans = (ans + BigInteger.valueOf(poly[i]).multiply(BigInteger.valueOf(INIT[i])).mod(BigInteger.valueOf(MOD))
                    .longValue()) % MOD;
        }
        return ans;
    }

    public static String solve() {
        long target = 1000000000000000000L;
        return String.valueOf(fMod(target));
    }

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