Problem 258: A Lagged Fibonacci Sequence

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 $g_0,g_1,g_2,\dots$ starts with $2000$ ones, $g_k=1$ for $0\le k\le1999$, and from then on each term is the sum of the two terms $2000$ and $1999$ places before it: $g_k=g_{k-2000}+g_{k-1999}$ for $k\ge2000$. We must find the remainder of $g_k$ after division by $20092010$ for $k=10^{18}$; the statement gives no checkpoint values. ^[1] With $2$ in place of $2000$ the rule would give the Fibonacci numbers $1,1,2,3,5,8,\dots$; rules that add two earlier numbers at fixed distances, the lags, modulo a fixed number are known as additive lagged Fibonacci generators and produce pseudorandom numbers. ^[2]

The whole solution in one sentence: in arithmetic with polynomials in which $x^{2000}$ may always be replaced by $x+1$, the power $x^{n}$ becomes a polynomial of degree below $2000$ whose coefficients add up to $g_n$, and $83$ multiplications of such polynomials, with every number reduced modulo $20092010$, produce it for $n=10^{18}$.

Mathematical Approach

We study the first terms, turn the rule into arithmetic with polynomials, prove it correct and then make it fast.

1. The first terms

For $2000\le k\le3998$ both $g_{k-2000}$ and $g_{k-1999}$ are among the first $2000$ ones, so $g_k=2$. Then $g_{3999}=g_{1999}+g_{2000}=3$, and $g_{4000}=g_{2000}+g_{2001}=4$ starts the next stretch of equal terms. The pattern goes on: for every $j$ with $0\le j\le1999$,

$$g_k=2^j\qquad\bigl(2000j\le k\le1999(j+1)\bigr).$$

For $j=0$ this is the definition, and if $k$ lies in the range for $j+1$, then $k-2000$ and $k-1999$ both lie in the range for $j$, so $g_k=2^j+2^j=2^{j+1}$. These plateaus have $2000-j$ terms, and a few terms climb between them: $g_{5998}=5$, $g_{5999}=7$, then $g_{7997}=9$, $g_{7998}=12$, $g_{7999}=15$. For example, $g_{25000}=2^{12}=4096$, because $24000\le25000\le25987$.

The sequence never decreases: for $k\ge2000$ we have $g_{k+1}-g_k=g_{k-1998}-g_{k-2000}$, which is not negative by induction. Hence $g_k\ge2g_{k-2000}$ and $g_k\le2g_{k-1999}$, and induction gives

$$2^{\lfloor k/2000\rfloor}\le g_k\le2^{\lceil k/1999\rceil}.$$

So $g_{10^{18}}$ has between $1.505\cdot10^{14}$ and $1.506\cdot10^{14}$ decimal digits, and $10^{18}$ steps of the rule are out of reach. From $g_{49986}$ on every term exceeds $20092010$, so the remainders differ from the terms.

2. A small model

Replacing $2000$ by $3$ gives a sequence small enough to follow by hand: $h_0=h_1=h_2=1$ and $h_k=h_{k-3}+h_{k-2}$ for $k\ge3$. This is the Padovan sequence $1,1,1,2,2,3,4,5,7,9,12,16,\dots$. ^[3]

3. Terms as combinations of the first terms

In the model, write each term through the starting values: $h_3=h_0+h_1$, $h_4=h_1+h_2$, $h_5=h_2+h_3=h_0+h_1+h_2$, and so on. By induction every term is a sum $c_0h_0+c_1h_1+c_2h_2$: each new term adds the multipliers of two earlier terms, so the multipliers never depend on the starting values; for $h_{10}$ they turn out to be $3$, $5$ and $4$ (Section 7). With all starting values equal to $1$ this gives $h_{10}=12$, and with the starting values $1,0,0$ it gives $h_{10}=3$, as running the rule confirms. In general, a recurrence of order $d$, that is, with longest lag $d$, is fixed by $d$ starting values. ^[4]

The multipliers can be found mechanically. Write $x^i$ as a placeholder for $h_i$. The rule $h_t=h_{t-3}+h_{t-2}$ then says that $x^t$ may be replaced by $x^{t-3}+x^{t-2}$, which is the identity $x^3=x+1$ multiplied by $x^{t-3}$. Starting from $x^{10}$ and replacing the highest power again and again until all exponents are below $3$ ends with $3+5x+4x^2$, the multipliers of $h_{10}$. For our sequence the rule $g_{t+2000}=g_t+g_{t+1}$ gives the replacement

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

The complex numbers work the same way: they are the polynomials in $i$ with $i^2$ replaced by $-1$. In algebra both are quotient rings of a polynomial ring, here by $P(x)=x^{2000}-x-1$. ^[5]

4. Why the remainder gives the term

Replacing the highest term $c\,x^t$, $t\ge2000$, by $c\,x^{t-2000}(x+1)$ subtracts $c\,x^{t-2000}P(x)$ and lowers the degree, so after finitely many replacements

$$\begin{aligned}x^n&=Q(x)\,P(x)+R(x),\\ R(x)&=c_0+c_1x+\dots+c_{1999}x^{1999},\end{aligned}$$

where $Q$ and $R$ have integer coefficients; no fractions occur because the leading coefficient of $P$ is $1$. $R$ is the remainder of $x^n$ after division by $P$.

Claim. $g_n=c_0g_0+c_1g_1+\dots+c_{1999}g_{1999}$. Since the first $2000$ terms are all $1$, $g_n$ is the sum of the coefficients of $R$.

Proof. For a polynomial $A(x)=\sum_ja_jx^j$ let $L(A)=\sum_ja_jg_j$: the map $L$ replaces each power $x^j$ by the term $g_j$. It is linear, $L(A+B)=L(A)+L(B)$ and $L(cA)=c\,L(A)$. For every $j\ge0$ the rule of the sequence with $k=j+2000$ gives

$$L\bigl(x^jP(x)\bigr)=g_{j+2000}-g_{j+1}-g_j=0.$$

Every multiple $Q(x)P(x)$ is a sum of such polynomials $x^jP(x)$ times numbers, so $L(QP)=0$, and therefore $g_n=L(x^n)=L(QP)+L(R)=L(R)$, which is the claim. Combining the coefficients of $x^n\bmod P$ with the starting values in this way is the basis of Fiduccia's algorithm for linear recurrences. ^[6]

The coefficients may only be added at the end. Setting $x=1$ in the rule $x^{2000}=x+1$ would give $1=2$, and indeed coefficient sums do not multiply: $x^{1999}$ and $x$ both have the sum $1$, but their product reduces to $1+x$, whose sum is $2=g_{2000}$. So the program keeps all $2000$ coefficients.

5. Multiplying remainders

For a huge $n$ we never write down $x^n$ itself. If $x^a=Q_1P+A$ and $x^b=Q_2P+B$, then

$$x^{a+b}=\bigl(Q_1Q_2P+Q_1B+Q_2A\bigr)P+AB,$$

so the remainder of $x^{a+b}$ is the remainder of the product $AB$ of the two remainders. If $a_i$ and $b_j$ are the coefficients of $A$ and $B$, the coefficient of $x^t$ in $AB$ is $d_t=\sum_{i+j=t}a_ib_j$, the schoolbook product with up to $2000^2$ terms $a_ib_j$. It has degree at most $3998$, and for $2000\le t\le3998$ the replacement $x^t\to x^{t-2000}+x^{t-1999}$ produces exponents of at most $1999$. So one sweep over the high coefficients finishes the reduction: with $d_{3999}=0$, the reduced product has the coefficients

$$r_i=d_i+d_{i+2000}+d_{i+1999}\qquad(1\le i\le1999)$$

and $r_0=d_0+d_{2000}$. Only additions occur, so the remainder of $x^n$ has non-negative coefficients, none of which exceeds their sum $g_n$.

6. Squaring instead of stepping

Exponentiation by squaring computes $x^n$ from the binary digits of $n$. ^[7] An accumulator starts at $1$ and a base at $x$. The digits are read from the lowest one: if the digit is $1$, the accumulator is multiplied by the base; then the base is squared, so that it runs through $x,x^2,x^4,x^8,\dots$. In the end the accumulator is the product of the powers $x^{2^j}$ over the digits $1$ of $n$, which is $x^n$. Every product is reduced as in Section 5. The number $10^{18}$ has $60$ binary digits, $24$ of them equal to $1$, so the loop needs $59$ squarings, as the base is not squared after the last digit, and $24$ multiplications, one more than necessary because the first one multiplies by $1$: $83$ products in place of $10^{18}$ steps.

7. Worked example in the small model

We compute $h_{10}$ as the program would, with the rule $x^3=x+1$. The number $10$ is $1010$ in binary.

Lowest digit $0$: the accumulator stays $1$; the base becomes $x^2$.
Digit $1$: the accumulator becomes $x^2$; the base becomes $x^4=x\cdot x^3$, which reduces to $x^2+x$.
Digit $0$: the base becomes $(x^2+x)^2=x^4+2x^3+x^2$, which reduces to $(x^2+x)+2(x+1)+x^2$${}=2x^2+3x+2$.
Digit $1$: the accumulator becomes $x^2\,(2x^2+3x+2)$${}=2x^4+3x^3+2x^2$, which reduces to $4x^2+5x+3$. This was the last digit, so the base is not squared again.

The coefficients add up to $3+5+4=12=h_{10}$. On the way, the base $2x^2+3x+2$ stood for $x^8$, and indeed $h_8=2+3+2=7$.

8. Checks with the real rule

The first remainders can be written down directly: $x^{2000}$ and $x^{2001}$ reduce to $1+x$ and $x+x^2$, so $g_{2000}=g_{2001}=2$. Next, $x^{3999}=x^{1999}x^{2000}$ reduces to $x^{2000}+x^{1999}$ and then to $1+x+x^{1999}$, so $g_{3999}=3$; and $x^{4000}=(x^{2000})^2$ reduces to $(1+x)^2=1+2x+x^2$, so $g_{4000}=4$. All four values agree with Section 1.

9. Working modulo $20092010$

Let $M=20092010$. Congruences modulo $M$ can be added and multiplied, ^[8] and the method uses nothing but additions and multiplications of integers, so every coefficient may be replaced by its remainder modulo $M$ as soon as it is formed; the sum of the final coefficients, reduced modulo $M$, is then $g_n\bmod M$. For example, modulo $5$ the remainder $3+5x+4x^2$ of Section 7 becomes $3+4x^2$, whose sum $7$ leaves the remainder $2$, as does $h_{10}=12$. The modulus $M=2\cdot5\cdot859\cdot2339$ is not prime, so modulo $M$ many numbers, for example all even ones, have no inverse, ^[8] but the method never divides. Every stored remainder is below $M\lt2^{25}$, so a product of two of them is below $2^{50}$: too large for $32$-bit integers, but comfortable in $64$ bits.

10. Matrices and faster products

The usual alternative turns the recurrence into a first-order recurrence for the vector of $2000$ consecutive terms, with a $2000\times2000$ matrix, ^[4] and raises that matrix to a power by squaring. One product of two such matrices costs $2000^3=8\cdot10^9$ multiplications, against $2000^2=4\cdot10^6$ for two polynomials. In fact the matrix is, up to transposition, the matrix of multiplication by $x$ in the quotient ring, so each of its powers is fixed by the coefficients of one remainder, and multiplying remainders instead of matrices lowers the cost for a recurrence of order $d$ from about $d^3\log N$ to about $d^2\log N$. ^[6] Products based on the Fast Fourier Transform would need far fewer than $d^2$ operations, ^[6] but for $d=2000$ schoolbook products are fast enough.

How the Code Works

The C++ program ^[9] fixes $K=2000$ and stores a remainder as a Poly, an array of $K$ values of type int whose entry $i$ is the coefficient of $x^i$. parse_arguments can replace the defaults of the Options record, the index $10^{18}$ and the modulus $20092010$, with --index= and --mod=; an unknown argument, or a modulus below $2$, ends the program with code 1. The digit parsers do not check for overflow: an index of $2^{64}$ or more wraps around, and a modulus beyond the range of int overflows.

multiply_mod(a, b, mod) is Section 5. It forms the schoolbook product in a temporary array tmp of $2K$ values of type long long, skipping zero coefficients of either factor and reducing an entry after every addition. Then, for $i$ from $2K-1$ down to $K$, it adds the entry at position $i$ to the positions $i-K$ and $i-K+1$, and it returns the low $K$ entries. Skipping zeros pays off at the start: the powers $x^{2^j}$ are single terms up to $x^{1024}$, $x^{2048}$ reduces to $x^{48}+x^{49}$, and only from $x^{2^{22}}$ on are all $2000$ coefficients non-zero.

solve(index, mod) returns $1$ for an index below $K$. Otherwise it sets acc to $1$ and base to $x$ and runs the loop of Section 6 on exp: it multiplies acc by base when the lowest bit is set, shifts exp right and squares base only while bits remain. Finally it adds the coefficients of acc modulo the modulus, as the claim of Section 4 says. This sum is kept in an int, which is safe for the default modulus but can overflow for a modulus above $2^{30}$ given with --mod=.

Unless --skip-checkpoints is given, run_checkpoints requires $g_0=g_{1999}=1$ and $g_{2000}=g_{2001}=2$ from solve, and compares solve(25000, 20092010) with solve_bruteforce, which runs the rule itself on a window of the last $2000$ terms, shifting the window by one place per step; by Section 1 both must give $2^{12}=4096$. A failed check ends the program with code 2; otherwise it prints the result of solve.

The Python and Java versions have the same multiplication, multiply_mod and multiplyMod, and the same loop, with the index and the modulus fixed and without options or checks. Python's integers have no size limit; Java keeps the coefficients in long arrays and also sums them in a long, reduced once at the end.

Complexity and Verification

The main computation takes $O(K^2\log N)$ time and $O(K)$ memory: for $N=10^{18}$ it makes $83$ products, each with at most $K^2=4\cdot10^6$ multiply-adds and one sweep over the high coefficients, and keeps only a few arrays of at most $2K$ numbers. Because zeros are skipped, it performs $235\,655\,314$ multiply-adds instead of the $332\,000\,000$ of $83$ full products. The C++ checkpoints add $53$ tiny products and the $23001$ steps of solve_bruteforce, which move about $4.6\cdot10^7$ window entries. In our runs the C++ and Java programs took about $0.15$ seconds and the Python program about $16$ seconds.

The checkpoints compare the polynomial arithmetic with the rule, but every number in them is small: up to index $25000$ the terms are at most $4096$, and no coefficient of a product in these runs exceeds $924$, so the reductions modulo $M$ never change anything there. For this article we therefore also verified:

for every $n\le200000$, far beyond index $49986$, the remainder of $x^n$, built step by step by multiplying by $x$ and reducing modulo $P$ and $M$, has coefficients that add up to $g_n\bmod M$ as given by the rule;
the published C++ program with --index= and --mod= agrees with the rule at $75$ indices up to $10^6$ for each of the moduli $20092010$, $10^9+7$, $2^{30}$, $1000$ and $2$; for the moduli $1500000007$ and $2^{31}-1$ the int sum overflowed and some printed values were wrong, while a copy with a 64-bit sum was right;
with $2000$ random starting values, the multipliers of Section 4 reproduce the terms at $7$ indices up to $120000$;
the plateau formula of Section 1 holds for all $0\le j\le1999$, checked modulo two large primes, and the bounds hold exactly for all $k\le60000$;
a separate implementation using a library convolution gives the same remainder for $N=10^{18}$ as the three programs, which all print the same value.

Footnotes and References

The references below supply background facts; the plateau formula, the bounds, the reduction rule for $x^{2000}$ and the operation counts are derived in this article.

Project Euler 258 — A Lagged Fibonacci Sequence. The official statement sets the first 2000 terms of the sequence to 1, makes every later term the sum of the terms 2000 and 1999 places before it, and asks for the term with index 10¹⁸ modulo 20092010.
Wikipedia — Lagged Fibonacci generator. A lagged Fibonacci generator is a pseudorandom number generator, named after the Fibonacci recurrence, that computes S(n) ≡ S(n − j) ⋆ S(n − k) (mod m) with 0 < j < k for a binary operation ⋆; with addition it is called an additive lagged Fibonacci generator.
Wikipedia — Padovan sequence. The Padovan numbers satisfy P(0) = P(1) = P(2) = 1 and P(n) = P(n − 2) + P(n − 3), and they begin 1, 1, 1, 2, 2, 3, 4, 5, 7, 9, 12, 16, ….
Wikipedia — Linear recurrence with constant coefficients. The order of a recurrence is its longest time lag, and that many initial values determine the sequence; a recurrence of order n can be rewritten as a first-order matrix equation.
Wikipedia — Quotient ring. Elements are congruent modulo an ideal when their difference lies in it, and the classes can be added and multiplied; ℝ[X]/(X² + 1) is isomorphic to the complex numbers, X playing the role of i.
A. Bostan and R. Mori — A Simple and Fast Algorithm for Computing the N-th Term of a Linearly Recurrent Sequence (arXiv:2008.08822). Section 1.2: matrix powering for a recurrence of order d costs O(d^θ log N) operations, θ between 2 and 3. Section 1.3, Fiduccia's algorithm (1985): the matrix is the transposed companion matrix of multiplication by x in R[x]/(Γ), Γ the characteristic polynomial, so the N-th term is the inner product of the coefficients of x^N mod Γ with the initial terms, found by binary powering in R[x]/(Γ) with O(log N) products; Fast Fourier Transform methods make each product nearly linear in d.
Wikipedia — Exponentiation by squaring. The method works for any element of a semigroup, such as a polynomial or a square matrix, with about ⌊log₂ n⌋ squarings and one multiplication fewer than the number of ones in the binary expansion of n.
Wikipedia — Modular arithmetic. Congruence modulo m is compatible with addition and multiplication; ℤ/mℤ is a commutative ring, a field only when m is prime, and a number has an inverse modulo m if and only if it is coprime with m.
C++ — Euler258.cpp. The linked immutable C++ revision contains multiply_mod, solve, solve_bruteforce and run_checkpoints; it is the implementation record for this article.

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());
    }
}