Problem 999: Alternating Recurrence
View on Project EulerProject Euler Problem 999 Solution
EulerSolve provides an optimized solution for Project Euler Problem 999, Alternating Recurrence, with C++, Python, Java, and a step-by-step mathematical explanation.
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....
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
- 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.
- 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\).
- 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.
Problem 999 source code
C++
#include <array>
#include <cassert>
#include <cstdint>
#include <iostream>
#include <vector>
namespace {
using i64 = std::int64_t;
using u64 = std::uint64_t;
using u128 = unsigned __int128;
constexpr u64 MOD = 1'234'567'891ULL;
constexpr u64 INV2 = (MOD + 1) / 2;
constexpr u64 TARGET_N = 1'000'000'000'000'000'003ULL;
u64 normalize(const i64 value) {
i64 result = value % static_cast<i64>(MOD);
if (result < 0) {
result += static_cast<i64>(MOD);
}
return static_cast<u64>(result);
}
u64 sub_mod(const u64 lhs, const u64 rhs) {
return lhs >= rhs ? lhs - rhs : lhs + MOD - rhs;
}
u64 mul_mod(const u64 lhs, const u64 rhs) {
return static_cast<u64>((static_cast<u128>(lhs) * rhs) % MOD);
}
u64 pow_mod(u64 base, u64 exponent) {
u64 result = 1;
while (exponent > 0) {
if ((exponent & 1) != 0) {
result = mul_mod(result, base);
}
base = mul_mod(base, base);
exponent >>= 1;
}
return result;
}
struct Field {
std::array<u64, 4> c{};
Field() = default;
explicit Field(const i64 c0) : c{normalize(c0), 0, 0, 0} {}
Field(const i64 c0, const i64 c1, const i64 c2, const i64 c3)
: c{normalize(c0), normalize(c1), normalize(c2), normalize(c3)} {}
bool is_base() const {
return c[1] == 0 && c[2] == 0 && c[3] == 0;
}
bool is_rho_multiple() const {
return c[0] == 0 && c[2] == 0 && c[3] == 0;
}
};
Field operator-(const Field& lhs, const Field& rhs) {
Field result;
for (int i = 0; i < 4; ++i) {
result.c[static_cast<std::size_t>(i)] =
sub_mod(lhs.c[static_cast<std::size_t>(i)], rhs.c[static_cast<std::size_t>(i)]);
}
return result;
}
Field operator-(const Field& value) {
Field result;
for (int i = 0; i < 4; ++i) {
const u64 current = value.c[static_cast<std::size_t>(i)];
result.c[static_cast<std::size_t>(i)] = current == 0 ? 0 : MOD - current;
}
return result;
}
Field operator*(const Field& lhs, const Field& rhs) {
std::array<u128, 7> product{};
for (int i = 0; i < 4; ++i) {
for (int j = 0; j < 4; ++j) {
product[static_cast<std::size_t>(i + j)] +=
static_cast<u128>(lhs.c[static_cast<std::size_t>(i)]) *
rhs.c[static_cast<std::size_t>(j)];
}
}
for (int i = 6; i >= 4; --i) {
product[static_cast<std::size_t>(i - 4)] += 2 * product[static_cast<std::size_t>(i)];
}
Field result;
for (int i = 0; i < 4; ++i) {
result.c[static_cast<std::size_t>(i)] =
static_cast<u64>(product[static_cast<std::size_t>(i)] % MOD);
}
return result;
}
bool operator==(const Field& lhs, const Field& rhs) {
return lhs.c == rhs.c;
}
Field square(const Field& value) {
return value * value;
}
Field cube(const Field& value) {
return value * value * value;
}
using Block = std::array<Field, 9>;
const Field RHO(0, 1, 0, 0);
const Field INV_RHO(0, 0, 0, static_cast<i64>(INV2));
Field get_value(const Block& block, const i64 center, const i64 index) {
const i64 offset = index - center;
assert(-4 <= offset && offset <= 4);
return block[static_cast<std::size_t>(offset + 4)];
}
Field odd_value(const Block& block, const i64 center, const i64 k) {
return get_value(block, center, k + 1) * cube(get_value(block, center, k - 1)) -
get_value(block, center, k - 2) * cube(get_value(block, center, k));
}
Field even_value(const Block& block, const i64 center, const i64 k) {
return get_value(block, center, k) * INV_RHO *
(get_value(block, center, k + 2) * square(get_value(block, center, k - 1)) -
get_value(block, center, k - 2) * square(get_value(block, center, k + 1)));
}
void advance(Block& block, i64& center, const int bit) {
const Block previous = block;
const i64 previous_center = center;
const i64 next_center = 2 * center + bit;
Block next{};
for (int offset = -4; offset <= 4; ++offset) {
const i64 index = next_center + offset;
if (index % 2 == 0) {
next[static_cast<std::size_t>(offset + 4)] = even_value(previous, previous_center, index / 2);
} else {
next[static_cast<std::size_t>(offset + 4)] = odd_value(previous, previous_center, (index + 1) / 2);
}
}
block = next;
center = next_center;
}
Field eds_value(const u64 n) {
assert(n > 0);
Block block{Field(1), -RHO, Field(-1), Field(0), Field(1), RHO, Field(-1), Field(0, -2, 0, 0), Field(-3)};
i64 center = 1;
int highest = 63;
while (((n >> highest) & 1ULL) == 0) {
--highest;
}
for (int bit_index = highest - 1; bit_index >= 0; --bit_index) {
advance(block, center, static_cast<int>((n >> bit_index) & 1ULL));
}
return block[4];
}
u64 sequence_value(const u64 n) {
const Field value = eds_value(n);
if ((n & 1ULL) != 0) {
assert(value.is_base());
if ((((n - 1) / 2) & 1ULL) == 0) {
return value.c[0];
}
return value.c[0] == 0 ? 0 : MOD - value.c[0];
}
assert(value.is_rho_multiple());
if ((((n - 2) / 2) & 1ULL) == 0) {
return value.c[1];
}
return value.c[1] == 0 ? 0 : MOD - value.c[1];
}
u64 direct_value(const int n) {
std::vector<u64> a(static_cast<std::size_t>(n + 5), 0);
a[1] = 1;
a[2] = 1;
a[3] = 1;
a[4] = 2;
for (int k = 3; k + 2 <= n; ++k) {
const u64 u = (k % 2 == 0) ? 1 : 2;
const u64 left = mul_mod(a[static_cast<std::size_t>(k)], a[static_cast<std::size_t>(k)]);
const u64 right = mul_mod(u, mul_mod(a[static_cast<std::size_t>(k + 1)], a[static_cast<std::size_t>(k - 1)]));
const u64 denominator = a[static_cast<std::size_t>(k - 2)];
assert(denominator != 0);
a[static_cast<std::size_t>(k + 2)] = mul_mod(sub_mod(left, right), pow_mod(denominator, MOD - 2));
}
return a[static_cast<std::size_t>(n)];
}
void run_checkpoints() {
assert(square(square(RHO)) == Field(2));
for (int n = 1; n <= 200; ++n) {
assert(sequence_value(static_cast<u64>(n)) == direct_value(n));
}
assert(sequence_value(13) == 23'321);
assert(direct_value(1003) == 231'906'014);
assert(sequence_value(1003) == 231'906'014);
}
} // namespace
int main() {
run_checkpoints();
std::cout << sequence_value(TARGET_N) << '\n';
return 0;
}
Python
MOD = 1_234_567_891
INV2 = (MOD + 1) // 2
TARGET_N = 1_000_000_000_000_000_003
def normalize(value):
return value % MOD
def sub_mod(lhs, rhs):
return lhs - rhs if lhs >= rhs else lhs + MOD - rhs
def mul_mod(lhs, rhs):
return (lhs * rhs) % MOD
class Field:
__slots__ = ("c",)
def __init__(self, c0=0, c1=0, c2=0, c3=0):
self.c = (
normalize(c0),
normalize(c1),
normalize(c2),
normalize(c3),
)
def is_base(self):
return self.c[1] == 0 and self.c[2] == 0 and self.c[3] == 0
def is_rho_multiple(self):
return self.c[0] == 0 and self.c[2] == 0 and self.c[3] == 0
def __neg__(self):
return Field(*(0 if x == 0 else MOD - x for x in self.c))
def __sub__(self, other):
return Field(*(sub_mod(self.c[i], other.c[i]) for i in range(4)))
def __mul__(self, other):
product = [0] * 7
for i in range(4):
for j in range(4):
product[i + j] = (product[i + j] + self.c[i] * other.c[j]) % MOD
for i in range(6, 3, -1):
product[i - 4] = (product[i - 4] + 2 * product[i]) % MOD
return Field(product[0], product[1], product[2], product[3])
def __eq__(self, other):
return self.c == other.c
RHO = Field(0, 1, 0, 0)
INV_RHO = Field(0, 0, 0, INV2)
def square(value):
return value * value
def cube(value):
return value * value * value
def get_value(block, center, index):
offset = index - center
assert -4 <= offset <= 4
return block[offset + 4]
def odd_value(block, center, k):
return (
get_value(block, center, k + 1) * cube(get_value(block, center, k - 1))
- get_value(block, center, k - 2) * cube(get_value(block, center, k))
)
def even_value(block, center, k):
return (
get_value(block, center, k)
* INV_RHO
* (
get_value(block, center, k + 2)
* square(get_value(block, center, k - 1))
- get_value(block, center, k - 2)
* square(get_value(block, center, k + 1))
)
)
def advance(block, center, bit):
previous = block
previous_center = center
next_center = 2 * center + bit
next_block = []
for offset in range(-4, 5):
index = next_center + offset
if index % 2 == 0:
next_block.append(even_value(previous, previous_center, index // 2))
else:
next_block.append(odd_value(previous, previous_center, (index + 1) // 2))
return next_block, next_center
def eds_value(n):
assert n > 0
block = [
Field(1),
-RHO,
Field(-1),
Field(0),
Field(1),
RHO,
Field(-1),
Field(0, -2, 0, 0),
Field(-3),
]
center = 1
for bit_index in range(n.bit_length() - 2, -1, -1):
block, center = advance(block, center, (n >> bit_index) & 1)
return block[4]
def sequence_value(n):
value = eds_value(n)
if n & 1:
assert value.is_base()
if (((n - 1) // 2) & 1) == 0:
return value.c[0]
return 0 if value.c[0] == 0 else MOD - value.c[0]
assert value.is_rho_multiple()
if (((n - 2) // 2) & 1) == 0:
return value.c[1]
return 0 if value.c[1] == 0 else MOD - value.c[1]
def direct_value(n):
a = [0] * (n + 5)
a[1] = 1
a[2] = 1
a[3] = 1
a[4] = 2
for k in range(3, n - 1):
u = 1 if k % 2 == 0 else 2
left = mul_mod(a[k], a[k])
right = mul_mod(u, mul_mod(a[k + 1], a[k - 1]))
denominator = a[k - 2]
assert denominator != 0
a[k + 2] = mul_mod(sub_mod(left, right), pow(denominator, MOD - 2, MOD))
return a[n]
def run_checkpoints():
assert square(square(RHO)) == Field(2)
for n in range(1, 201):
assert sequence_value(n) == direct_value(n)
assert sequence_value(13) == 23_321
assert direct_value(1003) == 231_906_014
assert sequence_value(1003) == 231_906_014
if __name__ == "__main__":
run_checkpoints()
print(sequence_value(TARGET_N))
Java
public class Euler999 {
private static final long MOD = 1_234_567_891L;
private static final long INV2 = (MOD + 1) / 2;
private static final long TARGET_N = 1_000_000_000_000_000_003L;
private static long normalize(long value) {
long result = value % MOD;
return result < 0 ? result + MOD : result;
}
private static long subMod(long lhs, long rhs) {
return lhs >= rhs ? lhs - rhs : lhs + MOD - rhs;
}
private static long mulMod(long lhs, long rhs) {
return (lhs * rhs) % MOD;
}
private static long powMod(long base, long exponent) {
long result = 1;
while (exponent > 0) {
if ((exponent & 1L) != 0) {
result = mulMod(result, base);
}
base = mulMod(base, base);
exponent >>= 1;
}
return result;
}
private static final class Field {
final long[] c;
Field(long c0) {
this(c0, 0, 0, 0);
}
Field(long c0, long c1, long c2, long c3) {
c = new long[] {normalize(c0), normalize(c1), normalize(c2), normalize(c3)};
}
boolean isBase() {
return c[1] == 0 && c[2] == 0 && c[3] == 0;
}
boolean isRhoMultiple() {
return c[0] == 0 && c[2] == 0 && c[3] == 0;
}
Field subtract(Field other) {
return new Field(
subMod(c[0], other.c[0]),
subMod(c[1], other.c[1]),
subMod(c[2], other.c[2]),
subMod(c[3], other.c[3])
);
}
Field negate() {
return new Field(
c[0] == 0 ? 0 : MOD - c[0],
c[1] == 0 ? 0 : MOD - c[1],
c[2] == 0 ? 0 : MOD - c[2],
c[3] == 0 ? 0 : MOD - c[3]
);
}
Field multiply(Field other) {
long[] product = new long[7];
for (int i = 0; i < 4; ++i) {
for (int j = 0; j < 4; ++j) {
product[i + j] = (product[i + j] + mulMod(c[i], other.c[j])) % MOD;
}
}
for (int i = 6; i >= 4; --i) {
product[i - 4] = (product[i - 4] + 2 * product[i]) % MOD;
}
return new Field(product[0], product[1], product[2], product[3]);
}
@Override
public boolean equals(Object obj) {
if (!(obj instanceof Field)) {
return false;
}
Field other = (Field) obj;
for (int i = 0; i < 4; ++i) {
if (c[i] != other.c[i]) {
return false;
}
}
return true;
}
@Override
public int hashCode() {
int result = 17;
for (long value : c) {
result = 31 * result + Long.hashCode(value);
}
return result;
}
}
private static final Field RHO = new Field(0, 1, 0, 0);
private static final Field INV_RHO = new Field(0, 0, 0, INV2);
private static Field square(Field value) {
return value.multiply(value);
}
private static Field cube(Field value) {
return value.multiply(value).multiply(value);
}
private static Field getValue(Field[] block, long center, long index) {
long offset = index - center;
if (offset < -4 || offset > 4) {
throw new AssertionError("Block lookup out of range");
}
return block[(int) offset + 4];
}
private static Field oddValue(Field[] block, long center, long k) {
return getValue(block, center, k + 1)
.multiply(cube(getValue(block, center, k - 1)))
.subtract(getValue(block, center, k - 2).multiply(cube(getValue(block, center, k))));
}
private static Field evenValue(Field[] block, long center, long k) {
return getValue(block, center, k)
.multiply(INV_RHO)
.multiply(
getValue(block, center, k + 2)
.multiply(square(getValue(block, center, k - 1)))
.subtract(
getValue(block, center, k - 2)
.multiply(square(getValue(block, center, k + 1)))
)
);
}
private static long advance(Field[] block, long center, int bit) {
Field[] previous = block.clone();
long previousCenter = center;
long nextCenter = 2 * center + bit;
for (int offset = -4; offset <= 4; ++offset) {
long index = nextCenter + offset;
if ((index & 1L) == 0) {
block[offset + 4] = evenValue(previous, previousCenter, index / 2);
} else {
block[offset + 4] = oddValue(previous, previousCenter, (index + 1) / 2);
}
}
return nextCenter;
}
private static Field edsValue(long n) {
if (n <= 0) {
throw new AssertionError("Index must be positive");
}
Field[] block = {
new Field(1),
RHO.negate(),
new Field(-1),
new Field(0),
new Field(1),
RHO,
new Field(-1),
new Field(0, -2, 0, 0),
new Field(-3)
};
long center = 1;
int highest = 63 - Long.numberOfLeadingZeros(n);
for (int bitIndex = highest - 1; bitIndex >= 0; --bitIndex) {
center = advance(block, center, (int) ((n >> bitIndex) & 1L));
}
return block[4];
}
private static long sequenceValue(long n) {
Field value = edsValue(n);
if ((n & 1L) != 0) {
if (!value.isBase()) {
throw new AssertionError("Odd index should be in the base field");
}
if ((((n - 1) / 2) & 1L) == 0) {
return value.c[0];
}
return value.c[0] == 0 ? 0 : MOD - value.c[0];
}
if (!value.isRhoMultiple()) {
throw new AssertionError("Even index should be a rho multiple");
}
if ((((n - 2) / 2) & 1L) == 0) {
return value.c[1];
}
return value.c[1] == 0 ? 0 : MOD - value.c[1];
}
private static long directValue(int n) {
long[] a = new long[n + 5];
a[1] = 1;
a[2] = 1;
a[3] = 1;
a[4] = 2;
for (int k = 3; k + 2 <= n; ++k) {
long u = (k % 2 == 0) ? 1 : 2;
long left = mulMod(a[k], a[k]);
long right = mulMod(u, mulMod(a[k + 1], a[k - 1]));
long denominator = a[k - 2];
if (denominator == 0) {
throw new AssertionError("Unexpected zero denominator");
}
a[k + 2] = mulMod(subMod(left, right), powMod(denominator, MOD - 2));
}
return a[n];
}
private static void check(boolean condition) {
if (!condition) {
throw new AssertionError();
}
}
private static void runCheckpoints() {
check(square(square(RHO)).equals(new Field(2)));
for (int n = 1; n <= 200; ++n) {
check(sequenceValue(n) == directValue(n));
}
check(sequenceValue(13) == 23_321);
check(directValue(1003) == 231_906_014);
check(sequenceValue(1003) == 231_906_014);
}
public static void main(String[] args) {
runCheckpoints();
System.out.println(sequenceValue(TARGET_N));
}
}