Detailed mathematical approach

Problem Summary

The sequence starts with \(a_1=a_2=a_3=1\) and \(a_4=2\). For every later step it is constrained by

\[ a_n^2=a_{n+2}a_{n-2}+u\,a_{n+1}a_{n-1}, \]

where \(u=1\) for even \(n\), and \(u=2\) for odd \(n\). The target index is enormous, \(n=10^{18}+3\), so iterating the recurrence is impossible. The task is to compute \(a_n\) modulo

\[ M=1234567891. \]

The implementation turns the alternating coefficient \(1,2,1,2,\dots\) into a constant-coefficient elliptic divisibility sequence over a small algebraic extension of \(\mathbb{F}_M\), then uses binary doubling to jump directly to the requested index.

Mathematical Approach: Removing the Alternation

The obstacle in the recurrence is the parity-dependent coefficient. Introduce a formal element \(\rho\) satisfying

\[ \rho^4=2. \]

All computations are done modulo \(M\), in the four-dimensional algebra

\[ R=\mathbb{F}_M[\rho]/(\rho^4-2). \]

An element of \(R\) is stored as

\[ c_0+c_1\rho+c_2\rho^2+c_3\rho^3. \]

Multiplication is ordinary polynomial multiplication followed by the reduction \(\rho^4=2\). Thus every term \(d_i\rho^i\) with \(i\ge4\) is folded into \(2d_i\rho^{i-4}\). This is exactly what the Field multiplication in the code implements.

The Normalized Sequence

Define a new sequence \(W_n\) by absorbing both the sign pattern and the parity factor:

\[ W_n= \begin{cases} (-1)^{(n-1)/2}a_n, & n \text{ odd},\\[2mm] (-1)^{(n-2)/2}a_n\rho, & n \text{ even}. \end{cases} \]

The first values become

\[ W_1=1,\qquad W_2=\rho,\qquad W_3=-1,\qquad W_4=-2\rho,\qquad W_5=-3. \]

With this normalization the original alternating recurrence is equivalent to the constant-coefficient relation

\[ W_{n+2}W_{n-2}=W_n^2+\rho^2W_{n+1}W_{n-1}. \]

For odd \(n\), the product \(W_{n+1}W_{n-1}\) contains \(\rho^2\), so the extra factor \(\rho^2\) contributes \(\rho^4=2\), matching the coefficient \(u=2\). For even \(n\), the common \(\rho^2\) factor cancels from the other terms, leaving coefficient \(u=1\). This is the key transformation: the alternation has not disappeared by approximation; it has been encoded exactly into the algebra.

Elliptic Divisibility Sequence Identities

The sequence \(W_n\) is handled as an elliptic divisibility sequence. The code uses the standard doubling identities

\[ W_{2k-1}=W_{k+1}W_{k-1}^3-W_{k-2}W_k^3, \]

and

\[ W_{2k}=\frac{W_k}{W_2}\left(W_{k+2}W_{k-1}^2-W_{k-2}W_{k+1}^2\right). \]

Here \(W_2=\rho\), and since \(\rho^4=2\),

\[ \rho^{-1}=\frac{\rho^3}{2}. \]

Modulo \(M\), division by \(2\) is multiplication by \((M+1)/2\), so the code stores \(\rho^{-1}\) as INV_RHO. These identities let us compute values around index \(2k\) from a small block of values around \(k\).

Maintaining a Nine-Term Window

The implementation keeps a block

\[ W_{c-4},W_{c-3},\dots,W_c,\dots,W_{c+4} \]

centered at an index \(c\). Initially \(c=1\), and the stored block is

\[ W_{-3}=1,\;W_{-2}=-\rho,\;W_{-1}=-1,\;W_0=0,\;W_1=1,\;W_2=\rho,\;W_3=-1,\;W_4=-2\rho,\;W_5=-3. \]

When the next binary bit of the target index is \(b\in\{0,1\}\), the center moves to

\[ c' = 2c+b. \]

For each offset \(-4\le r\le4\), the code computes \(W_{c'+r}\). If \(c'+r\) is odd it applies the \(W_{2k-1}\) formula; if it is even it applies the \(W_{2k}\) formula. Both formulas require only values with indices from \(k-2\) through \(k+2\), which are guaranteed to lie inside the previous nine-term block. This is why a constant-size window is enough for an index with 60 binary bits.

Recovering \(a_n\)

After all bits have been processed, the center of the block is exactly the requested \(n\), and block[4] contains \(W_n\). The original sequence value is recovered by undoing the normalization:

\[ a_n= \begin{cases} (-1)^{(n-1)/2}W_n, & n \text{ odd},\\[2mm] (-1)^{(n-2)/2}\,[\rho]W_n, & n \text{ even}, \end{cases} \]

where \([\rho]W_n\) denotes the coefficient of \(\rho\). The implementation asserts that odd-indexed values are in the base field and that even-indexed values are pure multiples of \(\rho\). These assertions are useful guards against an incorrect field operation or a broken doubling step.

Correctness Argument

Algebraic soundness. The definition of \(W_n\) converts the original parity-dependent recurrence into one recurrence in \(R\). Because \(\rho^4=2\), the coefficient is exactly \(2\) on odd steps and exactly \(1\) on even steps after the normalization is undone.

Doubling soundness. The two identities used by odd_value and even_value are elliptic divisibility sequence identities. Each call to advance computes the same \(W_j\) values that would be obtained by the recurrence, but for indices whose binary prefix has just been extended.

Window completeness. Every requested value in the next block depends only on five neighboring values in the previous block. Hence the nine-term window always contains all data needed for the next binary bit, and no earlier terms are needed.

Validation. The program checks \(\rho^4=2\), compares the fast method with a direct recurrence implementation for \(1\le n\le200\), and verifies the published checkpoints \(a_{13}=23321\) and \(a_{1003}\equiv231906014\pmod M\).

How the Code Works

Field represents \(c_0+c_1\rho+c_2\rho^2+c_3\rho^3\). Addition and subtraction are componentwise modulo \(M\); multiplication first builds seven polynomial coefficients and then folds the high-degree terms back with \(\rho^4=2\).

advance reads one bit of the target index. It builds the next nine-term block by choosing between the odd and even EDS doubling formulas. eds_value scans the binary expansion of \(n\) from most significant bit to least significant bit, so the number of block updates is \(O(\log n)\).

sequence_value converts \(W_n\) back into \(a_n\), applying the sign pattern and selecting either the scalar coefficient or the \(\rho\)-coefficient depending on parity.

Why Modular Computation Is Safe

The fast path never divides by a term of the original recurrence. It only divides by \(W_2=\rho\), whose inverse is known explicitly as \(\rho^3/2\). Since \(M\) is odd, \(2\) has the inverse \((M+1)/2\) modulo \(M\). Therefore the EDS doubling step is a sequence of additions, subtractions, multiplications, and one fixed multiplication by \(\rho^{-1}\).

The direct recurrence used in the checkpoint routine does divide by \(a_{k-2}\), but it is not part of the large-index algorithm. It is only a validation tool for small \(n\), and the code asserts that the denominator is nonzero before using Fermat inversion. The production computation for \(10^{18}+3\) does not rely on those divisions.

Binary Index Invariant

After processing a binary prefix \(p\) of the target index, the stored block is \(W_{p-4},W_{p-3},\ldots,W_{p+4}\); equivalently, its center is \(p\). The initial prefix is the leading bit \(1\), so the initial block is centered at \(1\). Appending a new bit \(b\) changes the represented integer from \(p\) to \(2p+b\), exactly matching the update \(c'=2c+b\). Thus, after the last bit, the center is the full target index. This is the same structural idea as binary exponentiation, but the object being doubled is a local EDS window instead of a single number.

Complexity Analysis

Each binary bit performs a constant number of field multiplications and subtractions on four coefficients. Therefore the running time is

\[ O(\log n). \]

The memory usage is \(O(1)\), since only two nine-term blocks and a few temporary field elements are needed. For \(n=10^{18}+3\), this means roughly sixty block updates rather than \(10^{18}\) recurrence steps.

References

1. Problem page: Project Euler 999. The statement supplies the alternating recurrence, the initial values, the two checkpoint values, and the modulus \(M=1234567891\) used throughout the solution.
2. Elliptic divisibility sequence: Wikipedia - Elliptic divisibility sequence. This is the source framework for the identities for \(W_{2k-1}\) and \(W_{2k}\), which allow a nine-term window centered at \(k\) to be transformed into one centered at \(2k\) or \(2k+1\).
3. Finite field: Wikipedia - Finite field. The algorithm performs all coefficient arithmetic modulo \(M\) inside the algebra \(R=\mathbb{F}_M[\rho]/(\rho^4-2)\), so finite-field arithmetic justifies reducing after every addition, subtraction, and multiplication.