Problem 258: A Lagged Fibonacci Sequence

View on Project Euler

Project Euler Problem 258 Solution

EulerSolve provides an optimized solution for Project Euler Problem 258, A Lagged Fibonacci Sequence, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The sequence is defined by $$G_0=G_1=\cdots=G_{1999}=1,$$ and for every \(n\ge 2000\), $$G_n=G_{n-2000}+G_{n-1999}\pmod M,$$ where in the Project Euler problem $$M=20092010.$$ The challenge is to compute \(G_N\) for an enormous index such as $$N=10^{18}.$$ A direct simulation is impossible. Even storing all terms up to \(10^{18}\) is absurd, and even an \(O(N)\) recurrence loop is far too slow. The code therefore converts the recurrence into polynomial arithmetic and uses binary exponentiation. Mathematical Approach 1. View the Recurrence as a Shift Operator Introduce a formal shift symbol \(x\). In a linear recurrence, multiplying by \(x\) means “move one step forward in the sequence”. So \(x^k\) corresponds to shifting by \(k\) positions. The recurrence $$G_n=G_{n-2000}+G_{n-1999}$$ can be rewritten as $$G_{n+2000}=G_n+G_{n+1}.$$ If \(x\) is the shift operator, then “shift by 2000” is the same as “shift by 0 plus shift by 1”, so formally $$x^{2000}=1+x.$$ This is the key algebraic identity behind the whole solution. 2....

Detailed mathematical approach

Problem Summary

The sequence is defined by

$$G_0=G_1=\cdots=G_{1999}=1,$$

and for every \(n\ge 2000\),

$$G_n=G_{n-2000}+G_{n-1999}\pmod M,$$

where in the Project Euler problem

$$M=20092010.$$

The challenge is to compute \(G_N\) for an enormous index such as

$$N=10^{18}.$$

A direct simulation is impossible. Even storing all terms up to \(10^{18}\) is absurd, and even an \(O(N)\) recurrence loop is far too slow. The code therefore converts the recurrence into polynomial arithmetic and uses binary exponentiation.

Mathematical Approach

1. View the Recurrence as a Shift Operator

Introduce a formal shift symbol \(x\). In a linear recurrence, multiplying by \(x\) means “move one step forward in the sequence”. So \(x^k\) corresponds to shifting by \(k\) positions.

The recurrence

$$G_n=G_{n-2000}+G_{n-1999}$$

can be rewritten as

$$G_{n+2000}=G_n+G_{n+1}.$$

If \(x\) is the shift operator, then “shift by 2000” is the same as “shift by 0 plus shift by 1”, so formally

$$x^{2000}=1+x.$$

This is the key algebraic identity behind the whole solution.

2. The Characteristic Polynomial

Move everything to one side:

$$x^{2000}-x-1=0.$$

So the natural modulus polynomial is

$$P(x)=x^{2000}-x-1.$$

The solver works in the quotient ring

$$\mathbb Z_M[x]\big/\bigl(P(x)\bigr),$$

which means that inside computations we are always allowed to replace

$$x^{2000}\equiv x+1\pmod{P(x)}.$$

Any higher power \(x^n\) can therefore be reduced to a polynomial of degree at most 1999.

3. Why Reducing \(x^N\) Solves the Recurrence

Because the recurrence has order 2000, every term \(G_N\) is a linear combination of the first 2000 values.

Concretely, there exist coefficients \(c_0,\dots,c_{1999}\) such that

$$G_N\equiv \sum_{i=0}^{1999} c_i\,G_i \pmod M.$$

Those coefficients are exactly the remainder of \(x^N\) modulo \(P(x)\):

$$x^N \equiv \sum_{i=0}^{1999} c_i x^i \pmod{P(x)}.$$

This statement is the polynomial version of companion-matrix exponentiation. Instead of exponentiating a \(2000\times 2000\) matrix, the code exponentiates the indeterminate \(x\) in the quotient ring.

4. Why the Final Answer Is Just the Sum of Coefficients

Normally we would still need to evaluate

$$\sum_{i=0}^{1999} c_i\,G_i.$$

But here the initial block is especially simple:

$$G_0=G_1=\cdots=G_{1999}=1.$$

Therefore

$$G_N\equiv \sum_{i=0}^{1999} c_i \pmod M.$$

This is why the implementation can finish by simply summing the coefficients of the reduced polynomial.

5. A First Example: \(N=2000\)

From the reduction rule,

$$x^{2000}\equiv x+1.$$

So the coefficient vector has only two nonzero entries, both equal to 1.

Hence

$$G_{2000}\equiv 1+1=2 \pmod M.$$

This matches the recurrence directly:

$$G_{2000}=G_0+G_1=1+1=2.$$

6. A Second Example: \(N=2001\)

Multiply the previous identity by \(x\):

$$x^{2001}\equiv x(x+1)=x^2+x.$$

Therefore

$$G_{2001}\equiv G_1+G_2=1+1=2.$$

Again this agrees with the recurrence:

$$G_{2001}=G_1+G_2=2.$$

7. A Larger Worked Example: \(N=4000\)

Reduce

$$x^{4000}=(x^{2000})^2\equiv (x+1)^2=x^2+2x+1.$$

So the reduced polynomial is

$$x^{4000}\equiv 1+2x+x^2.$$

Hence

$$G_{4000}\equiv G_0+2G_1+G_2=1+2+1=4 \pmod M.$$

And indeed the direct recurrence produces

$$G_{3999}=3,\qquad G_{4000}=G_{2000}+G_{2001}=2+2=4.$$

This example shows very clearly what the coefficient vector means: it tells us how \(G_N\) depends on the first 2000 terms.

8. How Polynomial Multiplication and Reduction Work

Suppose we already have two reduced polynomials

$$A(x)=\sum_{i=0}^{1999} a_i x^i,\qquad B(x)=\sum_{j=0}^{1999} b_j x^j.$$

Their ordinary product has degree at most 3998:

$$A(x)B(x)=\sum_{t=0}^{3998} d_t x^t.$$

Whenever a term \(x^t\) has \(t\ge 2000\), the reduction rule says

$$x^t=x^{t-2000}x^{2000}\equiv x^{t-2000}(x+1)=x^{t-1999}+x^{t-2000}.$$

So each high-degree coefficient \(d_t\) is pushed into exactly two lower positions:

$$d_t x^t \longmapsto d_t x^{t-2000}+d_t x^{t-1999}.$$

This is exactly what the function multiply_mod does in the C++ code. It first computes the raw convolution, then scans degrees from high to low and folds each coefficient back with the rule above.

9. Why Binary Exponentiation Applies

We need \(x^N\) modulo \(P(x)\), and exponentiation by squaring works in any associative multiplication structure. The quotient ring is associative, so we can use the standard binary method:

1. Start with the multiplicative identity \(1\).

2. Let the current base be \(x\).

3. Read the binary expansion of \(N\).

4. When a bit is 1, multiply the accumulator by the current base.

5. Square the base at every step.

6. Reduce after every multiplication.

Thus the number of polynomial multiplications is only

$$O(\log N),$$

not \(O(N)\).

10. Relation to Matrix Exponentiation

A standard method for linear recurrences is to build the \(2000\times 2000\) companion matrix and raise it to the \(N\)-th power.

The present approach is mathematically equivalent, but much lighter. Instead of matrix-vector multiplication, we work with one degree-1999 polynomial. Since the recurrence only contains two nonzero coefficients, the reduction rule \(x^{2000}\equiv x+1\) is extremely simple, and the polynomial method becomes especially attractive.

11. Validation Checkpoints

The implementation contains several checks:

$$G_0=1,\qquad G_{1999}=1,$$

$$G_{2000}=2,\qquad G_{2001}=2.$$

Then it cross-checks the fast solver against a direct sliding-window recurrence at

$$N=25000.$$

That shared value is

$$G_{25000}=4096 \pmod{20092010}.$$

This validates both the quotient-ring interpretation and the binary-exponentiation implementation.

How the Code Works

The code fixes

$$K=2000$$

and represents every reduced polynomial as an array of length 2000. Entry \(i\) is the coefficient of \(x^i\).

The function multiply_mod first performs an ordinary quadratic convolution into a temporary array of length \(2K\). It then reduces every coefficient of degree at least \(K\) by sending it to positions \(i-K\) and \(i-K+1\), which is the array form of

$$x^K\equiv x+1.$$

The function solve initializes

$$\text{acc}=1,\qquad \text{base}=x,$$

raises \(x\) to the target power by binary exponentiation, and finally returns the sum of all coefficients in the resulting remainder polynomial.

Complexity Analysis

Let

$$K=2000.$$

One polynomial multiplication with reduction costs

$$O(K^2),$$

because the raw convolution is quadratic. Binary exponentiation uses

$$O(\log N)$$

such multiplications, so the total time is

$$O(K^2\log N).$$

The working memory is only

$$O(K),$$

since we store a few arrays of length about \(K\) or \(2K\), not a huge matrix.

Footnotes and References

  1. Problem page: https://projecteuler.net/problem=258
  2. Linear recurrences with constant coefficients: Wikipedia - Linear recurrence
  3. Exponentiation by squaring: Wikipedia - Exponentiation by squaring
  4. Polynomial rings and quotient rings: Wikipedia - Polynomial ring
  5. Companion matrix viewpoint: Wikipedia - Companion matrix

Problem 258 source code

C++

#include <array>
#include <cstdint>
#include <iostream>
#include <string>

namespace {

using u64 = std::uint64_t;

constexpr int K = 2000;

struct Options {
    u64 index = 1000000000000000000ULL;
    int mod = 20092010;
    bool run_checkpoints = true;
};

bool parse_u64_after_prefix(const std::string& arg, const std::string& prefix, u64& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }

    u64 parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<u64>(c - '0');
    }
    value = parsed;
    return true;
}

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }

    int parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<int>(c - '0');
    }
    value = parsed;
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_u64_after_prefix(arg, "--index=", options.index) ||
            parse_int_after_prefix(arg, "--mod=", options.mod)) {
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }

    return options.mod >= 2;
}

using Poly = std::array<int, K>;

Poly multiply_mod(const Poly& a, const Poly& b, const int mod) {
    std::array<long long, 2 * K> tmp{};

    for (int i = 0; i < K; ++i) {
        const long long ai = a[static_cast<std::size_t>(i)];
        if (ai == 0) {
            continue;
        }
        for (int j = 0; j < K; ++j) {
            const int bj = b[static_cast<std::size_t>(j)];
            if (bj == 0) {
                continue;
            }
            tmp[static_cast<std::size_t>(i + j)] += ai * bj;
            tmp[static_cast<std::size_t>(i + j)] %= mod;
        }
    }

    for (int i = 2 * K - 1; i >= K; --i) {
        const long long v = tmp[static_cast<std::size_t>(i)] % mod;
        if (v == 0) {
            continue;
        }
        tmp[static_cast<std::size_t>(i - K)] += v;
        tmp[static_cast<std::size_t>(i - K + 1)] += v;
        tmp[static_cast<std::size_t>(i - K)] %= mod;
        tmp[static_cast<std::size_t>(i - K + 1)] %= mod;
    }

    Poly out{};
    for (int i = 0; i < K; ++i) {
        out[static_cast<std::size_t>(i)] = static_cast<int>(tmp[static_cast<std::size_t>(i)] % mod);
    }
    return out;
}

int solve(const u64 index, const int mod) {
    if (index < K) {
        return 1 % mod;
    }

    Poly acc{};
    Poly base{};
    acc[0] = 1;
    base[1] = 1;  // x

    u64 exp = index;
    while (exp > 0) {
        if (exp & 1ULL) {
            acc = multiply_mod(acc, base, mod);
        }
        exp >>= 1ULL;
        if (exp != 0) {
            base = multiply_mod(base, base, mod);
        }
    }

    int answer = 0;
    for (int c : acc) {
        answer += c;
        if (answer >= mod) {
            answer -= mod;
        }
    }
    return answer;
}

int solve_bruteforce(const int index, const int mod) {
    if (index < K) {
        return 1 % mod;
    }
    std::array<int, K> window{};
    window.fill(1 % mod);
    for (int i = K; i <= index; ++i) {
        const int next = (window[0] + window[1]) % mod;
        for (int j = 0; j < K - 1; ++j) {
            window[static_cast<std::size_t>(j)] = window[static_cast<std::size_t>(j + 1)];
        }
        window[K - 1] = next;
    }
    return window[K - 1];
}

bool run_checkpoints() {
    if (solve(0, 20092010) != 1 || solve(1999, 20092010) != 1) {
        std::cerr << "Checkpoint failed for initial terms" << '\n';
        return false;
    }
    if (solve(2000, 20092010) != 2 || solve(2001, 20092010) != 2) {
        std::cerr << "Checkpoint failed for first recursive terms" << '\n';
        return false;
    }

    const int test_index = 25000;
    const int fast = solve(test_index, 20092010);
    const int brute = solve_bruteforce(test_index, 20092010);
    if (fast != brute) {
        std::cerr << "Checkpoint failed for cross-check index 25000" << '\n';
        return false;
    }

    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }

    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }

    std::cout << solve(options.index, options.mod) << '\n';
    return 0;
}

Python

def solve():
    K = 2000
    index = 1_000_000_000_000_000_000
    mod = 20092010

    def multiply_mod(a, b, m):
        tmp = [0] * (2 * K)
        for i in range(K):
            if a[i] == 0:
                continue
            for j in range(K):
                if b[j] == 0:
                    continue
                tmp[i + j] = (tmp[i + j] + a[i] * b[j]) % m
        for i in range(2 * K - 1, K - 1, -1):
            v = tmp[i] % m
            if v == 0:
                continue
            tmp[i - K] = (tmp[i - K] + v) % m
            tmp[i - K + 1] = (tmp[i - K + 1] + v) % m
        return [tmp[i] % m for i in range(K)]

    if index < K:
        return str(1 % mod)

    acc = [0] * K
    acc[0] = 1
    base = [0] * K
    base[1] = 1

    exp = index
    while exp > 0:
        if exp & 1:
            acc = multiply_mod(acc, base, mod)
        exp >>= 1
        if exp != 0:
            base = multiply_mod(base, base, mod)

    answer = sum(acc) % mod
    return str(answer)

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

Java

public class Euler258 {
    static final int K = 2000;
    static final int MOD = 20092010;

    static long[] multiplyMod(long[] a, long[] b) {
        long[] tmp = new long[2 * K];

        for (int i = 0; i < K; ++i) {
            long ai = a[i];
            if (ai == 0)
                continue;
            for (int j = 0; j < K; ++j) {
                long bj = b[j];
                if (bj == 0)
                    continue;
                tmp[i + j] = (tmp[i + j] + ai * bj) % MOD;
            }
        }

        for (int i = 2 * K - 1; i >= K; --i) {
            long v = tmp[i];
            if (v == 0)
                continue;
            tmp[i - K] = (tmp[i - K] + v) % MOD;
            tmp[i - K + 1] = (tmp[i - K + 1] + v) % MOD;
        }

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

    public static String solve() {
        long index = 1_000_000_000_000_000_000L;
        if (index < K) {
            return String.valueOf(1 % MOD);
        }

        long[] acc = new long[K];
        acc[0] = 1;

        long[] base = new long[K];
        base[1] = 1;

        long exp = index;
        while (exp > 0) {
            if ((exp & 1) != 0) {
                acc = multiplyMod(acc, base);
            }
            exp >>= 1;
            if (exp > 0) {
                base = multiplyMod(base, base);
            }
        }

        long answer = 0;
        for (long c : acc) {
            answer += c;
        }
        answer %= MOD;

        return String.valueOf(answer);
    }

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