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