Problem 258: A Lagged Fibonacci Sequence

Mathematical Foundation

Theorem 1 (Characteristic Polynomial). The linear recurrence $g(k) = g(k-2000) + g(k-1999)$ has characteristic polynomial

Proof. Substituting the ansatz $g(k) = x^k$ into the recurrence:

Theorem 2 (Reduction to Polynomial Arithmetic). Let $M$ be the $2000 \times 2000$ companion matrix of the recurrence. By the Cayley—Hamilton theorem, $p(M) = 0$. Therefore, computing $g(10^{18})$ reduces to computing $x^{10^{18}} \bmod p(x)$ in the polynomial ring $\mathbb{Z}_m[x]$, where $m = 20092010$.

Proof. The state vector $\mathbf{v}_k = (g(k), g(k-1), \ldots, g(k-1999))^\top$ satisfies $\mathbf{v}_k = M \mathbf{v}_{k-1}$, so $\mathbf{v}_K = M^K \mathbf{v}_0$. By Cayley—Hamilton, $p(M) = 0$, so any polynomial in $M$ can be reduced modulo $p$. In particular, $M^K$ is determined by $x^K \bmod p(x)$. $\square$

Lemma 1 (Sparse Reduction). The reduction $x^{2000} \equiv x + 1 \pmod{p(x)}$ is extremely efficient: for a polynomial of degree $d \ge 2000$, each term $c \cdot x^d$ is replaced by $c \cdot x^{d-1999} + c \cdot x^{d-2000}$, requiring only 2 additions per excess degree.

Proof. From $x^{2000} = x + 1$ in the quotient ring, we get $x^d = x^{d-2000} \cdot x^{2000} = x^{d-2000}(x+1) = x^{d-1999} + x^{d-2000}$. Applying this rule iteratively from the highest-degree term downward reduces any polynomial to degree $< 2000$. Each application involves exactly 2 coefficient additions. $\square$

Theorem 3 (Extraction of $g(K)$). If $x^K \equiv \sum_{i=0}^{1999} a_i x^i \pmod{p(x)}$ in $\mathbb{Z}_m[x]$, then

since all initial values are $g(i) = 1$ for $0 \le i \le 1999$.

Proof. The general term of the recurrence is a linear combination of initial values: $g(K) = \sum_{i=0}^{1999} c_i g(i)$ where the coefficients $c_i$ are determined by $M^K$. The polynomial $x^K \bmod p(x)$ encodes exactly these coefficients (this follows from the isomorphism between the polynomial quotient ring and the matrix algebra generated by $M$). With $g(i) = 1$ for all $i$, the sum simplifies. $\square$

Editorial

g(k) = g(k-2000) + g(k-1999) for k >= 2000 g(k) = 1 for 0 <= k <= 1999 Find g(10^18) mod 20092010. Characteristic polynomial: p(x) = x^2000 - x - 1 So x^2000 ≡ x + 1 (mod p(x)) Compute x^(10^18) mod p(x) in Z_m[x], then g(10^18) = sum of coefficients (since all g(i) = 1 for i < 2000). We compute x^K mod p(x) in Z_m[x], where p(x) = x^2000 - x - 1. We then multiply polynomials A, B (each degree < 2000). Finally, result C has degree < 4000.

Pseudocode

K = 10^18, m = 20092010
Compute x^K mod p(x) in Z_m[x], where p(x) = x^2000 - x - 1
Multiply polynomials A, B (each degree < 2000)
Result C has degree < 4000
Reduce C modulo p(x): for d = deg(C) down to 2000:
C[d-1999] += C[d]
C[d-2000] += C[d]
C[d] = 0
All operations mod m
Binary exponentiation
if K is odd
Extract answer: sum of coefficients

Complexity Analysis

• Time: Binary exponentiation performs $O(\log K)$ polynomial multiplications. Each multiplication of two degree-$< 2000$ polynomials costs $O(n^2)$ (or $O(n \log n)$ with FFT) for the convolution, plus $O(n)$ for sparse reduction. Total: $O(n^2 \log K)$. With $n = 2000$ and $K = 10^{18}$: approximately $2000^2 \times 60 = 2.4 \times 10^8$ operations.
• Space: $O(n)$ for storing the polynomial coefficients ($n = 2000$).

Answer

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

/*
 * Problem 258: A Lagged Fibonacci Sequence
 *
 * g(k) = g(k-2000) + g(k-1999) for k >= 2000
 * g(k) = 1 for 0 <= k <= 1999
 *
 * Find g(10^18) mod 20092010.
 *
 * Characteristic polynomial: x^2000 - x - 1 = 0
 * So x^2000 ≡ x + 1 (mod p(x))
 *
 * We compute x^(10^18) mod p(x) in Z_m[x], then:
 * g(10^18) = sum of coefficients a_i * g(i) = sum of a_i (since g(i)=1 for i<2000)
 *
 * Answer: 12747994
 */

const int N = 2000;
const long long MOD = 20092010;

typedef vector<long long> Poly;

// Reduce polynomial modulo p(x) = x^2000 - x - 1
// i.e., x^2000 ≡ x + 1
void reduce(Poly& a) {
    while ((int)a.size() > N) {
        long long coeff = a.back();
        int deg = a.size() - 1;
        a.pop_back();
        // x^deg = coeff * x^(deg-2000) * x^2000 = coeff * x^(deg-2000) * (x + 1)
        // = coeff * (x^(deg-1999) + x^(deg-2000))
        int base = deg - N;
        a[base] = (a[base] + coeff) % MOD;
        a[base + 1] = (a[base + 1] + coeff) % MOD;
    }
}

// Multiply two polynomials mod p(x) mod MOD
Poly polymul(const Poly& a, const Poly& b) {
    int sa = a.size(), sb = b.size();
    Poly c(sa + sb - 1, 0);
    for (int i = 0; i < sa; i++) {
        if (a[i] == 0) continue;
        for (int j = 0; j < sb; j++) {
            c[i + j] = (c[i + j] + a[i] * b[j]) % MOD;
        }
    }
    reduce(c);
    return c;
}

int main() {
    // Compute x^(10^18) mod p(x) mod 20092010
    long long exp = 1000000000000000000LL;

    // result = 1 (the polynomial "1")
    Poly result(1, 1);
    // base = x
    Poly base(2, 0);
    base[1] = 1;

    while (exp > 0) {
        if (exp & 1) {
            result = polymul(result, base);
        }
        base = polymul(base, base);
        exp >>= 1;
    }

    // g(10^18) = sum of a_i * g(i) where g(i) = 1 for 0 <= i <= 1999
    long long ans = 0;
    for (int i = 0; i < (int)result.size(); i++) {
        ans = (ans + result[i]) % MOD;
    }

    cout << ans << endl;
    return 0;
}

"""
Problem 258: A Lagged Fibonacci Sequence

g(k) = g(k-2000) + g(k-1999) for k >= 2000
g(k) = 1 for 0 <= k <= 1999

Find g(10^18) mod 20092010.

Characteristic polynomial: p(x) = x^2000 - x - 1
So x^2000 ≡ x + 1 (mod p(x))

Compute x^(10^18) mod p(x) in Z_m[x], then
g(10^18) = sum of coefficients (since all g(i) = 1 for i < 2000).

Answer: 12747994
"""

N = 2000
MOD = 20092010

def reduce_poly(a):
    """Reduce polynomial modulo p(x) = x^2000 - x - 1.
    Uses: x^2000 ≡ x + 1 (mod p(x))."""
    while len(a) > N:
        coeff = a.pop()
        deg = len(a)  # after pop, this is the original degree
        base = deg - N
        a[base] = (a[base] + coeff) % MOD
        a[base + 1] = (a[base + 1] + coeff) % MOD
    return a

def poly_mul(a, b):
    """Multiply two polynomials mod p(x) mod MOD."""
    sa, sb = len(a), len(b)
    c = [0] * (sa + sb - 1)
    for i in range(sa):
        if a[i] == 0:
            continue
        ai = a[i]
        for j in range(sb):
            c[i + j] = (c[i + j] + ai * b[j]) % MOD
    return reduce_poly(c)

def poly_pow(exp_val):
    """Compute x^exp_val mod p(x) mod MOD using binary exponentiation."""
    result = [1]
    base = [0, 1]

    while exp_val > 0:
        if exp_val & 1:
            result = poly_mul(result, base)
        base = poly_mul(base, base)
        exp_val >>= 1

    return result

def main():
    exp_val = 10**18
    result = poly_pow(exp_val)

    # g(10^18) = sum of coefficients a_i * g(i), and g(i) = 1 for all i < 2000
    ans = sum(result) % MOD
    print(ans)

if __name__ == "__main__":
    main()