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.
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.
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.
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\).
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.
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.
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\).
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.
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.
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.
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.
Die Folge beginnt mit \(a_1=a_2=a_3=1\) und \(a_4=2\). Für alle weiteren Indizes gilt die Bedingung
\[ a_n^2=a_{n+2}a_{n-2}+u\,a_{n+1}a_{n-1}, \]
wobei \(u=1\) für gerade \(n\) und \(u=2\) für ungerade \(n\) ist. Der Zielindex ist riesig, nämlich \(n=10^{18}+3\), und gesucht ist \(a_n\) modulo
\[ M=1234567891. \]
Eine lineare Iteration der Rekurrenz ist damit ausgeschlossen. Die Lösung übersetzt den alternierenden Koeffizienten \(1,2,1,2,\dots\) exakt in eine konstante elliptische Divisibilitätsfolge über einer kleinen algebraischen Erweiterung von \(\mathbb{F}_M\), und springt dann mit binärer Verdopplung direkt zum Zielindex.
Das Hindernis ist der von der Parität abhängige Koeffizient. Dafür führen wir ein formales Element \(\rho\) ein, das
\[ \rho^4=2 \]
erfüllt. Alle Rechnungen finden modulo \(M\) in der vierdimensionalen Algebra
\[ R=\mathbb{F}_M[\rho]/(\rho^4-2). \]
statt. Ein Element wird als
\[ c_0+c_1\rho+c_2\rho^2+c_3\rho^3 \]
gespeichert. Multiplikation ist zuerst gewöhnliche Polynommultiplikation; danach werden alle Terme \(\rho^i\) mit \(i\ge4\) über \(\rho^4=2\) zurückgefaltet. Ein Term \(d_i\rho^i\) wird also zu \(2d_i\rho^{i-4}\). Genau diese Reduktion implementiert die Klasse Field im Code.
Wir definieren eine neue Folge \(W_n\), die sowohl das Vorzeichenmuster als auch den Paritätsfaktor aufnimmt:
\[ W_n= \begin{cases} (-1)^{(n-1)/2}a_n, & n \text{ ungerade},\\[2mm] (-1)^{(n-2)/2}a_n\rho, & n \text{ gerade}. \end{cases} \]
Die ersten Werte sind dann
\[ W_1=1,\qquad W_2=\rho,\qquad W_3=-1,\qquad W_4=-2\rho,\qquad W_5=-3. \]
Mit dieser Normalisierung wird die ursprüngliche alternierende Rekurrenz zur konstanten Relation
\[ W_{n+2}W_{n-2}=W_n^2+\rho^2W_{n+1}W_{n-1}. \]
Für ungerades \(n\) tragen \(W_{n+1}\) und \(W_{n-1}\) jeweils einen Faktor \(\rho\). Zusammen mit dem äußeren \(\rho^2\) entsteht \(\rho^4=2\), also genau der Koeffizient \(u=2\). Für gerades \(n\) hebt sich der gemeinsame \(\rho^2\)-Faktor in den anderen Termen heraus, und der Koeffizient bleibt \(u=1\). Die Alternation wird also nicht angenähert, sondern exakt in der Algebra kodiert.
Die Folge \(W_n\) wird als elliptische Divisibilitätsfolge behandelt. Der Code nutzt die Standardidentitäten
\[ W_{2k-1}=W_{k+1}W_{k-1}^3-W_{k-2}W_k^3, \]
und
\[ W_{2k}=\frac{W_k}{W_2}\left(W_{k+2}W_{k-1}^2-W_{k-2}W_{k+1}^2\right). \]
Hier ist \(W_2=\rho\). Aus \(\rho^4=2\) folgt
\[ \rho^{-1}=\frac{\rho^3}{2}. \]
Modulo \(M\) bedeutet Division durch \(2\) Multiplikation mit \((M+1)/2\). Deshalb speichert der Code \(\rho^{-1}\) als INV_RHO. Die beiden Identitäten berechnen Werte um \(2k\) aus einem kleinen Block von Werten um \(k\).
Die Implementierung hält stets den Block
\[ W_{c-4},W_{c-3},\ldots,W_c,\ldots,W_{c+4} \]
um ein Zentrum \(c\). Anfangs ist \(c=1\), und gespeichert werden
\[ 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. \]
Wenn das nächste Bit des Zielindex \(b\in\{0,1\}\) ist, wandert das Zentrum nach
\[ c'=2c+b. \]
Für jeden Offset \(-4\le r\le4\) berechnet der Code \(W_{c'+r}\). Ist \(c'+r\) ungerade, wird die Formel für \(W_{2k-1}\) benutzt; ist der Index gerade, wird die Formel für \(W_{2k}\) benutzt. Beide Formeln benötigen nur die Werte \(k-2,\ldots,k+2\), und diese liegen garantiert im alten Neun-Terme-Fenster. Daher reicht ein Fenster konstanter Größe für einen Index mit etwa sechzig Binärstellen.
Nachdem alle Bits verarbeitet wurden, ist das Zentrum des Blocks genau der gesuchte Index \(n\), und block[4] enthält \(W_n\). Der ursprüngliche Wert wird durch Umkehrung der Normalisierung gewonnen:
\[ a_n= \begin{cases} (-1)^{(n-1)/2}W_n, & n \text{ ungerade},\\[2mm] (-1)^{(n-2)/2}[\rho]W_n, & n \text{ gerade}. \end{cases} \]
Dabei bezeichnet \([\rho]W_n\) den Koeffizienten von \(\rho\). Die Implementierung prüft, dass ungerade Indizes wirklich im Grundkörper liegen und gerade Indizes reine \(\rho\)-Vielfache sind. Diese Assertions sind nützliche Schutzmechanismen gegen Fehler in der Ringarithmetik oder in der Verdopplung.
Algebraische Korrektheit. Die Definition von \(W_n\) verwandelt die paritätsabhängige Rekurrenz in eine einzige Rekurrenz in \(R\). Wegen \(\rho^4=2\) entsteht nach Rückumwandlung genau der Faktor \(2\) bei ungeraden und der Faktor \(1\) bei geraden Indizes.
Korrektheit der Verdopplung. Die Funktionen odd_value und even_value implementieren Identitäten elliptischer Divisibilitätsfolgen. Jeder Aufruf von advance erzeugt daher dieselben \(W_j\)-Werte wie die Rekurrenz, nur für den gerade erweiterten binären Präfix.
Vollständigkeit des Fensters. Jeder neue Wert hängt nur von fünf benachbarten alten Werten ab. Das Neun-Terme-Fenster enthält also immer alle Daten, die für das nächste Bit benötigt werden; frühere Terme müssen nicht gespeichert werden.
Validierung. Das Programm prüft \(\rho^4=2\), vergleicht die schnelle Methode für \(1\le n\le200\) mit einer direkten Rekurrenz und bestätigt \(a_{13}=23321\) sowie \(a_{1003}\equiv231906014\pmod M\).
Field repräsentiert \(c_0+c_1\rho+c_2\rho^2+c_3\rho^3\). Addition und Subtraktion erfolgen komponentenweise modulo \(M\). Bei der Multiplikation entstehen zunächst sieben Polynomkoeffizienten, anschließend werden die hohen Grade mit \(\rho^4=2\) zurückgefaltet.
advance liest ein Bit des Zielindex und baut den nächsten Neun-Terme-Block mit der ungeraden oder geraden EDS-Formel. eds_value durchläuft die Binärdarstellung von \(n\) von links nach rechts, sodass nur \(O(\log n)\) Blockupdates nötig sind. sequence_value wandelt \(W_n\) zurück in \(a_n\), korrigiert das Vorzeichen und wählt je nach Parität den Skalar- oder \(\rho\)-Koeffizienten.
Der schnelle Pfad teilt nie durch einen variablen Term der ursprünglichen Folge. Die einzige Division ist die feste Division durch \(W_2=\rho\), deren Inverse \(\rho^3/2\) explizit bekannt ist. Da \(M\) ungerade ist, besitzt \(2\) die Inverse \((M+1)/2\) modulo \(M\). Die große Rechnung besteht daher nur aus Additionen, Subtraktionen, Multiplikationen und einer festen Multiplikation mit \(\rho^{-1}\).
Die direkte Rekurrenz im Testteil teilt zwar durch \(a_{k-2}\), ist aber nicht Teil des Algorithmus für den großen Index. Sie dient nur als unabhängige Kontrolle für kleine \(n\), und der Code prüft vor jeder Fermat-Inversion, dass der Nenner nicht null ist.
Nach dem Lesen eines binären Präfixes \(p\) des Zielindex ist der gespeicherte Block \(W_{p-4},W_{p-3},\ldots,W_{p+4}\); sein Zentrum ist also \(p\). Der Anfangspräfix ist das führende Bit \(1\), und der Anfangsblock ist bei \(1\) zentriert. Das Anhängen eines Bits \(b\) ändert den dargestellten Wert von \(p\) zu \(2p+b\), genau wie die Aktualisierung \(c'=2c+b\). Nach dem letzten Bit ist das Zentrum deshalb der vollständige Zielindex.
Jedes Binärbit verursacht nur konstant viele Multiplikationen und Subtraktionen in einer vierkoeffizientigen Algebra. Die Laufzeit ist daher
\[ O(\log n), \]
und der Speicherverbrauch ist \(O(1)\), weil nur zwei Neun-Terme-Blöcke und einige temporäre Elemente gehalten werden. Für \(n=10^{18}+3\) bedeutet das ungefähr sechzig Blockupdates statt \(10^{18}\) Rekurrenzschritten.
Dizi \(a_1=a_2=a_3=1\) ve \(a_4=2\) değerleriyle başlar. Sonraki terimler
\[ a_n^2=a_{n+2}a_{n-2}+u\,a_{n+1}a_{n-1} \]
bağıntısıyla belirlenir; burada \(n\) çiftse \(u=1\), tekse \(u=2\). Hedef \(n=10^{18}+3\) için \(a_n\) değerini
\[ M=1234567891 \]
modunda bulmaktır. Bu kadar büyük bir indekse sırayla ilerlemek mümkün değildir. Çözüm, katsayının tek-çift durumuna göre değişmesini \(\rho^4=2\) koşullu küçük bir cebirsel uzantıya gömer ve sonra elliptic divisibility sequence ikiye katlama formülleriyle doğrudan hedef indekse atlar.
Temel fikir, \(\rho\) adında biçimsel bir eleman tanımlayıp
\[ \rho^4=2 \]
koşulunu dayatmaktır. Hesaplar şu halkada yapılır:
\[ R=\mathbb{F}_M[\rho]/(\rho^4-2). \]
Bu halkadaki her eleman
\[ c_0+c_1\rho+c_2\rho^2+c_3\rho^3 \]
şeklinde saklanır. İki eleman çarpıldığında önce normal polinom çarpımı yapılır, sonra \(\rho^4=2\) kullanılarak derece \(4\) ve üzerindeki terimler geri katlanır. Kodda Field sınıfının yaptığı işlem tam olarak budur.
Asıl \(a_n\) dizisi yerine
\[ W_n= \begin{cases} (-1)^{(n-1)/2}a_n, & n \text{ tek},\\[2mm] (-1)^{(n-2)/2}a_n\rho, & n \text{ çift} \end{cases} \]
dizisini tanımlayalım. Başlangıç değerleri
\[ W_1=1,\qquad W_2=\rho,\qquad W_3=-1,\qquad W_4=-2\rho,\qquad W_5=-3 \]
olur. Bu dönüşümden sonra orijinal bağıntı tek bir sabit katsayılı bağıntıya dönüşür:
\[ W_{n+2}W_{n-2}=W_n^2+\rho^2W_{n+1}W_{n-1}. \]
\(n\) tek olduğunda \(W_{n+1}\) ve \(W_{n-1}\) çift indekslidir, yani ikisi de \(\rho\) katsayısı taşır. Dışarıdaki \(\rho^2\) ile birlikte \(\rho^4=2\) oluşur; bu da orijinal bağıntıdaki \(u=2\)'dir. \(n\) çift olduğunda ortak \(\rho^2\) faktörü sadeleşir ve katsayı \(u=1\) kalır. Yani tek-çift değişimi yaklaşık değil, birebir cebirsel olarak kodlanmıştır.
\(W_n\) bir elliptic divisibility sequence gibi işlenir. Kod iki standart ikiye katlama formülünü kullanır:
\[ W_{2k-1}=W_{k+1}W_{k-1}^3-W_{k-2}W_k^3, \]
ve
\[ W_{2k}=\frac{W_k}{W_2}\left(W_{k+2}W_{k-1}^2-W_{k-2}W_{k+1}^2\right). \]
Burada \(W_2=\rho\). \(\rho^4=2\) olduğundan
\[ \rho^{-1}=\frac{\rho^3}{2}. \]
Mod \(M\)'de \(2\)'ye bölmek, \((M+1)/2\) ile çarpmaktır; bu yüzden kod INV_RHO değerini \(\rho^3/2\) olarak saklar.
Algoritma her anda merkez \(c\) etrafındaki dokuz değeri tutar:
\[ W_{c-4},W_{c-3},\dots,W_c,\dots,W_{c+4}. \]
Başlangıçta \(c=1\)'dir ve pencere \(W_{-3}\)'ten \(W_5\)'e kadar gider. Hedef indeksin sıradaki ikili basamağı \(b\) ise yeni merkez
\[ c'=2c+b \]
olur. Yeni pencerenin her elemanı için indeks tekse \(W_{2k-1}\), çiftse \(W_{2k}\) formülü uygulanır. Bu formüller yalnızca \(k-2,k-1,k,k+1,k+2\) değerlerine ihtiyaç duyar; bunlar da eski dokuzlu pencerenin içindedir. Bu yüzden algoritma \(10^{18}\) adım yerine hedefin yaklaşık 60 ikili basamağı kadar güncelleme yapar.
Bitler bittiğinde pencerenin ortası \(W_n\)'dir. Eğer \(n\) tekse bu değer yalnızca sabit katsayı taşır; eğer \(n\) çiftse yalnızca \(\rho\) katsayısı taşır. Kod bunu assert eder ve sonra işaret desenini geri alır:
\[ a_n=(-1)^{(n-1)/2}W_n \quad(n\text{ tek}), \]
\[ a_n=(-1)^{(n-2)/2}[\rho]W_n \quad(n\text{ çift}). \]
Burada \([\rho]W_n\), \(W_n\)'nin \(\rho\) katsayısıdır.
Cebirsel sağlamlık. \(W_n\) tanımı orijinal tek-çift katsayısını \(\rho^4=2\) üzerinden sabit katsayılı bir bağıntıya dönüştürür. Dönüşüm geri alındığında tek indekslerde katsayı \(2\), çift indekslerde katsayı \(1\) elde edilir.
İkiye katlama sağlamlığı. odd_value ve even_value fonksiyonları EDS ikiye katlama kimlikleridir. advance her ikili basamakta aynı \(W_j\) değerlerini üretir; sadece indeksi ikili gösterim üzerinden büyütür.
Pencere tamlığı. Yeni penceredeki her değer eski penceredeki beş komşu değere bağlıdır. Dolayısıyla dokuz değerlik pencere hiçbir bilgiyi eksik bırakmaz.
Kontroller. Program \(\rho^4=2\)'yi doğrular, hızlı yöntemi \(1\le n\le200\) aralığında doğrudan hesapla karşılaştırır ve verilen \(a_{13}=23321\), \(a_{1003}\equiv231906014\pmod M\) kontrollerini çalıştırır.
Hızlı yol, orijinal bağıntıdaki gibi değişken bir \(a_{k-2}\) terimine bölme yapmaz. Tek bölme \(W_2=\rho\)'ya karşılık gelir; onun tersi de \(\rho^{-1}=\rho^3/2\) olarak sabittir. \(M\) tek olduğu için \(2\)'nin modüler tersi \((M+1)/2\)'dir. Bu nedenle hedef indeks için yapılan bütün işlemler toplama, çıkarma, çarpma ve sabit \(\rho^{-1}\) çarpımından oluşur.
Doğrudan recurrence sadece test amacıyla kullanılır. Küçük indekslerde hızlı yöntemi bağımsız bir hesapla karşılaştırır; her bölmeden önce paydanın sıfır olmadığını assert eder. Büyük hedef için bu doğrudan bölmeli yöntem kullanılmaz.
Hedef indeksin ikili gösterimi soldan sağa okunur. Okunmuş önek \(p\) ise tutulan pencere \(W_{p-4},W_{p-3},\ldots,W_{p+4}\) değerlerinden oluşur; yani pencerenin merkezi \(p\)'dir. Başlangıç öneki en soldaki \(1\) bitidir ve pencere zaten \(1\) merkezlidir. Yeni bit \(b\) eklendiğinde önek \(2p+b\) olur; advance fonksiyonunun merkezi \(2c+b\)'ye taşıması tam olarak bu yüzden doğrudur. Son bit de işlendiğinde pencerenin merkezi hedef indeksin kendisidir.
Field, dört katsayılı alan elemanını temsil eder. advance, hedef indeksin bir bitini okuyup merkezi \(2c+b\)'ye taşır ve yeni dokuzlu pencereyi kurar. eds_value, \(n\)'nin ikili gösterimini soldan sağa işleyerek \(W_n\)'ye ulaşır. sequence_value ise pariteye göre sabit katsayıyı veya \(\rho\) katsayısını seçip işareti düzeltir.
Her ikili basamakta sabit sayıda halka çarpımı ve çıkarması yapılır. Bu nedenle çalışma süresi
\[ O(\log n) \]
ve bellek kullanımı \(O(1)\)'dir. \(n=10^{18}+3\) için bu, yaklaşık altmış pencere güncellemesi anlamına gelir.
La sucesión empieza con \(a_1=a_2=a_3=1\) y \(a_4=2\). Para los términos posteriores se impone la relación
\[ a_n^2=a_{n+2}a_{n-2}+u\,a_{n+1}a_{n-1}, \]
donde \(u=1\) si \(n\) es par y \(u=2\) si \(n\) es impar. El índice objetivo es enorme, \(n=10^{18}+3\), por lo que iterar la recurrencia es imposible. Se pide calcular \(a_n\) módulo
\[ M=1234567891. \]
La implementación convierte el coeficiente alternante \(1,2,1,2,\dots\) en una elliptic divisibility sequence de coeficiente constante sobre una pequeña extensión algebraica de \(\mathbb{F}_M\), y después usa duplicación binaria para llegar directamente al índice pedido.
El obstáculo real es el coeficiente dependiente de la paridad. Introducimos un elemento formal \(\rho\) que satisface
\[ \rho^4=2. \]
Todos los cálculos se realizan módulo \(M\), dentro del álgebra de dimensión cuatro
\[ R=\mathbb{F}_M[\rho]/(\rho^4-2). \]
Un elemento de \(R\) se almacena como
\[ c_0+c_1\rho+c_2\rho^2+c_3\rho^3. \]
La multiplicación es una multiplicación polinómica ordinaria seguida de la reducción \(\rho^4=2\). Por tanto, cada término \(d_i\rho^i\) con \(i\ge4\) se pliega como \(2d_i\rho^{i-4}\). Esto es exactamente lo que implementa la multiplicación de Field.
Definimos una nueva sucesión \(W_n\) que absorbe tanto el patrón de signos como el factor de paridad:
\[ W_n= \begin{cases} (-1)^{(n-1)/2}a_n, & n \text{ impar},\\[2mm] (-1)^{(n-2)/2}a_n\rho, & n \text{ par}. \end{cases} \]
Los primeros valores pasan a ser
\[ W_1=1,\qquad W_2=\rho,\qquad W_3=-1,\qquad W_4=-2\rho,\qquad W_5=-3. \]
Con esta normalización, la recurrencia alternante original se transforma en la relación de coeficiente constante
\[ W_{n+2}W_{n-2}=W_n^2+\rho^2W_{n+1}W_{n-1}. \]
Si \(n\) es impar, los dos vecinos \(W_{n+1}\) y \(W_{n-1}\) tienen índice par y contienen un factor \(\rho\). Junto con el factor exterior \(\rho^2\) producen \(\rho^4=2\), que coincide con \(u=2\). Si \(n\) es par, el factor común \(\rho^2\) se cancela de los otros términos y queda \(u=1\). Así, la alternancia no se aproxima; queda codificada exactamente en el álgebra.
La sucesión \(W_n\) se maneja como una elliptic divisibility sequence. El código usa las identidades estándar de duplicación
\[ W_{2k-1}=W_{k+1}W_{k-1}^3-W_{k-2}W_k^3, \]
y
\[ W_{2k}=\frac{W_k}{W_2}\left(W_{k+2}W_{k-1}^2-W_{k-2}W_{k+1}^2\right). \]
Aquí \(W_2=\rho\). Como \(\rho^4=2\), se tiene
\[ \rho^{-1}=\frac{\rho^3}{2}. \]
Módulo \(M\), dividir por \(2\) equivale a multiplicar por \((M+1)/2\), por lo que el código almacena \(\rho^{-1}\) como INV_RHO. Estas identidades permiten calcular valores alrededor de \(2k\) a partir de un pequeño bloque alrededor de \(k\).
La implementación mantiene un bloque
\[ W_{c-4},W_{c-3},\ldots,W_c,\ldots,W_{c+4} \]
centrado en un índice \(c\). Al inicio \(c=1\), con la ventana
\[ 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. \]
Si el siguiente bit del índice objetivo es \(b\in\{0,1\}\), el centro pasa a
\[ c'=2c+b. \]
Para cada desplazamiento \(-4\le r\le4\), el código calcula \(W_{c'+r}\). Si \(c'+r\) es impar se aplica la fórmula de \(W_{2k-1}\); si es par se aplica la fórmula de \(W_{2k}\). Ambas fórmulas sólo necesitan valores con índices desde \(k-2\) hasta \(k+2\), que siempre caen dentro de la ventana anterior. Por eso una ventana de tamaño constante basta para un índice con unas sesenta cifras binarias.
Después de procesar todos los bits, el centro del bloque es exactamente el índice pedido \(n\), y block[4] contiene \(W_n\). El valor original se recupera deshaciendo la normalización:
\[ a_n= \begin{cases} (-1)^{(n-1)/2}W_n, & n \text{ impar},\\[2mm] (-1)^{(n-2)/2}[\rho]W_n, & n \text{ par}. \end{cases} \]
Aquí \([\rho]W_n\) denota el coeficiente de \(\rho\). La implementación comprueba que los valores de índice impar están en el campo base y que los de índice par son múltiplos puros de \(\rho\). Esas comprobaciones ayudan a detectar errores en la operación de campo o en el paso de duplicación.
Solidez algebraica. La definición de \(W_n\) transforma la recurrencia dependiente de la paridad en una sola recurrencia dentro de \(R\). Como \(\rho^4=2\), al deshacer la normalización se obtiene exactamente el coeficiente \(2\) en los pasos impares y el coeficiente \(1\) en los pasos pares.
Solidez de la duplicación. Las dos identidades usadas por odd_value y even_value son identidades de elliptic divisibility sequences. Cada llamada a advance calcula los mismos valores \(W_j\) que produciría la recurrencia, pero para el prefijo binario recién extendido.
Completitud de la ventana. Cada valor nuevo depende sólo de cinco vecinos de la ventana anterior. Por tanto, la ventana de nueve términos contiene siempre toda la información necesaria para el siguiente bit.
Validación. El programa verifica \(\rho^4=2\), compara el método rápido con una recurrencia directa para \(1\le n\le200\), y confirma los puntos de control \(a_{13}=23321\) y \(a_{1003}\equiv231906014\pmod M\).
Field representa \(c_0+c_1\rho+c_2\rho^2+c_3\rho^3\). La suma y la resta son componente a componente módulo \(M\); la multiplicación crea primero siete coeficientes polinómicos y luego reduce los grados altos usando \(\rho^4=2\).
advance lee un bit del índice objetivo. Construye la siguiente ventana de nueve términos eligiendo entre las fórmulas EDS impar y par. eds_value recorre la expansión binaria de \(n\) desde el bit más significativo al menos significativo, de modo que sólo hay \(O(\log n)\) actualizaciones de ventana. sequence_value convierte \(W_n\) de vuelta en \(a_n\), aplicando el signo y eligiendo el coeficiente adecuado según la paridad.
La ruta rápida nunca divide por un término variable de la recurrencia original. Sólo divide por \(W_2=\rho\), cuya inversa se conoce explícitamente como \(\rho^3/2\). Como \(M\) es impar, \(2\) tiene inversa \((M+1)/2\) módulo \(M\). Así, el paso de duplicación EDS sólo usa sumas, restas, multiplicaciones y una multiplicación fija por \(\rho^{-1}\).
La recurrencia directa usada para los controles sí divide por \(a_{k-2}\), pero no forma parte del algoritmo de índice grande. Sólo se usa como herramienta de validación para \(n\) pequeño, y el código comprueba que el denominador no sea cero antes de aplicar la inversión por Fermat.
Después de procesar un prefijo binario \(p\) del índice objetivo, el bloque almacenado es \(W_{p-4},W_{p-3},\ldots,W_{p+4}\); equivalentemente, su centro es \(p\). El prefijo inicial es el bit principal \(1\), así que la ventana inicial está centrada en \(1\). Al añadir un bit \(b\), el entero representado pasa de \(p\) a \(2p+b\), exactamente igual que la actualización \(c'=2c+b\). Tras el último bit, el centro coincide con el índice objetivo completo.
Cada bit requiere una cantidad constante de multiplicaciones y restas en una estructura de cuatro coeficientes. Por tanto el tiempo es
\[ O(\log n), \]
y la memoria es \(O(1)\), porque sólo se conservan dos ventanas de nueve términos y unos pocos temporales. Para \(n=10^{18}+3\), esto reemplaza \(10^{18}\) pasos de recurrencia por unas sesenta actualizaciones de ventana.
数列从 \(a_1=a_2=a_3=1\)、\(a_4=2\) 开始。之后的项满足
\[ a_n^2=a_{n+2}a_{n-2}+u\,a_{n+1}a_{n-1}, \]
其中 \(n\) 为偶数时 \(u=1\)，为奇数时 \(u=2\)。目标下标极大，\(n=10^{18}+3\)，所以不能逐项递推。要求计算 \(a_n\) 模
\[ M=1234567891. \]
实现方法把交替出现的系数 \(1,2,1,2,\dots\) 精确地嵌入到 \(\mathbb{F}_M\) 的一个小代数扩张中，使其变成一条固定形式的 elliptic divisibility sequence，然后用二进制倍增直接跳到目标下标。
困难在于系数依赖下标奇偶性。引入形式元素 \(\rho\)，令它满足
\[ \rho^4=2. \]
所有计算都在模 \(M\) 下进行，并且位于四维代数
\[ R=\mathbb{F}_M[\rho]/(\rho^4-2) \]
中。一个元素保存为
\[ c_0+c_1\rho+c_2\rho^2+c_3\rho^3. \]
乘法先按普通多项式相乘，然后用 \(\rho^4=2\) 约化。也就是说，任意 \(i\ge4\) 的项 \(d_i\rho^i\) 会被折回成 \(2d_i\rho^{i-4}\)。代码中 Field 的乘法正是这一过程。
定义新数列 \(W_n\)，同时吸收符号模式和奇偶系数：
\[ W_n= \begin{cases} (-1)^{(n-1)/2}a_n, & n \text{ 为奇数},\\[2mm] (-1)^{(n-2)/2}a_n\rho, & n \text{ 为偶数}. \end{cases} \]
最初的值变为
\[ W_1=1,\qquad W_2=\rho,\qquad W_3=-1,\qquad W_4=-2\rho,\qquad W_5=-3. \]
在这个归一化之后，原来的交替递推变成固定系数关系
\[ W_{n+2}W_{n-2}=W_n^2+\rho^2W_{n+1}W_{n-1}. \]
当 \(n\) 为奇数时，\(W_{n+1}\) 和 \(W_{n-1}\) 都是偶下标项，各自带有一个 \(\rho\)。再乘外面的 \(\rho^2\)，得到 \(\rho^4=2\)，正好对应原式中的 \(u=2\)。当 \(n\) 为偶数时，公共的 \(\rho^2\) 因子从其他项中消去，只留下 \(u=1\)。因此交替系数不是被近似处理，而是被精确编码进代数结构。
数列 \(W_n\) 作为 elliptic divisibility sequence 来处理。代码使用标准倍增恒等式
\[ W_{2k-1}=W_{k+1}W_{k-1}^3-W_{k-2}W_k^3, \]
以及
\[ W_{2k}=\frac{W_k}{W_2}\left(W_{k+2}W_{k-1}^2-W_{k-2}W_{k+1}^2\right). \]
这里 \(W_2=\rho\)。由 \(\rho^4=2\) 可得
\[ \rho^{-1}=\frac{\rho^3}{2}. \]
在模 \(M\) 下，除以 \(2\) 等价于乘以 \((M+1)/2\)，所以代码把 \(\rho^{-1}\) 存为 INV_RHO。这些恒等式允许我们从 \(k\) 附近的一小块值计算出 \(2k\) 附近的值。
实现始终维护以 \(c\) 为中心的九项窗口
\[ W_{c-4},W_{c-3},\ldots,W_c,\ldots,W_{c+4}. \]
初始中心是 \(c=1\)，窗口为
\[ 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. \]
若目标下标的下一位二进制位是 \(b\in\{0,1\}\)，新的中心为
\[ c'=2c+b. \]
对每个偏移 \(-4\le r\le4\)，代码计算 \(W_{c'+r}\)。如果 \(c'+r\) 为奇数，就使用 \(W_{2k-1}\) 公式；如果为偶数，就使用 \(W_{2k}\) 公式。两种公式都只需要 \(k-2\) 到 \(k+2\) 的值，而这些值必定位于旧的九项窗口内。因此，约六十个二进制位的目标下标只需要常数大小的窗口。
所有二进制位处理完后，窗口中心正好是目标 \(n\)，而 block[4] 中存放 \(W_n\)。原数列的值通过撤销归一化得到：
\[ a_n= \begin{cases} (-1)^{(n-1)/2}W_n, & n \text{ 为奇数},\\[2mm] (-1)^{(n-2)/2}[\rho]W_n, & n \text{ 为偶数}. \end{cases} \]
其中 \([\rho]W_n\) 表示 \(W_n\) 的 \(\rho\) 系数。代码断言奇数下标的值在基域中，偶数下标的值是纯 \(\rho\) 倍数。这些检查可以及时发现域运算或倍增步骤的错误。
代数正确性。 \(W_n\) 的定义把原来依赖奇偶性的递推变成 \(R\) 中的一条统一递推。由于 \(\rho^4=2\)，撤销归一化后，奇数步骤得到的系数正是 \(2\)，偶数步骤得到的系数正是 \(1\)。
倍增正确性。 odd_value 与 even_value 使用的是 elliptic divisibility sequence 的恒等式。每次 advance 都计算与递推相同的 \(W_j\) 值，只是下标由刚扩展的二进制前缀给出。
窗口完备性。 新窗口中的每个值只依赖旧窗口中的五个相邻值。因此九项窗口始终包含处理下一位所需的全部数据，不需要保存更早的项。
验证。 程序检查 \(\rho^4=2\)，把快速方法与 \(1\le n\le200\) 的直接递推结果比较，并验证题目给出的 \(a_{13}=23321\) 与 \(a_{1003}\equiv231906014\pmod M\)。
Field 表示 \(c_0+c_1\rho+c_2\rho^2+c_3\rho^3\)。加法和减法逐系数模 \(M\) 进行；乘法先生成七个多项式系数，再利用 \(\rho^4=2\) 折回高次项。
advance 读取目标下标的一位，使用奇数或偶数的 EDS 倍增公式构造下一组九项窗口。eds_value 从最高位到最低位扫描 \(n\) 的二进制展开，所以窗口更新次数是 \(O(\log n)\)。sequence_value 根据奇偶性选取标量系数或 \(\rho\) 系数，并恢复符号，从 \(W_n\) 得到 \(a_n\)。
快速路径从不除以原始递推中的可变项。唯一的除法是除以固定的 \(W_2=\rho\)，其逆元显式为 \(\rho^3/2\)。由于 \(M\) 是奇数，\(2\) 在模 \(M\) 下的逆元是 \((M+1)/2\)。因此大的计算只包含加、减、乘以及固定的 \(\rho^{-1}\) 乘法。
用于检查的小规模直接递推会除以 \(a_{k-2}\)，但它不是大下标算法的一部分。它只作为验证工具使用，并且代码在费马求逆前断言分母非零。
处理目标下标的二进制前缀 \(p\) 后，保存的窗口是 \(W_{p-4},W_{p-3},\ldots,W_{p+4}\)，也就是窗口中心为 \(p\)。初始前缀是最高位 \(1\)，初始窗口也以 \(1\) 为中心。追加一位 \(b\) 时，前缀从 \(p\) 变成 \(2p+b\)，这正好对应更新 \(c'=2c+b\)。最后一位处理完后，中心就是完整的目标下标。
每个二进制位只执行常数次四系数乘法和减法。因此运行时间为
\[ O(\log n), \]
空间复杂度为 \(O(1)\)，因为只需要两个九项窗口和少量临时元素。对于 \(n=10^{18}+3\)，这意味着约六十次窗口更新，而不是 \(10^{18}\) 次递推。
Последовательность начинается с \(a_1=a_2=a_3=1\) и \(a_4=2\). Для дальнейших членов задано условие
\[ a_n^2=a_{n+2}a_{n-2}+u\,a_{n+1}a_{n-1}, \]
где \(u=1\) для четного \(n\) и \(u=2\) для нечетного \(n\). Целевой индекс огромен: \(n=10^{18}+3\). Поэтому требуется найти \(a_n\) по модулю
\[ M=1234567891, \]
не выполняя линейную итерацию. Решение переводит чередующийся коэффициент \(1,2,1,2,\dots\) в рекурсию с постоянным коэффициентом над малым алгебраическим расширением \(\mathbb{F}_M\), а затем использует двоичное удвоение.
Главная трудность состоит в коэффициенте, зависящем от четности. Введем формальный элемент \(\rho\), удовлетворяющий
\[ \rho^4=2. \]
Все вычисления выполняются по модулю \(M\) в четырехмерной алгебре
\[ R=\mathbb{F}_M[\rho]/(\rho^4-2). \]
Элемент \(R\) хранится в виде
\[ c_0+c_1\rho+c_2\rho^2+c_3\rho^3. \]
Умножение сначала является обычным умножением многочленов, после чего все степени \(\rho^i\) при \(i\ge4\) сворачиваются через \(\rho^4=2\). То есть член \(d_i\rho^i\) превращается в \(2d_i\rho^{i-4}\). Именно это делает умножение в классе Field.
Определим новую последовательность \(W_n\), которая поглощает и знаковый шаблон, и фактор четности:
\[ W_n= \begin{cases} (-1)^{(n-1)/2}a_n, & n \text{ нечетно},\\[2mm] (-1)^{(n-2)/2}a_n\rho, & n \text{ четно}. \end{cases} \]
Первые значения становятся равными
\[ W_1=1,\qquad W_2=\rho,\qquad W_3=-1,\qquad W_4=-2\rho,\qquad W_5=-3. \]
После этой нормализации исходная чередующаяся рекурсия принимает вид отношения с постоянным коэффициентом:
\[ W_{n+2}W_{n-2}=W_n^2+\rho^2W_{n+1}W_{n-1}. \]
Если \(n\) нечетно, то \(W_{n+1}\) и \(W_{n-1}\) имеют четные индексы и содержат по одному фактору \(\rho\). Вместе с внешним множителем \(\rho^2\) получается \(\rho^4=2\), то есть ровно коэффициент \(u=2\). Если \(n\) четно, общий фактор \(\rho^2\) сокращается из остальных членов и остается коэффициент \(u=1\). Следовательно, чередование не приближено, а точно закодировано в алгебре.
Последовательность \(W_n\) обрабатывается как elliptic divisibility sequence. Код использует стандартные тождества удвоения
\[ W_{2k-1}=W_{k+1}W_{k-1}^3-W_{k-2}W_k^3, \]
и
\[ W_{2k}=\frac{W_k}{W_2}\left(W_{k+2}W_{k-1}^2-W_{k-2}W_{k+1}^2\right). \]
Здесь \(W_2=\rho\). Из \(\rho^4=2\) следует
\[ \rho^{-1}=\frac{\rho^3}{2}. \]
По модулю \(M\) деление на \(2\) реализуется умножением на \((M+1)/2\), поэтому код хранит \(\rho^{-1}\) как INV_RHO. Эти тождества позволяют вычислять значения около \(2k\) из небольшого блока значений около \(k\).
Реализация хранит блок
\[ W_{c-4},W_{c-3},\ldots,W_c,\ldots,W_{c+4} \]
с центром в индексе \(c\). Изначально \(c=1\), а сохраненный блок равен
\[ 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. \]
Если следующий двоичный бит целевого индекса равен \(b\in\{0,1\}\), центр переходит в
\[ c'=2c+b. \]
Для каждого смещения \(-4\le r\le4\) код вычисляет \(W_{c'+r}\). Если \(c'+r\) нечетно, применяется формула для \(W_{2k-1}\); если четно, применяется формула для \(W_{2k}\). Обе формулы требуют только значения с индексами от \(k-2\) до \(k+2\), которые гарантированно находятся в предыдущем окне. Поэтому окно постоянного размера достаточно даже для индекса с примерно шестьюдесятью двоичными битами.
После обработки всех битов центр блока точно равен требуемому \(n\), и block[4] содержит \(W_n\). Исходное значение восстанавливается обратной нормализацией:
\[ a_n= \begin{cases} (-1)^{(n-1)/2}W_n, & n \text{ нечетно},\\[2mm] (-1)^{(n-2)/2}[\rho]W_n, & n \text{ четно}. \end{cases} \]
Здесь \([\rho]W_n\) означает коэффициент при \(\rho\). Реализация проверяет, что значения с нечетными индексами лежат в базовом поле, а значения с четными индексами являются чистыми кратными \(\rho\). Эти проверки защищают от ошибок в арифметике расширения или в шаге удвоения.
Алгебраическая корректность. Определение \(W_n\) переводит исходную рекурсию с зависимостью от четности в одну рекурсию в \(R\). Благодаря \(\rho^4=2\) после обратного преобразования коэффициент равен \(2\) на нечетных шагах и \(1\) на четных.
Корректность удвоения. Две формулы, используемые в odd_value и even_value, являются тождествами elliptic divisibility sequence. Каждый вызов advance вычисляет те же \(W_j\), что и рекурсия, но для индексов, соответствующих расширенному двоичному префиксу.
Полнота окна. Каждый новый элемент зависит только от пяти соседних элементов предыдущего окна. Поэтому окно из девяти значений всегда содержит все данные, необходимые для следующего бита, и более ранние члены не нужны.
Проверка. Программа проверяет \(\rho^4=2\), сравнивает быстрый метод с прямой рекурсией для \(1\le n\le200\), а также подтверждает контрольные значения \(a_{13}=23321\) и \(a_{1003}\equiv231906014\pmod M\).
Field представляет \(c_0+c_1\rho+c_2\rho^2+c_3\rho^3\). Сложение и вычитание выполняются покомпонентно по модулю \(M\). Умножение сначала строит семь коэффициентов многочлена, а затем сворачивает старшие степени с помощью \(\rho^4=2\).
advance читает один бит целевого индекса и строит следующее окно из девяти членов, выбирая четную или нечетную EDS-формулу. eds_value просматривает двоичную запись \(n\) от старшего бита к младшему, поэтому число обновлений равно \(O(\log n)\). sequence_value переводит \(W_n\) обратно в \(a_n\), применяя знак и выбирая скалярный или \(\rho\)-коэффициент по четности.
Быстрый путь никогда не делит на переменный член исходной последовательности. Единственное деление — фиксированное деление на \(W_2=\rho\), обратный элемент которого явно равен \(\rho^3/2\). Так как \(M\) нечетно, число \(2\) имеет обратный элемент \((M+1)/2\) по модулю \(M\). Поэтому большой расчет состоит только из сложений, вычитаний, умножений и одного фиксированного умножения на \(\rho^{-1}\).
Прямая рекурсия, используемая в проверке, действительно делит на \(a_{k-2}\), но она не является частью алгоритма для большого индекса. Она служит только независимой проверкой для малых \(n\), и код проверяет, что знаменатель ненулевой, прежде чем применять инверсию по теореме Ферма.
После обработки двоичного префикса \(p\) целевого индекса сохраненный блок равен \(W_{p-4},W_{p-3},\ldots,W_{p+4}\); то есть его центр находится в \(p\). Начальный префикс — ведущий бит \(1\), и исходное окно центрировано в \(1\). Добавление бита \(b\) меняет представленное число с \(p\) на \(2p+b\), что точно совпадает с обновлением \(c'=2c+b\). После последнего бита центр равен полному целевому индексу.
На каждый двоичный бит приходится постоянное число умножений и вычитаний над четырьмя коэффициентами. Поэтому время работы равно
\[ O(\log n), \]
а память равна \(O(1)\), поскольку хранятся только два окна из девяти членов и несколько временных элементов. Для \(n=10^{18}+3\) это около шестидесяти обновлений вместо \(10^{18}\) шагов рекурсии.
تبدأ المتتالية بالقيم \(a_1=a_2=a_3=1\) و \(a_4=2\). وبعد ذلك تحقق الحدود العلاقة
\[ a_n^2=a_{n+2}a_{n-2}+u\,a_{n+1}a_{n-1}, \]
حيث \(u=1\) إذا كان \(n\) زوجياً و \(u=2\) إذا كان فردياً. الفهرس المطلوب ضخم، \(n=10^{18}+3\)، ولذلك يستحيل السير في العلاقة خطوة بخطوة. المطلوب هو حساب \(a_n\) بترديد
\[ M=1234567891. \]
تحول الخوارزمية المعامل المتناوب \(1,2,1,2,\dots\) إلى elliptic divisibility sequence ذات معامل ثابت داخل امتداد جبري صغير للحقل \(\mathbb{F}_M\)، ثم تستخدم التضعيف الثنائي للوصول مباشرة إلى الفهرس المطلوب.
العائق الأساسي هو أن المعامل يعتمد على زوجية الفهرس. ندخل عنصراً شكلياً \(\rho\) يحقق
\[ \rho^4=2 \]
ونجري جميع الحسابات بترديد \(M\) داخل الجبر رباعي الأبعاد
\[ R=\mathbb{F}_M[\rho]/(\rho^4-2). \]
يمثل كل عنصر بالصورة
\[ c_0+c_1\rho+c_2\rho^2+c_3\rho^3. \]
الضرب هو ضرب كثيرات حدود عادي أولاً، ثم اختزال باستعمال \(\rho^4=2\). لذلك فإن أي حد \(d_i\rho^i\) مع \(i\ge4\) يطوى إلى \(2d_i\rho^{i-4}\). هذا هو بالضبط ما تنفذه عملية الضرب في Field.
نعرف متتالية جديدة \(W_n\) تمتص نمط الإشارة وعامل الزوجية معاً:
\[ 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} \]
فتصبح القيم الأولى
\[ W_1=1,\qquad W_2=\rho,\qquad W_3=-1,\qquad W_4=-2\rho,\qquad W_5=-3. \]
وبهذا التطبيع تتحول العلاقة الأصلية المتناوبة إلى علاقة ثابتة الشكل:
\[ W_{n+2}W_{n-2}=W_n^2+\rho^2W_{n+1}W_{n-1}. \]
إذا كان \(n\) فردياً، فإن \(W_{n+1}\) و \(W_{n-1}\) لهما فهرسان زوجيان ويحمل كل منهما عاملاً \(\rho\). ومع العامل الخارجي \(\rho^2\) نحصل على \(\rho^4=2\)، وهو معامل \(u=2\) المطلوب. وإذا كان \(n\) زوجياً، فإن عامل \(\rho^2\) المشترك يلغى من الحدود الأخرى ويبقى \(u=1\). لذلك لا يتم تقريب التناوب، بل يرمز بدقة داخل الجبر.
تتعامل الخوارزمية مع \(W_n\) كمتتالية elliptic divisibility sequence. وتستخدم هويتي التضعيف القياسيتين
\[ W_{2k-1}=W_{k+1}W_{k-1}^3-W_{k-2}W_k^3, \]
و
\[ W_{2k}=\frac{W_k}{W_2}\left(W_{k+2}W_{k-1}^2-W_{k-2}W_{k+1}^2\right). \]
هنا \(W_2=\rho\). وبما أن \(\rho^4=2\)، فإن
\[ \rho^{-1}=\frac{\rho^3}{2}. \]
القسمة على \(2\) بترديد \(M\) هي الضرب في \((M+1)/2\)، ولذلك يخزن الكود \(\rho^{-1}\) في INV_RHO. تسمح هذه الهويات بحساب القيم حول \(2k\) انطلاقاً من كتلة صغيرة حول \(k\).
يحفظ البرنامج دائماً الكتلة
\[ W_{c-4},W_{c-3},\ldots,W_c,\ldots,W_{c+4} \]
متمركزة عند فهرس \(c\). في البداية \(c=1\)، والكتلة المخزنة هي
\[ 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. \]
عند قراءة البت التالي \(b\in\{0,1\}\) من الفهرس المطلوب، ينتقل المركز إلى
\[ c'=2c+b. \]
لكل إزاحة \(-4\le r\le4\)، يحسب الكود \(W_{c'+r}\). إذا كان الفهرس فردياً يستعمل صيغة \(W_{2k-1}\)، وإذا كان زوجياً يستعمل صيغة \(W_{2k}\). وكلتا الصيغتين تحتاج فقط إلى القيم من \(k-2\) إلى \(k+2\)، وهي موجودة دائماً في النافذة السابقة. ولهذا تكفي نافذة ثابتة الحجم حتى لفهرس له نحو ستين بتاً ثنائياً.
بعد معالجة جميع البتات، يكون مركز الكتلة هو الفهرس المطلوب \(n\)، وتحتوي block[4] على \(W_n\). نسترجع القيمة الأصلية بعكس التطبيع:
\[ 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} \]
الرمز \([\rho]W_n\) يعني معامل \(\rho\) في \(W_n\). ويتحقق الكود من أن القيم ذات الفهارس الفردية تقع في الحقل الأساسي، وأن القيم ذات الفهارس الزوجية مضاعفات صافية لـ \(\rho\). هذه الفحوص مفيدة لكشف أي خطأ في حسابات الحقل أو في خطوة التضعيف.
الصحة الجبرية. تعريف \(W_n\) يحول العلاقة الأصلية ذات المعامل المعتمد على الزوجية إلى علاقة واحدة داخل \(R\). وبسبب \(\rho^4=2\)، يظهر المعامل \(2\) تماماً في الخطوات الفردية والمعامل \(1\) تماماً في الخطوات الزوجية عند الرجوع إلى \(a_n\).
صحة التضعيف. الدالتان odd_value و even_value تستخدمان هويات elliptic divisibility sequence. لذلك فإن كل استدعاء لـ advance يحسب قيم \(W_j\) نفسها التي تعطيها العلاقة، لكن للفهرس الذي يمثله السابق الثنائي بعد إضافة بت جديد.
اكتمال النافذة. كل قيمة جديدة تعتمد فقط على خمس قيم متجاورة من النافذة السابقة. إذن النافذة ذات التسع قيم تحتوي دائماً على كل البيانات اللازمة للبت التالي، ولا حاجة إلى الحدود الأقدم.
التحقق. يتحقق البرنامج من \(\rho^4=2\)، ويقارن الطريقة السريعة مع علاقة مباشرة لكل \(1\le n\le200\)، ويتحقق من نقطتي الاختبار \(a_{13}=23321\) و \(a_{1003}\equiv231906014\pmod M\).
Field يمثل \(c_0+c_1\rho+c_2\rho^2+c_3\rho^3\). الجمع والطرح يتمان على المعاملات بترديد \(M\). أما الضرب فيبني أولاً سبعة معاملات لكثير الحدود، ثم يطوي الدرجات العالية باستعمال \(\rho^4=2\).
advance يقرأ بتاً واحداً من الفهرس المطلوب، ويبني النافذة التالية باستخدام صيغة EDS الفردية أو الزوجية. eds_value يمسح التمثيل الثنائي لـ \(n\) من البت الأعلى إلى الأدنى، ولذلك يكون عدد تحديثات النافذة \(O(\log n)\). ثم يحول sequence_value القيمة \(W_n\) إلى \(a_n\) بتصحيح الإشارة واختيار المعامل المناسب حسب الزوجية.
المسار السريع لا يقسم أبداً على حد متغير من المتتالية الأصلية. القسمة الوحيدة هي القسمة الثابتة على \(W_2=\rho\)، ومعكوسها معروف صراحةً: \(\rho^3/2\). وبما أن \(M\) فردي، فإن \(2\) له معكوس \((M+1)/2\) بترديد \(M\). لذلك تتكون خطوة التضعيف الكبيرة من جمع وطرح وضرب وضرب ثابت في \(\rho^{-1}\).
العلاقة المباشرة المستخدمة في الاختبارات تقسم فعلاً على \(a_{k-2}\)، لكنها ليست جزءاً من خوارزمية الفهرس الكبير. هي أداة تحقق للقيم الصغيرة فقط، والكود يتأكد من أن المقام غير صفري قبل استخدام عكس فيرما.
بعد معالجة سابق ثنائي \(p\) من الفهرس المطلوب، تكون الكتلة المخزنة هي \(W_{p-4},W_{p-3},\ldots,W_{p+4}\)، أي إن مركزها هو \(p\). السابق الابتدائي هو البت الرائد \(1\)، والنافذة الابتدائية متمركزة عند \(1\). عند إضافة بت \(b\)، ينتقل العدد الممثل من \(p\) إلى \(2p+b\)، وهذا يطابق تماماً تحديث المركز \(c'=2c+b\). بعد آخر بت يكون المركز هو الفهرس المطلوب كاملاً.
كل بت ثنائي ينفذ عدداً ثابتاً من عمليات الضرب والطرح على أربعة معاملات. لذلك زمن التنفيذ هو
\[ O(\log n), \]
واستخدام الذاكرة \(O(1)\)، لأن الخوارزمية تحتفظ بنافذتين من تسع قيم وبعض العناصر المؤقتة فقط. بالنسبة إلى \(n=10^{18}+3\)، فهذا يعني نحو ستين تحديثاً للنافذة بدلاً من \(10^{18}\) خطوة تكرار.