Solution

Answer: 12747994

Let

$$g_k = \begin{cases} 1,&0\le k\le 1999,\[4pt] g_{k-2000}+g_{k-1999},&k\ge 2000. \end{cases}$$

We must compute

$$g_{10^{18}} \pmod{20092010}.$$

1. Mathematical analysis

This is a linear recurrence of order $2000$.

Define the characteristic relation:

$$g_n-g_{n-1999}-g_{n-2000}=0.$$

The characteristic polynomial is therefore

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

The key fact for linear recurrences is:

Every term $g_n$ can be expressed as a linear combination of the first $2000$ terms.

Since all initial values are

$$g_0=g_1=\cdots=g_{1999}=1,$$

if we can compute the coefficients of $x^n$ modulo the characteristic polynomial $P(x)$, then the answer is simply the sum of those coefficients.

Polynomial viewpoint

Inside the quotient ring

$$\mathbb Z_m[x]/(x^{2000}-x-1),$$

we have the reduction rule

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

Thus any power $x^n$ can be reduced to a polynomial of degree $<2000$:

$$x^n \equiv \sum_{i=0}^{1999} c_i x^i.$$

Then

$$g_n=\sum_{i=0}^{1999} c_i g_i.$$

But every $g_i=1$, so

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

Therefore the task reduces to:

Compute $x^{10^{18}} \bmod (x^{2000}-x-1)$.
Sum its coefficients modulo $20092010$.

2. Efficient algorithm

We use binary exponentiation on polynomials.

When multiplying two polynomials, degrees can go up to $3998$.

Whenever a term $x^d$ with $d\ge 2000$ appears, reduce using

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

Thus

$$x^d = x^{d-2000}x^{2000} \equiv x^{d-1999}+x^{d-2000}.$$

This gives an $O(K^2\log n)$ algorithm with $K=2000$, which is easily fast enough.

3. Python implementation

MOD = 20092010
K = 2000
N = 10**18

def combine(a, b):
    """
    Multiply two polynomials modulo:
        x^2000 - x - 1
    Result has degree < 2000.
    """
    tmp = [0] * (2 * K - 1)

    # Ordinary polynomial multiplication
    for i, ai in enumerate(a):
        if ai:
            for j, bj in enumerate(b):
                if bj:
                    tmp[i + j] = (tmp[i + j] + ai * bj) % MOD

    # Reduce high-degree terms using:
    # x^2000 = x + 1
    for d in range(2 * K - 2, K - 1, -1):
        v = tmp[d]
        if v:
            tmp[d - K] = (tmp[d - K] + v) % MOD
            tmp[d - K + 1] = (tmp[d - K + 1] + v) % MOD

    return tmp[:K]

def power_x(n):
    """
    Compute x^n modulo the characteristic polynomial.
    Return coefficients of degree < 2000 polynomial.
    """

    # res = 1
    res = [1] + [0] * (K - 1)

    # base = x
    base = [0, 1] + [0] * (K - 2)

    while n > 0:
        if n & 1:
            res = combine(res, base)

        base = combine(base, base)
        n >>= 1

    return res

# coefficients for x^N
coeffs = power_x(N)

# since g_0 ... g_1999 are all 1,
# g_N is the sum of coefficients
answer = sum(coeffs) % MOD

print(answer)

4. Code walkthrough

`combine(a, b)`

This function multiplies two reduced polynomials.

First:

tmp[i + j] += ai * bj

performs ordinary polynomial multiplication.

Then we reduce terms of degree $\ge 2000$.

Since

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

we replace

$$x^d \equiv x^{d-1999}+x^{d-2000}.$$

So coefficient v of $x^d$ is distributed to:

tmp[d - K]
tmp[d - K + 1]

which correspond to

$$x^{d-2000},\quad x^{d-1999}.$$

`power_x(n)`

This computes

$$x^n \bmod (x^{2000}-x-1)$$

using binary exponentiation.

res starts as $1$
base starts as $x$

Whenever a binary bit of $n$ is set:

res = combine(res, base)

Then square the base each step.

5. Small sanity checks

From the recurrence:

$$g_{2000}=g_0+g_1=2,$$

$$g_{2001}=g_1+g_2=2.$$

The program reproduces these correctly.

All arithmetic is performed modulo $20092010$, exactly as required.

After carrying out the computation for $n=10^{18}$, we obtain:

Answer: 12747994

Project Euler Problem 258