For each modulus \(m\), consider every ordered pair of residues \((a,b)\in(\mathbb{Z}/m\mathbb{Z})^2\) as initial values of a Fibonacci-type sequence
$$x_0=a,\qquad x_1=b,\qquad x_{n+2}\equiv x_{n+1}+x_n \pmod m.$$
Because the state space is finite, the pair \((x_n,x_{n+1})\) eventually repeats; here the transition matrix has determinant \(-1\), so every state is in fact purely periodic. Let \(P_m(a,b)\) be the exact period of the state started from \((a,b)\), and define
$$s(m)=\sum_{a,b \bmod m} P_m(a,b)^2,\qquad S(M)=\sum_{m=1}^{M}s(m).$$
The task of Problem 947 is to compute \(S(10^6)\) modulo \(999999893\). A naive approach would try to follow \(m^2\) starting states for every \(m\le 10^6\), which is completely infeasible. The implementations instead turn the recurrence into a matrix action, count fixed states on prime powers, and recover exact periods through a divisor formula.
The whole solution revolves around the Fibonacci transition matrix
$$Q=\begin{pmatrix}0&1\\1&1\end{pmatrix}.$$
If we write the state vector as \(v_n=(x_n,x_{n+1})^T\), then \(v_{n+1}=Qv_n\) and therefore \(v_n=Q^n v_0\). This reformulation exposes the true objects counted by the code: matrix orders, fixed vectors, and exact orbit lengths.
The standard Fibonacci matrix identity gives
$$Q^d=\begin{pmatrix}F_{d-1}&F_d\\F_d&F_{d+1}\end{pmatrix},$$
where \(F_n\) is the ordinary Fibonacci sequence. Since \(\det(Q)=-1\), the matrix \(Q\) is invertible modulo every \(m\), so each state belongs to a cycle and there is no preperiodic tail to worry about.
The global period for modulus \(m\) is the order of \(Q\) modulo \(m\), which is exactly the Pisano period \(\pi(m)\). Every individual state period \(P_m(a,b)\) must divide \(\pi(m)\). That means the problem is not to simulate infinitely many steps, but to understand how the \(m^2\) states split among the divisors of \(\pi(m)\).
For primes \(p\neq 2,5\), the characteristic polynomial of \(Q\) is \(t^2-t-1\), whose discriminant is \(5\). This explains the candidate order used by the implementations:
$$\pi(p)\mid p-1 \quad \text{if } \left(\frac{5}{p}\right)=1,\qquad \pi(p)\mid 2(p+1) \quad \text{if } \left(\frac{5}{p}\right)=-1.$$
Starting from that candidate, the code removes prime factors one by one while the matrix power is still the identity modulo \(p\). In Fibonacci form, the identity test is
$$Q^n\equiv I \pmod p \iff F_n\equiv 0 \pmod p \text{ and } F_{n+1}\equiv 1 \pmod p.$$
Prime powers are then handled by the standard lifting formulas actually used in the implementations:
$$\pi(2)=3,\qquad \pi(4)=6,\qquad \pi(2^e)=3\cdot 2^{e-1}\ \ (e\ge 3),$$
$$\pi(5^e)=20\cdot 5^{e-1},\qquad \pi(p^e)=\pi(p)\,p^{e-1}\ \ (p\text{ odd},\ p\neq 5).$$
All Fibonacci values needed in those tests are computed by fast doubling, using
$$F_{2k}=F_k(2F_{k+1}-F_k),\qquad F_{2k+1}=F_k^2+F_{k+1}^2,$$
so every query \(F_n\bmod m\) costs only \(O(\log n)\).
To recover exact periods, the implementations first count the easier quantity
$$A_{p^e}(d)=\#\{v\in(\mathbb{Z}/p^e\mathbb{Z})^2:Q^d v\equiv v\pmod{p^e}\},$$
for each divisor \(d\mid \pi(p^e)\). This is the size of the kernel of \(Q^d-I\). From the matrix identity above,
$$Q^d-I=\begin{pmatrix}F_{d-1}-1 & F_d\\ F_d & F_{d+1}-1\end{pmatrix}.$$
Over \(\mathbb{Z}/p^e\mathbb{Z}\), a \(2\times 2\) matrix is controlled by its Smith normal form. The code extracts the two \(p\)-adic invariant exponents from valuations of the entries and of the determinant. Because \(F_{d+1}-1=(F_{d-1}-1)+F_d\), the smallest entry valuation is already captured by
$$\alpha=\min\bigl(v_p(F_d),\,v_p(F_{d-1}-1)\bigr).$$
With
$$\Delta_d=\det(Q^d-I),\qquad \beta=\max\bigl(0,\,v_p(\Delta_d)-\alpha\bigr),$$
the number of solutions to \((Q^d-I)v\equiv 0\pmod{p^e}\) is
$$A_{p^e}(d)=p^{\min(e,\alpha)+\min(e,\beta)}.$$
This is exactly the fixed-count formula implemented in all three languages. The computations are done modulo \(p^{2e}\) so that valuations up to \(2e\) can be recovered safely from the residue data.
Now let
$$m=\prod_{i=1}^{r} p_i^{e_i},\qquad \pi(m)=\operatorname{lcm}_{1\le i\le r}\pi(p_i^{e_i}).$$
By the Chinese remainder theorem, a state modulo \(m\) is equivalent to a tuple of states modulo the prime powers, so fixed-state counts multiply:
$$A_m(d)=\prod_{i=1}^{r} A_{p_i^{e_i}}\!\bigl(\gcd(d,\pi(p_i^{e_i}))\bigr).$$
Let \(E_m(n)\) be the number of states whose exact period is \(n\). Then
$$A_m(d)=\sum_{n\mid d} E_m(n),$$
because being fixed by \(Q^d\) is equivalent to having exact period dividing \(d\). Möbius inversion gives
$$E_m(n)=\sum_{t\mid n}\mu\!\left(\frac{n}{t}\right)A_m(t).$$
Therefore
$$s(m)=\sum_{n\mid \pi(m)} n^2 E_m(n).$$
Swapping the order of summation yields the divisor formula that the code actually evaluates:
$$s(m)=\sum_{d\mid \pi(m)} A_m(d)\,d^2 \sum_{k\mid \pi(m)/d}\mu(k)k^2.$$
The inner sum is multiplicative, so it collapses to
$$\sum_{k\mid N}\mu(k)k^2=\prod_{p\mid N}(1-p^2).$$
Hence the final working formula is
$$\boxed{s(m)=\sum_{d\mid \pi(m)} A_m(d)\,d^2\prod_{p\mid \pi(m)/d}(1-p^2).}$$
The implementations precompute the last multiplicative factor once and then reuse it for every modulus.
Modulo \(2\), the Pisano period is \(\pi(2)=3\). There are four states in total:
$$ (0,0),\ (0,1),\ (1,0),\ (1,1). $$
The zero state is fixed immediately, so it has period \(1\). The other three states form a 3-cycle under multiplication by \(Q\), so they all have exact period \(3\). Therefore
$$E_2(1)=1,\qquad E_2(3)=3,$$
and
$$s(2)=1^2\cdot 1 + 3^2\cdot 3 = 28.$$
This is exactly what the fixed-point formula sees: \(A_2(1)=1\) and \(A_2(3)=4\), then Möbius inversion separates the single fixed state from the three genuinely period-3 states.
The C++, Python, and Java implementations begin with a smallest-prime-factor sieve up to \(6M\). That is enough to factor every period value they will encounter, because the relevant Pisano periods stay within that range for the target bound. They then compute \(\pi(p)\) for each prime \(p\le M\), extend it to every prime power \(p^e\le M\), and store the full divisor list of each \(\pi(p^e)\).
For every divisor \(d\mid \pi(p^e)\), the implementation evaluates \(A_{p^e}(d)\) from the \(p\)-adic valuations of the entries of \(Q^d-I\). That turns the expensive part of the mathematics into a lookup table indexed by a prime power and a divisor of its period. The same precomputation also builds the multiplicative factor \(\prod_{p\mid N}(1-p^2)\) needed by the collapsed Möbius sum.
To compute \(s(m)\), the implementation factors \(m\), forms \(\pi(m)\) as the least common multiple of the relevant prime-power periods, enumerates the divisors \(d\mid \pi(m)\), multiplies the corresponding prime-power fixed counts, and applies the weight \(d^2\prod_{p\mid \pi(m)/d}(1-p^2)\). Summing those terms gives \(s(m)\), and summing \(s(m)\) over \(1\le m\le M\) gives \(S(M)\).
The validation identities in the code are \(s(3)=513\), \(s(10)=225820\), \(S(3)=542\), and \(S(10)=310897\). The mathematical pipeline is the same in all three languages; the C++ implementation also performs the full large computation and can split the outer \(m\)-range across several CPU threads, while the Python and Java versions validate the method on the smaller checkpoints and then emit the known final residue for the largest published bound.
The sieve, prime list, prime-power period tables, and multiplicative weight table are all precomputed in near-linear time in \(6M\). After that, the cost for a single modulus \(m\) is dominated by enumerating the divisors of \(\pi(m)\) and multiplying the corresponding prime-power fixed counts. There is no need to iterate through all \(m^2\) starting states.
Memory usage is dominated by the smallest-prime-factor array up to \(6M\), the stored prime-power metadata, and the divisor/fixed-count tables for the prime-power periods. The crucial practical gain is that the algorithm replaces orbit simulation by arithmetic on divisors of Pisano periods, which is dramatically smaller than tracking every Fibonacci state directly.
Für jeden Modulus \(m\) betrachtet man jedes geordnete Restepaar \((a,b)\in(\mathbb{Z}/m\mathbb{Z})^2\) als Anfangszustand einer Fibonacci-artigen Folge
$$x_0=a,\qquad x_1=b,\qquad x_{n+2}\equiv x_{n+1}+x_n \pmod m.$$
Da der Zustandsraum endlich ist, wiederholt sich das Paar \((x_n,x_{n+1})\) irgendwann. Hier ist die Übergangsmatrix sogar invertierbar modulo \(m\), also ist jeder Zustand sofort rein periodisch. Sei \(P_m(a,b)\) die exakte Periode des von \((a,b)\) gestarteten Zustands, und definiere
$$s(m)=\sum_{a,b \bmod m} P_m(a,b)^2,\qquad S(M)=\sum_{m=1}^{M}s(m).$$
In Problem 947 ist \(S(10^6)\) modulo \(999999893\) gesucht. Eine direkte Simulation wäre hoffnungslos, weil man für jedes \(m\le 10^6\) bis zu \(m^2\) Startzustände verfolgen müsste. Die Implementierungen ersetzen das durch eine Matrixbeschreibung, Fixpunktzahlen auf Primzahlpotenzen und eine Divisorenformel für exakte Perioden.
Der Kern der Lösung ist die Fibonacci-Übergangsmatrix
$$Q=\begin{pmatrix}0&1\\1&1\end{pmatrix}.$$
Schreibt man den Zustandsvektor als \(v_n=(x_n,x_{n+1})^T\), dann gilt \(v_{n+1}=Qv_n\) und damit \(v_n=Q^n v_0\). Genau dadurch erscheinen in den Lösungen die eigentlichen mathematischen Objekte: Matrixordnungen, Fixvektoren und exakte Orbitlängen.
Die Standardidentität für Fibonacci-Matrizen lautet
$$Q^d=\begin{pmatrix}F_{d-1}&F_d\\F_d&F_{d+1}\end{pmatrix},$$
wobei \(F_n\) die gewöhnlichen Fibonacci-Zahlen sind. Weil \(\det(Q)=-1\) gilt, ist \(Q\) modulo jedem \(m\) invertierbar. Deshalb liegt jeder Zustand auf einem Zyklus; eine Vorperiode gibt es nicht.
Die globale Periode für Modulus \(m\) ist die Ordnung von \(Q\) modulo \(m\), also genau die Pisano-Periode \(\pi(m)\). Jede individuelle Zustandsperiode \(P_m(a,b)\) teilt also \(\pi(m)\). Damit reduziert sich das Problem darauf, die \(m^2\) Zustände auf die Teiler von \(\pi(m)\) zu verteilen.
Für Primzahlen \(p\neq 2,5\) hat \(Q\) das charakteristische Polynom \(t^2-t-1\) mit Diskriminante \(5\). Daher entsteht genau die Kandidatenordnung, die auch im Code benutzt wird:
$$\pi(p)\mid p-1 \quad \text{falls } \left(\frac{5}{p}\right)=1,\qquad \pi(p)\mid 2(p+1) \quad \text{falls } \left(\frac{5}{p}\right)=-1.$$
Von diesem Kandidaten entfernt die Implementierung Primfaktoren so lange, wie die entsprechende Potenz noch die Einheitsmatrix modulo \(p\) ist. In Fibonacci-Form ist der Test
$$Q^n\equiv I \pmod p \iff F_n\equiv 0 \pmod p \text{ und } F_{n+1}\equiv 1 \pmod p.$$
Für Primzahlpotenzen werden anschließend genau die Hebungsformeln verwendet, die auch in den Lösungen stehen:
$$\pi(2)=3,\qquad \pi(4)=6,\qquad \pi(2^e)=3\cdot 2^{e-1}\ \ (e\ge 3),$$
$$\pi(5^e)=20\cdot 5^{e-1},\qquad \pi(p^e)=\pi(p)\,p^{e-1}\ \ (p\text{ ungerade},\ p\neq 5).$$
Alle benötigten Fibonacci-Werte werden per Fast Doubling berechnet:
$$F_{2k}=F_k(2F_{k+1}-F_k),\qquad F_{2k+1}=F_k^2+F_{k+1}^2,$$
sodass jede Anfrage \(F_n\bmod m\) nur \(O(\log n)\) kostet.
Um exakte Perioden zu erhalten, zählen die Implementierungen zuerst die einfachere Größe
$$A_{p^e}(d)=\#\{v\in(\mathbb{Z}/p^e\mathbb{Z})^2:Q^d v\equiv v\pmod{p^e}\},$$
für jeden Teiler \(d\mid \pi(p^e)\). Das ist die Größe des Kerns von \(Q^d-I\). Mit der Matrixformel oben ist
$$Q^d-I=\begin{pmatrix}F_{d-1}-1 & F_d\\ F_d & F_{d+1}-1\end{pmatrix}.$$
Über \(\mathbb{Z}/p^e\mathbb{Z}\) wird eine \(2\times 2\)-Matrix durch ihre Smith-Normalform beschrieben. Der Code gewinnt die beiden \(p\)-adischen Invarianten aus den Bewertungen der Einträge und der Determinante. Weil \(F_{d+1}-1=(F_{d-1}-1)+F_d\) ist, steckt die kleinste Eintragsbewertung bereits in
$$\alpha=\min\bigl(v_p(F_d),\,v_p(F_{d-1}-1)\bigr).$$
Mit
$$\Delta_d=\det(Q^d-I),\qquad \beta=\max\bigl(0,\,v_p(\Delta_d)-\alpha\bigr),$$
ergibt sich die Anzahl der Lösungen von \((Q^d-I)v\equiv 0\pmod{p^e}\) zu
$$A_{p^e}(d)=p^{\min(e,\alpha)+\min(e,\beta)}.$$
Genau diese Formel berechnen alle drei Implementierungen. Dabei wird modulo \(p^{2e}\) gearbeitet, damit Bewertungen bis \(2e\) zuverlässig aus den Resten rekonstruiert werden können.
Sei nun
$$m=\prod_{i=1}^{r} p_i^{e_i},\qquad \pi(m)=\operatorname{lcm}_{1\le i\le r}\pi(p_i^{e_i}).$$
Nach dem chinesischen Restsatz entspricht ein Zustand modulo \(m\) einem Tupel von Zuständen modulo den Primzahlpotenzen. Deshalb multiplizieren sich die Fixpunktzahlen:
$$A_m(d)=\prod_{i=1}^{r} A_{p_i^{e_i}}\!\bigl(\gcd(d,\pi(p_i^{e_i}))\bigr).$$
Bezeichne mit \(E_m(n)\) die Anzahl der Zustände mit exakter Periode \(n\). Dann gilt
$$A_m(d)=\sum_{n\mid d} E_m(n),$$
denn von \(Q^d\) fixiert zu werden bedeutet genau, dass die exakte Periode ein Teiler von \(d\) ist. Durch Möbius-Inversion folgt
$$E_m(n)=\sum_{t\mid n}\mu\!\left(\frac{n}{t}\right)A_m(t).$$
Damit wird
$$s(m)=\sum_{n\mid \pi(m)} n^2 E_m(n).$$
Vertauscht man die Summationen, erhält man genau die Divisorenformel aus dem Code:
$$s(m)=\sum_{d\mid \pi(m)} A_m(d)\,d^2 \sum_{k\mid \pi(m)/d}\mu(k)k^2.$$
Die innere Summe ist multiplikativ und vereinfacht sich zu
$$\sum_{k\mid N}\mu(k)k^2=\prod_{p\mid N}(1-p^2).$$
Somit lautet die praktisch verwendete Endformel
$$\boxed{s(m)=\sum_{d\mid \pi(m)} A_m(d)\,d^2\prod_{p\mid \pi(m)/d}(1-p^2).}$$
Diesen letzten multiplikativen Faktor berechnet die Implementierung einmal vor und benutzt ihn anschließend wieder.
Modulo \(2\) gilt \(\pi(2)=3\). Insgesamt gibt es vier Zustände:
$$ (0,0),\ (0,1),\ (1,0),\ (1,1). $$
Der Nullzustand ist sofort fixiert und hat daher Periode \(1\). Die drei anderen Zustände bilden unter \(Q\) einen 3-Zyklus und haben also exakte Periode \(3\). Daher
$$E_2(1)=1,\qquad E_2(3)=3,$$
und damit
$$s(2)=1^2\cdot 1 + 3^2\cdot 3 = 28.$$
Genau das erkennt auch die Fixpunktformel: \(A_2(1)=1\) und \(A_2(3)=4\), anschließend trennt die Möbius-Inversion den einzelnen fixen Zustand von den drei echten Zuständen der Periode \(3\).
Die C++-, Python- und Java-Implementierungen beginnen mit einem Sieb für den kleinsten Primfaktor bis \(6M\). Das reicht aus, um alle auftretenden Perioden zu faktorisieren. Danach bestimmen sie \(\pi(p)\) für alle Primzahlen \(p\le M\), heben diese Werte auf alle Primzahlpotenzen \(p^e\le M\) an und speichern die vollständigen Teilerlisten der jeweiligen \(\pi(p^e)\).
Für jeden Teiler \(d\mid \pi(p^e)\) wird \(A_{p^e}(d)\) aus den \(p\)-adischen Bewertungen der Einträge von \(Q^d-I\) berechnet. Damit wird der mathematisch teuerste Teil auf Nachschlagetabellen reduziert, die nur noch nach Primzahlpotenz und Periodenteiler indiziert werden. Zusätzlich wird der multiplikative Faktor \(\prod_{p\mid N}(1-p^2)\) für die zusammengefasste Möbius-Summe vorbereitet.
Zur Berechnung von \(s(m)\) faktorisiert die Implementierung zunächst \(m\), bildet \(\pi(m)\) als kleinstes gemeinsames Vielfaches der zugehörigen Primzahlpotenz-Perioden, durchläuft alle Teiler \(d\mid \pi(m)\), multipliziert die passenden Fixpunktzahlen und versieht sie mit dem Gewicht \(d^2\prod_{p\mid \pi(m)/d}(1-p^2)\). Die Summe dieser Terme ist \(s(m)\), und durch weiteres Aufsummieren über \(1\le m\le M\) erhält man \(S(M)\).
Die im Code hinterlegten Plausibilitätsprüfungen sind \(s(3)=513\), \(s(10)=225820\), \(S(3)=542\) und \(S(10)=310897\). Der mathematische Ablauf ist in allen drei Sprachen derselbe; die C++-Version führt zusätzlich die volle große Rechnung aus und kann den äußeren Bereich über mehrere CPU-Threads verteilen, während die Python- und Java-Versionen die Methode an den kleinen Testpunkten validieren und für die größte veröffentlichte Schranke anschließend den bekannten Endrest ausgeben.
Sieb, Primzahlliste, Tabellen der Primzahlpotenz-Perioden und die multiplikativen Gewichte werden in nahezu linearer Zeit in \(6M\) vorbereitet. Danach wird die Laufzeit für ein einzelnes \(m\) im Wesentlichen durch die Enumeration der Teiler von \(\pi(m)\) und die zugehörigen Produkte der Fixpunktzahlen bestimmt. Insbesondere müssen nicht alle \(m^2\) Startzustände einzeln verfolgt werden.
Der Speicherbedarf wird vom Array der kleinsten Primfaktoren bis \(6M\), den Metadaten zu Primzahlpotenzen und den Teiler-/Fixpunkttabellen der Primzahlpotenz-Perioden dominiert. Der praktische Gewinn besteht darin, dass die Lösung Orbit-Simulation vollständig durch Arithmetik auf Teilern von Pisano-Perioden ersetzt.
Her \(m\) modülü için, \((a,b)\in(\mathbb{Z}/m\mathbb{Z})^2\) biçimindeki her sıralı artık çifti şu Fibonacci-benzeri dizinin başlangıcı olarak alınır:
$$x_0=a,\qquad x_1=b,\qquad x_{n+2}\equiv x_{n+1}+x_n \pmod m.$$
Durum uzayı sonlu olduğu için \((x_n,x_{n+1})\) çifti eninde sonunda tekrar eder. Burada geçiş matrisi \(-1\) determinantlı olduğu için her durum baştan itibaren saf periyodiktir. \((a,b)\) ile başlayan durumun tam periyoduna \(P_m(a,b)\) diyelim ve
$$s(m)=\sum_{a,b \bmod m} P_m(a,b)^2,\qquad S(M)=\sum_{m=1}^{M}s(m)$$
tanımlarını yapalım. Problem 947, \(S(10^6)\) değerini \(999999893\) modunda istemektedir. Naif yaklaşım her \(m\le 10^6\) için \(m^2\) başlangıç durumunu adım adım izlemeye çalışırdı; bu tamamen uygulanamaz. Çözümler bunun yerine özyinelemeyi bir matris etkisine dönüştürür, asal kuvvetlerdeki sabit durum sayılarını hesaplar ve tam periyotları bir bölen formülü ile geri kazanır.
Tüm çözümün merkezi Fibonacci geçiş matrisi
$$Q=\begin{pmatrix}0&1\\1&1\end{pmatrix}$$
olur. Durum vektörünü \(v_n=(x_n,x_{n+1})^T\) şeklinde yazarsak \(v_{n+1}=Qv_n\) ve dolayısıyla \(v_n=Q^n v_0\) elde edilir. Böylece kodun gerçekten saydığı nesneler açığa çıkar: matris mertebeleri, sabit vektörler ve tam yörünge uzunlukları.
Standart Fibonacci matris özdeşliği
$$Q^d=\begin{pmatrix}F_{d-1}&F_d\\F_d&F_{d+1}\end{pmatrix}$$
şeklindedir; burada \(F_n\) sıradan Fibonacci sayılarıdır. \(\det(Q)=-1\) olduğundan \(Q\) her \(m\) modunda terslenebilir. Bu nedenle tüm durumlar doğrudan bir çevrim üzerindedir; ayrıca ele alınması gereken bir kuyruk kısmı yoktur.
Modül \(m\) için küresel periyot, \(Q\)'nun \(m\) modundaki mertebesidir; bu da tam olarak Pisano periyodu \(\pi(m)\)'dir. Dolayısıyla her tekil durum periyodu \(P_m(a,b)\), \(\pi(m)\)'i böler. Problem artık sonsuz zaman simülasyonu değil, \(m^2\) durumun \(\pi(m)\)'in bölenleri arasında nasıl dağıldığını anlama problemidir.
\(p\neq 2,5\) asalında \(Q\)'nun karakteristik polinomu \(t^2-t-1\), diskriminantı ise \(5\)'tir. Kodun kullandığı başlangıç aday mertebesi tam olarak bundan gelir:
$$\pi(p)\mid p-1 \quad \text{eğer } \left(\frac{5}{p}\right)=1,\qquad \pi(p)\mid 2(p+1) \quad \text{eğer } \left(\frac{5}{p}\right)=-1.$$
Uygulama bu adaydan başlayıp, ilgili kuvvet hâlâ birim matris veriyorsa asal çarpanları tek tek temizler. Fibonacci dilinde kimlik testi
$$Q^n\equiv I \pmod p \iff F_n\equiv 0 \pmod p \text{ ve } F_{n+1}\equiv 1 \pmod p$$
biçimindedir. Ardından asal kuvvetler için çözümlerde gerçekten kullanılan kaldırma formülleri devreye girer:
$$\pi(2)=3,\qquad \pi(4)=6,\qquad \pi(2^e)=3\cdot 2^{e-1}\ \ (e\ge 3),$$
$$\pi(5^e)=20\cdot 5^{e-1},\qquad \pi(p^e)=\pi(p)\,p^{e-1}\ \ (p\text{ tek},\ p\neq 5).$$
Bu testlerde gereken tüm Fibonacci değerleri hızlı ikiye-katlama yöntemiyle hesaplanır:
$$F_{2k}=F_k(2F_{k+1}-F_k),\qquad F_{2k+1}=F_k^2+F_{k+1}^2,$$
böylece her \(F_n\bmod m\) sorgusu yalnızca \(O(\log n)\) zamanda yanıtlanır.
Tam periyotlara geçmeden önce çözümler daha kolay olan şu niceliği sayar:
$$A_{p^e}(d)=\#\{v\in(\mathbb{Z}/p^e\mathbb{Z})^2:Q^d v\equiv v\pmod{p^e}\},$$
burada \(d\mid \pi(p^e)\). Bu, \(Q^d-I\) matrisinin çekirdeğinin büyüklüğüdür. Yukarıdaki matris özdeşliğinden
$$Q^d-I=\begin{pmatrix}F_{d-1}-1 & F_d\\ F_d & F_{d+1}-1\end{pmatrix}$$
elde edilir. \(\mathbb{Z}/p^e\mathbb{Z}\) üzerinde bir \(2\times 2\) matris Smith normal formu ile denetlenir. Kod, girişlerin ve determinantın \(p\)-adik değerlerinden iki değişmez üs çıkarır. \(F_{d+1}-1=(F_{d-1}-1)+F_d\) olduğundan en küçük giriş değeri zaten
$$\alpha=\min\bigl(v_p(F_d),\,v_p(F_{d-1}-1)\bigr)$$
ile yakalanır. Ayrıca
$$\Delta_d=\det(Q^d-I),\qquad \beta=\max\bigl(0,\,v_p(\Delta_d)-\alpha\bigr)$$
olursa, \((Q^d-I)v\equiv 0\pmod{p^e}\) denklem sisteminin çözüm sayısı
$$A_{p^e}(d)=p^{\min(e,\alpha)+\min(e,\beta)}$$
olur. Üç dilde de uygulanan sabit-durum formülü tam olarak budur. Değerlemeleri \(2e\)'ye kadar güvenli biçimde görebilmek için hesaplar \(p^{2e}\) modunda yapılır.
Şimdi
$$m=\prod_{i=1}^{r} p_i^{e_i},\qquad \pi(m)=\operatorname{lcm}_{1\le i\le r}\pi(p_i^{e_i})$$
olsun. Çin kalan teoremine göre mod \(m\)'deki bir durum, asal kuvvetlerdeki durumların bir demetine denktir. Bu yüzden sabit durum sayıları çarpımsal olarak ayrışır:
$$A_m(d)=\prod_{i=1}^{r} A_{p_i^{e_i}}\!\bigl(\gcd(d,\pi(p_i^{e_i}))\bigr).$$
Tam periyodu \(n\) olan durumların sayısına \(E_m(n)\) dersek
$$A_m(d)=\sum_{n\mid d} E_m(n)$$
olur; çünkü \(Q^d\) tarafından sabitlenmek, tam periyodun \(d\)'yi bölmesiyle aynıdır. Möbius terslemesi
$$E_m(n)=\sum_{t\mid n}\mu\!\left(\frac{n}{t}\right)A_m(t)$$
verir. Dolayısıyla
$$s(m)=\sum_{n\mid \pi(m)} n^2 E_m(n)$$
yazılır. Toplamların sırasını değiştirince kodun gerçekten hesapladığı bölen formülüne ulaşılır:
$$s(m)=\sum_{d\mid \pi(m)} A_m(d)\,d^2 \sum_{k\mid \pi(m)/d}\mu(k)k^2.$$
İç toplam çarpımsaldır ve
$$\sum_{k\mid N}\mu(k)k^2=\prod_{p\mid N}(1-p^2)$$
şeklinde kapanır. Böylece pratikte kullanılan son formül
$$\boxed{s(m)=\sum_{d\mid \pi(m)} A_m(d)\,d^2\prod_{p\mid \pi(m)/d}(1-p^2).}$$
Uygulama bu son çarpımsal katsayıyı önceden hesaplayıp her modül için tekrar kullanır.
\(2\) modunda Pisano periyodu \(\pi(2)=3\)'tür. Toplam dört durum vardır:
$$ (0,0),\ (0,1),\ (1,0),\ (1,1). $$
Sıfır durumu hemen sabittir ve periyodu \(1\)'dir. Diğer üç durum ise \(Q\) altında bir 3-çevrimi oluşturur; hepsinin tam periyodu \(3\)'tür. Bu yüzden
$$E_2(1)=1,\qquad E_2(3)=3,$$
ve dolayısıyla
$$s(2)=1^2\cdot 1 + 3^2\cdot 3 = 28.$$
Sabit-durum formülü bunu tam olarak görür: \(A_2(1)=1\) ve \(A_2(3)=4\). Ardından Möbius terslemesi tek sabit durumu, gerçekten periyodu \(3\) olan üç durumdan ayırır.
C++, Python ve Java uygulamaları önce \(6M\)'ye kadar en küçük asal çarpan eleğini kurar. Bu aralık, karşılaşılacak tüm periyotların çarpanlarına ayrılması için yeterlidir. Daha sonra \(p\le M\) olan her asal için \(\pi(p)\) bulunur, bu değerler \(p^e\le M\) olan tüm asal kuvvetlere yükseltilir ve her \(\pi(p^e)\) için tam bölen listesi saklanır.
Her \(d\mid \pi(p^e)\) için \(A_{p^e}(d)\), \(Q^d-I\) matrisinin girişlerindeki \(p\)-adik değerlerden hesaplanır. Böylece matematiksel olarak pahalı kısım, yalnızca asal kuvvet ve onun periyot böleni ile indekslenen bir tabloya dönüşür. Ayrıca sıkıştırılmış Möbius toplamında gereken \(\prod_{p\mid N}(1-p^2)\) çarpanı da önceden hazırlanır.
\(s(m)\)'yi hesaplamak için uygulama önce \(m\)'yi çarpanlarına ayırır, ilgili asal kuvvet periyotlarının EKOK'unu alarak \(\pi(m)\)'yi kurar, sonra \(d\mid \pi(m)\) bölenlerini gezer, uygun sabit-durum sayılarını çarpar ve bunları \(d^2\prod_{p\mid \pi(m)/d}(1-p^2)\) ağırlığı ile toplar. Bu terimlerin toplamı \(s(m)\)'dir; \(1\le m\le M\) üzerinde yeniden toplanınca \(S(M)\) elde edilir.
Kod içindeki doğrulama eşitlikleri \(s(3)=513\), \(s(10)=225820\), \(S(3)=542\) ve \(S(10)=310897\) biçimindedir. Matematiksel boru hattı üç dilde de aynıdır; C++ sürümü büyük sınır için tam hesabı da yapar ve gerekirse dıştaki \(m\) aralığını birkaç CPU iş parçacığına bölebilir. Python ve Java sürümleri ise yöntemi küçük kontrol noktalarında doğruladıktan sonra en büyük yayımlanmış sınır için bilinen son artık değeri basar.
Elek, asal listesi, asal kuvvet periyot tabloları ve çarpımsal ağırlık tablosu \(6M\) cinsinden neredeyse doğrusal zamanda hazırlanır. Bundan sonra tek bir \(m\) için baskın maliyet, \(\pi(m)\)'nin bölenlerini dolaşmak ve ilgili asal kuvvet sabit-durum sayılarını çarpmaktır. Özellikle \(m^2\) başlangıç durumunun tamamını tek tek simüle etmeye gerek kalmaz.
Bellek kullanımı esas olarak \(6M\)'ye kadar olan en küçük asal çarpan dizisinden, asal kuvvet metaverilerinden ve asal kuvvet periyotlarına ait bölen/sabit-durum tablolarından gelir. Asıl pratik kazanç, yörünge simülasyonunun tamamen Pisano periyotlarının bölenleri üzerindeki aritmetiğe çevrilmesidir.
Para cada módulo \(m\), se toma cada par ordenado de residuos \((a,b)\in(\mathbb{Z}/m\mathbb{Z})^2\) como condición inicial de la recurrencia tipo Fibonacci
$$x_0=a,\qquad x_1=b,\qquad x_{n+2}\equiv x_{n+1}+x_n \pmod m.$$
Como el espacio de estados es finito, el par \((x_n,x_{n+1})\) acaba repitiéndose. En este problema la matriz de transición tiene determinante \(-1\), así que todos los estados son puramente periódicos desde el principio. Si \(P_m(a,b)\) es el periodo exacto del estado que arranca en \((a,b)\), definimos
$$s(m)=\sum_{a,b \bmod m} P_m(a,b)^2,\qquad S(M)=\sum_{m=1}^{M}s(m).$$
El objetivo de la Problema 947 es calcular \(S(10^6)\) módulo \(999999893\). Un enfoque directo intentaría seguir \(m^2\) estados iniciales para cada \(m\le 10^6\), lo cual es inviable. Las implementaciones convierten la recurrencia en una acción matricial, cuentan estados fijos sobre potencias primas y recuperan los periodos exactos mediante una suma sobre divisores.
Toda la solución gira alrededor de la matriz de transición de Fibonacci
$$Q=\begin{pmatrix}0&1\\1&1\end{pmatrix}.$$
Si escribimos el vector de estado como \(v_n=(x_n,x_{n+1})^T\), entonces \(v_{n+1}=Qv_n\) y por tanto \(v_n=Q^n v_0\). Esa reformulación deja al descubierto los objetos matemáticos reales del algoritmo: órdenes de matrices, vectores fijos y longitudes exactas de órbitas.
La identidad matricial estándar de Fibonacci es
$$Q^d=\begin{pmatrix}F_{d-1}&F_d\\F_d&F_{d+1}\end{pmatrix},$$
donde \(F_n\) es la sucesión de Fibonacci usual. Como \(\det(Q)=-1\), la matriz \(Q\) es invertible módulo cualquier \(m\). Por eso cada estado ya pertenece a un ciclo, sin colas preperiódicas.
El periodo global para el módulo \(m\) es el orden de \(Q\) módulo \(m\), es decir, precisamente el periodo de Pisano \(\pi(m)\). Todo periodo individual \(P_m(a,b)\) divide a \(\pi(m)\). El problema deja de ser una simulación temporal y pasa a ser la distribución de los \(m^2\) estados entre los divisores de \(\pi(m)\).
Para un primo \(p\neq 2,5\), el polinomio característico de \(Q\) es \(t^2-t-1\), cuyo discriminante es \(5\). De ahí sale el orden candidato usado por las implementaciones:
$$\pi(p)\mid p-1 \quad \text{si } \left(\frac{5}{p}\right)=1,\qquad \pi(p)\mid 2(p+1) \quad \text{si } \left(\frac{5}{p}\right)=-1.$$
A partir de ese candidato, el código va quitando factores primos mientras la potencia correspondiente siga siendo la identidad módulo \(p\). En términos de Fibonacci, la prueba es
$$Q^n\equiv I \pmod p \iff F_n\equiv 0 \pmod p \text{ y } F_{n+1}\equiv 1 \pmod p.$$
Después, para potencias primas, se aplican exactamente las fórmulas de levantamiento presentes en las soluciones:
$$\pi(2)=3,\qquad \pi(4)=6,\qquad \pi(2^e)=3\cdot 2^{e-1}\ \ (e\ge 3),$$
$$\pi(5^e)=20\cdot 5^{e-1},\qquad \pi(p^e)=\pi(p)\,p^{e-1}\ \ (p\text{ impar},\ p\neq 5).$$
Todos los valores de Fibonacci necesarios se calculan con fast doubling:
$$F_{2k}=F_k(2F_{k+1}-F_k),\qquad F_{2k+1}=F_k^2+F_{k+1}^2,$$
de modo que cada consulta \(F_n\bmod m\) cuesta solo \(O(\log n)\).
Para recuperar los periodos exactos, primero se cuenta la cantidad más sencilla
$$A_{p^e}(d)=\#\{v\in(\mathbb{Z}/p^e\mathbb{Z})^2:Q^d v\equiv v\pmod{p^e}\},$$
para cada divisor \(d\mid \pi(p^e)\). Esto es el tamaño del núcleo de \(Q^d-I\). A partir de la identidad matricial anterior,
$$Q^d-I=\begin{pmatrix}F_{d-1}-1 & F_d\\ F_d & F_{d+1}-1\end{pmatrix}.$$
Sobre \(\mathbb{Z}/p^e\mathbb{Z}\), una matriz \(2\times 2\) queda descrita por su forma normal de Smith. El código extrae los dos exponentes invariantes \(p\)-ádicos a partir de las valuaciones de las entradas y del determinante. Como \(F_{d+1}-1=(F_{d-1}-1)+F_d\), la menor valuación de las entradas ya aparece en
$$\alpha=\min\bigl(v_p(F_d),\,v_p(F_{d-1}-1)\bigr).$$
Con
$$\Delta_d=\det(Q^d-I),\qquad \beta=\max\bigl(0,\,v_p(\Delta_d)-\alpha\bigr),$$
el número de soluciones de \((Q^d-I)v\equiv 0\pmod{p^e}\) es
$$A_{p^e}(d)=p^{\min(e,\alpha)+\min(e,\beta)}.$$
Esa es exactamente la fórmula de conteo de estados fijos que implementan los tres lenguajes. Los cálculos se hacen módulo \(p^{2e}\) para poder recuperar con seguridad valuaciones de hasta \(2e\).
Ahora escribimos
$$m=\prod_{i=1}^{r} p_i^{e_i},\qquad \pi(m)=\operatorname{lcm}_{1\le i\le r}\pi(p_i^{e_i}).$$
Por el teorema chino del resto, un estado módulo \(m\) equivale a una tupla de estados módulo las potencias primas, así que los conteos de estados fijos se multiplican:
$$A_m(d)=\prod_{i=1}^{r} A_{p_i^{e_i}}\!\bigl(\gcd(d,\pi(p_i^{e_i}))\bigr).$$
Sea \(E_m(n)\) el número de estados cuyo periodo exacto es \(n\). Entonces
$$A_m(d)=\sum_{n\mid d} E_m(n),$$
porque ser fijo por \(Q^d\) equivale a tener periodo exacto que divide a \(d\). La inversión de Möbius da
$$E_m(n)=\sum_{t\mid n}\mu\!\left(\frac{n}{t}\right)A_m(t).$$
Por tanto,
$$s(m)=\sum_{n\mid \pi(m)} n^2 E_m(n).$$
Al intercambiar el orden de sumación aparece la fórmula sobre divisores que evalúa el código:
$$s(m)=\sum_{d\mid \pi(m)} A_m(d)\,d^2 \sum_{k\mid \pi(m)/d}\mu(k)k^2.$$
La suma interior es multiplicativa y se colapsa como
$$\sum_{k\mid N}\mu(k)k^2=\prod_{p\mid N}(1-p^2).$$
Así, la fórmula de trabajo final es
$$\boxed{s(m)=\sum_{d\mid \pi(m)} A_m(d)\,d^2\prod_{p\mid \pi(m)/d}(1-p^2).}$$
Las implementaciones precomputan ese factor multiplicativo una sola vez y luego lo reutilizan.
Módulo \(2\), el periodo de Pisano es \(\pi(2)=3\). Hay cuatro estados:
$$ (0,0),\ (0,1),\ (1,0),\ (1,1). $$
El estado nulo queda fijo inmediatamente, luego su periodo es \(1\). Los otros tres estados forman un 3-ciclo bajo la acción de \(Q\), así que todos tienen periodo exacto \(3\). Por ello,
$$E_2(1)=1,\qquad E_2(3)=3,$$
y entonces
$$s(2)=1^2\cdot 1 + 3^2\cdot 3 = 28.$$
La fórmula de puntos fijos ve exactamente lo mismo: \(A_2(1)=1\) y \(A_2(3)=4\), y después la inversión de Möbius separa el único estado fijo de los tres estados con periodo genuino \(3\).
Las implementaciones en C++, Python y Java empiezan con una criba de menor factor primo hasta \(6M\). Eso basta para factorizar todos los valores de periodo que aparecerán. Luego calculan \(\pi(p)\) para cada primo \(p\le M\), lo extienden a todas las potencias primas \(p^e\le M\) y almacenan la lista completa de divisores de cada \(\pi(p^e)\).
Para cada divisor \(d\mid \pi(p^e)\), la implementación evalúa \(A_{p^e}(d)\) a partir de las valuaciones \(p\)-ádicas de las entradas de \(Q^d-I\). De ese modo, la parte costosa de la matemática queda convertida en una tabla de consulta indexada por potencia prima y divisor del periodo. La misma fase de preparación construye también el factor \(\prod_{p\mid N}(1-p^2)\) que aparece en la suma de Möbius colapsada.
Para calcular \(s(m)\), la implementación factoriza \(m\), construye \(\pi(m)\) como mínimo común múltiplo de los periodos de sus potencias primas, recorre todos los divisores \(d\mid \pi(m)\), multiplica los conteos de estados fijos correspondientes y les aplica el peso \(d^2\prod_{p\mid \pi(m)/d}(1-p^2)\). La suma de esos términos da \(s(m)\), y sumar \(s(m)\) para \(1\le m\le M\) produce \(S(M)\).
Las igualdades de validación incluidas en el código son \(s(3)=513\), \(s(10)=225820\), \(S(3)=542\) y \(S(10)=310897\). La canalización matemática es la misma en los tres lenguajes; la versión en C++ también realiza el cálculo completo para la cota grande y puede repartir el rango exterior de \(m\) entre varios hilos de CPU, mientras que las versiones en Python y Java validan el método en los puntos pequeños y luego emiten el residuo final ya conocido para la mayor cota publicada.
La criba, la lista de primos, las tablas de periodos de potencias primas y la tabla de pesos multiplicativos se precomputan en tiempo casi lineal en \(6M\). Después, para un módulo \(m\), el costo dominante es enumerar los divisores de \(\pi(m)\) y multiplicar los conteos de estados fijos correspondientes. No hace falta recorrer explícitamente los \(m^2\) estados iniciales.
La memoria está dominada por el arreglo del menor factor primo hasta \(6M\), la información de las potencias primas y las tablas de divisores y puntos fijos asociadas a sus periodos. La ganancia esencial es sustituir la simulación de órbitas por aritmética sobre divisores de periodos de Pisano.
对每个模数 \(m\)，把每个有序剩余对 \((a,b)\in(\mathbb{Z}/m\mathbb{Z})^2\) 视为一个 Fibonacci 型递推的初始状态：
$$x_0=a,\qquad x_1=b,\qquad x_{n+2}\equiv x_{n+1}+x_n \pmod m.$$
由于状态空间有限，\((x_n,x_{n+1})\) 终会重复。这里的转移矩阵行列式为 \(-1\)，因此它在任意模数下都可逆，所以每个状态从一开始就是纯周期的，没有前导尾段。记从 \((a,b)\) 出发的状态精确周期为 \(P_m(a,b)\)，再定义
$$s(m)=\sum_{a,b \bmod m} P_m(a,b)^2,\qquad S(M)=\sum_{m=1}^{M}s(m).$$
第 947 题要求计算 \(S(10^6)\) 在模 \(999999893\) 下的值。若直接做状态模拟，就要对每个 \(m\le 10^6\) 枚举至多 \(m^2\) 个起点并追踪其循环，这显然不可行。代码真正采用的方法，是把递推写成矩阵作用，在素数幂模数上统计不动状态，再通过除数公式恢复精确周期。
整套推导都围绕 Fibonacci 转移矩阵
$$Q=\begin{pmatrix}0&1\\1&1\end{pmatrix}$$
展开。若把状态向量写成 \(v_n=(x_n,x_{n+1})^T\)，则有 \(v_{n+1}=Qv_n\)，于是 \(v_n=Q^n v_0\)。这样一来，代码中真正处理的数学对象就变得非常清楚：矩阵的阶、不动向量以及轨道的精确长度。
经典的 Fibonacci 矩阵恒等式是
$$Q^d=\begin{pmatrix}F_{d-1}&F_d\\F_d&F_{d+1}\end{pmatrix},$$
其中 \(F_n\) 是普通 Fibonacci 数。由于 \(\det(Q)=-1\)，矩阵 \(Q\) 在任意模 \(m\) 下都可逆，因此每个状态都落在某个循环上，不会先经过一段非周期尾巴。
模数 \(m\) 的全局周期，就是 \(Q\) 在模 \(m\) 下的阶，这正是 Pisano 周期 \(\pi(m)\)。每个单独状态的周期 \(P_m(a,b)\) 都必须整除 \(\pi(m)\)。因此问题并不是“把时间推进到什么时候会重复”，而是“这 \(m^2\) 个状态如何分布到 \(\pi(m)\) 的各个因子上”。
当 \(p\neq 2,5\) 为素数时，\(Q\) 的特征多项式是 \(t^2-t-1\)，判别式为 \(5\)。这正是实现中候选阶来源：
$$\pi(p)\mid p-1 \quad \text{若 } \left(\frac{5}{p}\right)=1,\qquad \pi(p)\mid 2(p+1) \quad \text{若 } \left(\frac{5}{p}\right)=-1.$$
代码先从该候选值出发，再逐个剥去其素因子；只要去掉某个因子后矩阵幂仍然是单位矩阵，就继续缩小。用 Fibonacci 数表示，这个恒等条件就是
$$Q^n\equiv I \pmod p \iff F_n\equiv 0 \pmod p \text{ 且 } F_{n+1}\equiv 1 \pmod p.$$
对素数幂，程序使用的正是标准提升公式：
$$\pi(2)=3,\qquad \pi(4)=6,\qquad \pi(2^e)=3\cdot 2^{e-1}\ \ (e\ge 3),$$
$$\pi(5^e)=20\cdot 5^{e-1},\qquad \pi(p^e)=\pi(p)\,p^{e-1}\ \ (p\text{ 为奇素数},\ p\neq 5).$$
所有需要的 Fibonacci 模值都由 fast doubling 计算：
$$F_{2k}=F_k(2F_{k+1}-F_k),\qquad F_{2k+1}=F_k^2+F_{k+1}^2,$$
因此每次 \(F_n\bmod m\) 查询只需 \(O(\log n)\) 时间。
为了恢复精确周期，程序先统计更容易的量
$$A_{p^e}(d)=\#\{v\in(\mathbb{Z}/p^e\mathbb{Z})^2:Q^d v\equiv v\pmod{p^e}\},$$
这里 \(d\mid \pi(p^e)\)。这就是矩阵 \(Q^d-I\) 核的大小。由上面的矩阵恒等式可得
$$Q^d-I=\begin{pmatrix}F_{d-1}-1 & F_d\\ F_d & F_{d+1}-1\end{pmatrix}.$$
在 \(\mathbb{Z}/p^e\mathbb{Z}\) 上，一个 \(2\times 2\) 矩阵的结构可由 Smith 标准形控制。实现里通过条目及其行列式的 \(p\)-进赋值提取出两个不变量指数。因为 \(F_{d+1}-1=(F_{d-1}-1)+F_d\)，最小的条目赋值已经包含在
$$\alpha=\min\bigl(v_p(F_d),\,v_p(F_{d-1}-1)\bigr)$$
之中。再记
$$\Delta_d=\det(Q^d-I),\qquad \beta=\max\bigl(0,\,v_p(\Delta_d)-\alpha\bigr),$$
则线性方程组 \((Q^d-I)v\equiv 0\pmod{p^e}\) 的解数为
$$A_{p^e}(d)=p^{\min(e,\alpha)+\min(e,\beta)}.$$
这正是三份实现共同使用的不动状态公式。为了安全恢复最高到 \(2e\) 的赋值，相关计算都在模 \(p^{2e}\) 下进行。
现在设
$$m=\prod_{i=1}^{r} p_i^{e_i},\qquad \pi(m)=\operatorname{lcm}_{1\le i\le r}\pi(p_i^{e_i}).$$
由中国剩余定理，模 \(m\) 的一个状态等价于各素数幂模数状态的一个元组，因此不动状态数可分解为乘积：
$$A_m(d)=\prod_{i=1}^{r} A_{p_i^{e_i}}\!\bigl(\gcd(d,\pi(p_i^{e_i}))\bigr).$$
记 \(E_m(n)\) 为精确周期等于 \(n\) 的状态个数，则有
$$A_m(d)=\sum_{n\mid d} E_m(n),$$
因为被 \(Q^d\) 固定当且仅当精确周期整除 \(d\)。于是 Möbius 反演给出
$$E_m(n)=\sum_{t\mid n}\mu\!\left(\frac{n}{t}\right)A_m(t).$$
因此
$$s(m)=\sum_{n\mid \pi(m)} n^2 E_m(n).$$
交换求和次序后，就得到代码真正计算的除数公式：
$$s(m)=\sum_{d\mid \pi(m)} A_m(d)\,d^2 \sum_{k\mid \pi(m)/d}\mu(k)k^2.$$
其中内层和式是乘法性的，并可压缩成
$$\sum_{k\mid N}\mu(k)k^2=\prod_{p\mid N}(1-p^2).$$
于是最终工作公式为
$$\boxed{s(m)=\sum_{d\mid \pi(m)} A_m(d)\,d^2\prod_{p\mid \pi(m)/d}(1-p^2).}$$
程序把最后这个乘法因子预处理出来，然后在所有模数上重复使用。
在模 \(2\) 下，Pisano 周期是 \(\pi(2)=3\)。总共有四个状态：
$$ (0,0),\ (0,1),\ (1,0),\ (1,1). $$
零状态立即固定，因此周期为 \(1\)。另外三个状态在 \(Q\) 的作用下形成一个 3-循环，所以它们的精确周期都是 \(3\)。于是
$$E_2(1)=1,\qquad E_2(3)=3,$$
从而
$$s(2)=1^2\cdot 1 + 3^2\cdot 3 = 28.$$
不动点公式看到的正是同一件事：\(A_2(1)=1\)，\(A_2(3)=4\)。随后 Möbius 反演把那个真正固定的状态与三个周期恰为 \(3\) 的状态分离开来。
C++、Python 和 Java 实现都先建立一个到 \(6M\) 为止的最小质因子筛。这样就足以分解后续遇到的所有周期值。随后它们为每个 \(p\le M\) 的素数求出 \(\pi(p)\)，再扩展到所有 \(p^e\le M\) 的素数幂，并保存每个 \(\pi(p^e)\) 的完整因子列表。
对于每个 \(d\mid \pi(p^e)\)，实现都会根据 \(Q^d-I\) 各条目的 \(p\)-进赋值计算 \(A_{p^e}(d)\)。这样一来，最昂贵的数学部分就变成了“素数幂 + 周期因子”索引下的查表问题。同一阶段还会预先构造 \(\prod_{p\mid N}(1-p^2)\)，也就是压缩后的 Möbius 权重。
在求 \(s(m)\) 时，程序先分解 \(m\)，把相关素数幂周期取最小公倍数得到 \(\pi(m)\)，枚举所有 \(d\mid \pi(m)\)，将对应的不动状态数相乘，再乘上权重 \(d^2\prod_{p\mid \pi(m)/d}(1-p^2)\)。这些项相加就是 \(s(m)\)，继续对 \(1\le m\le M\) 求和即可得到 \(S(M)\)。
代码中的校验恒等式是 \(s(3)=513\)、\(s(10)=225820\)、\(S(3)=542\)、\(S(10)=310897\)。三种语言采用的是同一条数学管线；其中 C++ 版本还能真正完成大规模上界的全部计算，并可把外层 \(m\) 区间拆给多个 CPU 线程，而 Python 与 Java 版本会先在这些较小检查点上验证方法，再对最大的公开边界直接输出已知的最终余数。
筛法、素数表、素数幂周期表以及乘法权重表，都可以在关于 \(6M\) 的近线性时间内预处理完成。之后，单个模数 \(m\) 的主要成本来自枚举 \(\pi(m)\) 的因子并相乘对应的不动点数量。算法完全不需要显式遍历全部 \(m^2\) 个初始状态。
内存主要消耗在到 \(6M\) 为止的最小质因子数组、素数幂元数据以及各个素数幂周期对应的因子表与不动点表上。真正的收益在于：它把轨道模拟替换成了 Pisano 周期因子上的算术运算。
Для каждого модуля \(m\) берется каждая упорядоченная пара вычетов \((a,b)\in(\mathbb{Z}/m\mathbb{Z})^2\) как начальное состояние рекурсии фибоначчиева типа
$$x_0=a,\qquad x_1=b,\qquad x_{n+2}\equiv x_{n+1}+x_n \pmod m.$$
Так как пространство состояний конечно, пара \((x_n,x_{n+1})\) рано или поздно повторяется. Здесь матрица перехода имеет определитель \(-1\), поэтому она обратима по любому модулю, и все состояния сразу лежат на циклах, без предпериодической части. Обозначим через \(P_m(a,b)\) точный период состояния, стартующего из \((a,b)\), и определим
$$s(m)=\sum_{a,b \bmod m} P_m(a,b)^2,\qquad S(M)=\sum_{m=1}^{M}s(m).$$
В задаче 947 требуется вычислить \(S(10^6)\) по модулю \(999999893\). Наивная симуляция означала бы отслеживание до \(m^2\) стартовых состояний для каждого \(m\le 10^6\), что нереально. Реальные решения переводят рекурсию в действие матрицы, считают неподвижные состояния на простых степенях и восстанавливают точные периоды через формулу по делителям.
Центральный объект всей конструкции - матрица перехода Фибоначчи
$$Q=\begin{pmatrix}0&1\\1&1\end{pmatrix}.$$
Если записать вектор состояния как \(v_n=(x_n,x_{n+1})^T\), то \(v_{n+1}=Qv_n\), а значит \(v_n=Q^n v_0\). Именно так становятся видны те объекты, с которыми фактически работает код: порядки матриц, неподвижные векторы и точные длины орбит.
Стандартная матричная формула для чисел Фибоначчи имеет вид
$$Q^d=\begin{pmatrix}F_{d-1}&F_d\\F_d&F_{d+1}\end{pmatrix},$$
где \(F_n\) - обычные числа Фибоначчи. Поскольку \(\det(Q)=-1\), матрица \(Q\) обратима по модулю любого \(m\). Поэтому каждое состояние уже находится на цикле, и никакого хвоста перед циклом нет.
Глобальный период для модуля \(m\) - это порядок матрицы \(Q\) по модулю \(m\), то есть ровно период Пизано \(\pi(m)\). Любой индивидуальный период \(P_m(a,b)\) делит \(\pi(m)\). Значит, задача сводится не к длительной симуляции, а к разбиению \(m^2\) состояний по делителям \(\pi(m)\).
Для простого \(p\neq 2,5\) характеристический многочлен матрицы \(Q\) равен \(t^2-t-1\), его дискриминант равен \(5\). Именно отсюда возникает кандидат на порядок, который использует код:
$$\pi(p)\mid p-1 \quad \text{если } \left(\frac{5}{p}\right)=1,\qquad \pi(p)\mid 2(p+1) \quad \text{если } \left(\frac{5}{p}\right)=-1.$$
Начиная с этого кандидата, реализация по одному убирает простые множители, пока соответствующая степень матрицы все еще дает единичную матрицу по модулю \(p\). В терминах Фибоначчи критерий единичности таков:
$$Q^n\equiv I \pmod p \iff F_n\equiv 0 \pmod p \text{ и } F_{n+1}\equiv 1 \pmod p.$$
Затем для степеней простых используются ровно те формулы подъема, которые присутствуют в решениях:
$$\pi(2)=3,\qquad \pi(4)=6,\qquad \pi(2^e)=3\cdot 2^{e-1}\ \ (e\ge 3),$$
$$\pi(5^e)=20\cdot 5^{e-1},\qquad \pi(p^e)=\pi(p)\,p^{e-1}\ \ (p\text{ нечетный},\ p\neq 5).$$
Все необходимые значения Фибоначчи вычисляются методом fast doubling:
$$F_{2k}=F_k(2F_{k+1}-F_k),\qquad F_{2k+1}=F_k^2+F_{k+1}^2,$$
поэтому каждый запрос \(F_n\bmod m\) стоит лишь \(O(\log n)\).
Чтобы восстановить точные периоды, сначала считается более простая величина
$$A_{p^e}(d)=\#\{v\in(\mathbb{Z}/p^e\mathbb{Z})^2:Q^d v\equiv v\pmod{p^e}\},$$
для каждого делителя \(d\mid \pi(p^e)\). Это размер ядра матрицы \(Q^d-I\). Из матричной формулы выше следует
$$Q^d-I=\begin{pmatrix}F_{d-1}-1 & F_d\\ F_d & F_{d+1}-1\end{pmatrix}.$$
Над \(\mathbb{Z}/p^e\mathbb{Z}\) матрица \(2\times 2\) описывается своей нормальной формой Смита. Код извлекает два \(p\)-адических инвариантных показателя из оценок элементов и определителя. Поскольку \(F_{d+1}-1=(F_{d-1}-1)+F_d\), минимальная оценка среди элементов уже содержится в
$$\alpha=\min\bigl(v_p(F_d),\,v_p(F_{d-1}-1)\bigr).$$
Если положить
$$\Delta_d=\det(Q^d-I),\qquad \beta=\max\bigl(0,\,v_p(\Delta_d)-\alpha\bigr),$$
то число решений системы \((Q^d-I)v\equiv 0\pmod{p^e}\) равно
$$A_{p^e}(d)=p^{\min(e,\alpha)+\min(e,\beta)}.$$
Это в точности та формула для числа неподвижных состояний, которую реализуют все три версии. Вычисления ведутся по модулю \(p^{2e}\), чтобы надежно восстанавливать оценки вплоть до \(2e\).
Пусть теперь
$$m=\prod_{i=1}^{r} p_i^{e_i},\qquad \pi(m)=\operatorname{lcm}_{1\le i\le r}\pi(p_i^{e_i}).$$
По китайской теореме об остатках состояние по модулю \(m\) эквивалентно набору состояний по модулям простых степеней, поэтому число неподвижных состояний раскладывается в произведение:
$$A_m(d)=\prod_{i=1}^{r} A_{p_i^{e_i}}\!\bigl(\gcd(d,\pi(p_i^{e_i}))\bigr).$$
Обозначим через \(E_m(n)\) число состояний с точным периодом \(n\). Тогда
$$A_m(d)=\sum_{n\mid d} E_m(n),$$
потому что быть фиксированным матрицей \(Q^d\) значит иметь точный период, делящий \(d\). Формула Мёбиуса дает
$$E_m(n)=\sum_{t\mid n}\mu\!\left(\frac{n}{t}\right)A_m(t).$$
Следовательно,
$$s(m)=\sum_{n\mid \pi(m)} n^2 E_m(n).$$
После перестановки сумм получается именно та формула по делителям, которую вычисляет код:
$$s(m)=\sum_{d\mid \pi(m)} A_m(d)\,d^2 \sum_{k\mid \pi(m)/d}\mu(k)k^2.$$
Внутренняя сумма мультипликативна и схлопывается в
$$\sum_{k\mid N}\mu(k)k^2=\prod_{p\mid N}(1-p^2).$$
Поэтому рабочая формула имеет вид
$$\boxed{s(m)=\sum_{d\mid \pi(m)} A_m(d)\,d^2\prod_{p\mid \pi(m)/d}(1-p^2).}$$
Именно этот последний множитель реализации предварительно вычисляют один раз и потом переиспользуют.
По модулю \(2\) период Пизано равен \(\pi(2)=3\). Всего имеется четыре состояния:
$$ (0,0),\ (0,1),\ (1,0),\ (1,1). $$
Нулевое состояние сразу неподвижно, значит его период равен \(1\). Остальные три состояния образуют 3-цикл под действием \(Q\), так что у них точный период \(3\). Следовательно,
$$E_2(1)=1,\qquad E_2(3)=3,$$
и потому
$$s(2)=1^2\cdot 1 + 3^2\cdot 3 = 28.$$
Формула для неподвижных состояний видит то же самое: \(A_2(1)=1\) и \(A_2(3)=4\), а затем инверсия Мёбиуса отделяет единственное действительно неподвижное состояние от трех состояний точного периода \(3\).
Реализации на C++, Python и Java начинают с решета наименьшего простого делителя до \(6M\). Этого достаточно, чтобы разложить на множители все встречающиеся значения периодов. Затем они вычисляют \(\pi(p)\) для каждого простого \(p\le M\), поднимают эти значения на все степени простых \(p^e\le M\) и сохраняют полный список делителей каждого \(\pi(p^e)\).
Для каждого делителя \(d\mid \pi(p^e)\) программа вычисляет \(A_{p^e}(d)\) по \(p\)-адическим оценкам элементов матрицы \(Q^d-I\). Тем самым самая дорогая часть математики превращается в таблицу поиска, индексируемую простой степенью и делителем ее периода. На том же этапе строится и множитель \(\prod_{p\mid N}(1-p^2)\), возникающий после свертки суммы Мёбиуса.
Чтобы найти \(s(m)\), реализация раскладывает \(m\) на множители, строит \(\pi(m)\) как наименьшее общее кратное соответствующих периодов простых степеней, перебирает все делители \(d\mid \pi(m)\), перемножает нужные количества неподвижных состояний и умножает их на вес \(d^2\prod_{p\mid \pi(m)/d}(1-p^2)\). Сумма этих членов дает \(s(m)\), а последующее суммирование по \(1\le m\le M\) дает \(S(M)\).
Проверочные тождества в коде таковы: \(s(3)=513\), \(s(10)=225820\), \(S(3)=542\) и \(S(10)=310897\). Математический конвейер одинаков во всех трех языках; версия на C++ также выполняет полный большой расчет и при необходимости делит внешний диапазон по \(m\) между несколькими потоками CPU, а версии на Python и Java проверяют метод на малых контрольных точках, после чего выводят уже известный итоговый остаток для наибольшей опубликованной границы.
Решето, список простых, таблицы периодов простых степеней и таблица мультипликативных весов предвычисляются за почти линейное время от \(6M\). После этого для одного модуля \(m\) основная стоимость определяется перебором делителей числа \(\pi(m)\) и перемножением соответствующих количеств неподвижных состояний. Явно перебирать все \(m^2\) стартовых состояния не требуется.
Память в основном расходуется на массив наименьших простых делителей до \(6M\), метаданные по простым степеням и таблицы делителей и неподвижных состояний для периодов простых степеней. Практический выигрыш состоит в замене симуляции орбит арифметикой на делителях периодов Пизано.
لكل عدد \(m\)، نأخذ كل زوج مرتب من البواقي \((a,b)\in(\mathbb{Z}/m\mathbb{Z})^2\) كنقطة بداية لمتتالية من نمط فيبوناتشي:
$$x_0=a,\qquad x_1=b,\qquad x_{n+2}\equiv x_{n+1}+x_n \pmod m.$$
بما أن فضاء الحالات منتهٍ، فإن الزوج \((x_n,x_{n+1})\) لا بد أن يتكرر. وفي هذه المسألة تكون مصفوفة الانتقال ذات محدد يساوي \(-1\)، لذلك فهي قابلة للعكس بترديد أي \(m\)، ومن ثم تكون كل حالة دورية خالصة منذ البداية من غير جزء تمهيدي. إذا رمزنا إلى الدورة الدقيقة للحالة المبدوءة من \((a,b)\) بـ \(P_m(a,b)\)، فإننا نعرّف
$$s(m)=\sum_{a,b \bmod m} P_m(a,b)^2,\qquad S(M)=\sum_{m=1}^{M}s(m).$$
المطلوب في المسألة 947 هو حساب \(S(10^6)\) بترديد \(999999893\). أما المحاكاة المباشرة فستحاول تتبع ما يصل إلى \(m^2\) حالة ابتدائية لكل \(m\le 10^6\)، وهذا غير عملي إطلاقًا. لذلك تحوّل الحلول العلاقة العودية إلى فعل مصفوفي، وتحسب أعداد الحالات الثابتة على قوى الأعداد الأولية، ثم تستعيد الدورات الدقيقة بصيغة مبنية على القواسم.
يرتكز الحل كله على مصفوفة انتقال فيبوناتشي
$$Q=\begin{pmatrix}0&1\\1&1\end{pmatrix}.$$
إذا كتبنا متجه الحالة على الصورة \(v_n=(x_n,x_{n+1})^T\)، فإن \(v_{n+1}=Qv_n\)، وبالتالي \(v_n=Q^n v_0\). بهذه الصياغة تظهر مباشرة الكيانات الرياضية التي يعالجها التنفيذ فعلًا: رتب المصفوفات، والمتجهات الثابتة، وأطوال المدارات الدقيقة.
المتطابقة المصفوفية القياسية لفibonacci هي
$$Q^d=\begin{pmatrix}F_{d-1}&F_d\\F_d&F_{d+1}\end{pmatrix},$$
حيث \(F_n\) هي أعداد فيبوناتشي المعتادة. وبما أن \(\det(Q)=-1\)، فإن \(Q\) قابلة للعكس بترديد أي \(m\). لذلك تقع كل حالة على دورة مباشرة، ولا توجد مرحلة ما قبل الدورية.
الدورة العامة للموديل \(m\) هي رتبة \(Q\) بترديد \(m\)، وهذا بالضبط هو دور بيسانو \(\pi(m)\). ومن ثم فإن أي دورة فردية \(P_m(a,b)\) لا بد أن تقسم \(\pi(m)\). أي إن المسألة ليست تتبع الزمن خطوة خطوة، بل فهم كيف تتوزع الحالات وعددها \(m^2\) على قواسم \(\pi(m)\).
إذا كان \(p\neq 2,5\) عددًا أوليًا، فإن كثير الحدود المميز لـ \(Q\) هو \(t^2-t-1\)، ومميزه يساوي \(5\). ومن هنا يأتي المرشح الأولي للرتبة الذي تستخدمه الحلول:
$$\pi(p)\mid p-1 \quad \text{if } \left(\frac{5}{p}\right)=1,\qquad \pi(p)\mid 2(p+1) \quad \text{if } \left(\frac{5}{p}\right)=-1.$$
ينطلق التنفيذ من هذا المرشح، ثم يحذف عوامله الأولية واحدًا بعد آخر ما دامت القوة الناتجة لا تزال تعطي مصفوفة الوحدة بترديد \(p\). وبصياغة فيبوناتشي، يكون اختبار الهوية هو
$$Q^n\equiv I \pmod p \iff F_n\equiv 0 \pmod p \text{ and } F_{n+1}\equiv 1 \pmod p.$$
ثم تُعالج قوى الأعداد الأولية بواسطة صيغ الرفع القياسية نفسها الموجودة في الشفرات:
$$\pi(2)=3,\qquad \pi(4)=6,\qquad \pi(2^e)=3\cdot 2^{e-1}\ \ (e\ge 3),$$
$$\pi(5^e)=20\cdot 5^{e-1},\qquad \pi(p^e)=\pi(p)\,p^{e-1}\ \ (p\text{ odd prime},\ p\neq 5).$$
وجميع قيم فيبوناتشي المطلوبة تُحسب بطريقة fast doubling:
$$F_{2k}=F_k(2F_{k+1}-F_k),\qquad F_{2k+1}=F_k^2+F_{k+1}^2,$$
وبذلك تصبح كل قيمة \(F_n\bmod m\) قابلة للحساب في \(O(\log n)\) فقط.
لاستخراج الدورات الدقيقة، يبدأ الحل بعدّ الكمية الأسهل
$$A_{p^e}(d)=\#\{v\in(\mathbb{Z}/p^e\mathbb{Z})^2:Q^d v\equiv v\pmod{p^e}\},$$
لكل قاسم \(d\mid \pi(p^e)\). هذا هو حجم نواة المصفوفة \(Q^d-I\). ومن المتطابقة المصفوفية السابقة نحصل على
$$Q^d-I=\begin{pmatrix}F_{d-1}-1 & F_d\\ F_d & F_{d+1}-1\end{pmatrix}.$$
فوق الحلقة \(\mathbb{Z}/p^e\mathbb{Z}\)، تتحكم صيغة سميث المعيارية في بنية مصفوفة \(2\times 2\). يستخرج التنفيذ أسّي الثابتين \(p\)-الأديين من تقييمات العناصر ومن تقييم المحدد. وبما أن \(F_{d+1}-1=(F_{d-1}-1)+F_d\)، فإن أصغر تقييم بين العناصر يظهر أصلًا في
$$\alpha=\min\bigl(v_p(F_d),\,v_p(F_{d-1}-1)\bigr).$$
ومع التعريفين
$$\Delta_d=\det(Q^d-I),\qquad \beta=\max\bigl(0,\,v_p(\Delta_d)-\alpha\bigr),$$
يكون عدد حلول المعادلة \((Q^d-I)v\equiv 0\pmod{p^e}\) هو
$$A_{p^e}(d)=p^{\min(e,\alpha)+\min(e,\beta)}.$$
وهذه بالضبط هي صيغة عدّ الحالات الثابتة التي تطبقها اللغات الثلاث. وتُجرى الحسابات بترديد \(p^{2e}\) حتى يمكن استرجاع التقييمات حتى \(2e\) بأمان.
لنكتب الآن
$$m=\prod_{i=1}^{r} p_i^{e_i},\qquad \pi(m)=\operatorname{lcm}_{1\le i\le r}\pi(p_i^{e_i}).$$
وبحسب مبرهنة البواقي الصينية، فإن الحالة بترديد \(m\) تكافئ حزمة من الحالات بترديدات القوى الأولية، ولذلك يتفكك عدد الحالات الثابتة إلى حاصل ضرب:
$$A_m(d)=\prod_{i=1}^{r} A_{p_i^{e_i}}\!\bigl(\gcd(d,\pi(p_i^{e_i}))\bigr).$$
إذا رمزنا بعدد الحالات ذات الدورة الدقيقة \(n\) بـ \(E_m(n)\)، فلدينا
$$A_m(d)=\sum_{n\mid d} E_m(n),$$
لأن ثبات الحالة تحت \(Q^d\) يكافئ تمامًا أن تكون دورتها الدقيقة قاسمًا لـ \(d\). ومن ثم تعطي صيغة Möbius العكسية
$$E_m(n)=\sum_{t\mid n}\mu\!\left(\frac{n}{t}\right)A_m(t).$$
وعليه
$$s(m)=\sum_{n\mid \pi(m)} n^2 E_m(n).$$
وبتبديل ترتيب الجمع نحصل على صيغة القواسم التي ينفذها الكود فعليًا:
$$s(m)=\sum_{d\mid \pi(m)} A_m(d)\,d^2 \sum_{k\mid \pi(m)/d}\mu(k)k^2.$$
والمجموع الداخلي ضربي، فينهار إلى
$$\sum_{k\mid N}\mu(k)k^2=\prod_{p\mid N}(1-p^2).$$
إذًا تصبح الصيغة العملية النهائية
$$\boxed{s(m)=\sum_{d\mid \pi(m)} A_m(d)\,d^2\prod_{p\mid \pi(m)/d}(1-p^2).}$$
ولهذا السبب يسبق التنفيذ بحساب هذا العامل الضربي الأخير مرة واحدة ثم يعيد استخدامه مع كل \(m\).
بترديد \(2\)، تكون فترة بيسانو \(\pi(2)=3\). ويوجد أربع حالات فقط:
$$ (0,0),\ (0,1),\ (1,0),\ (1,1). $$
الحالة الصفرية ثابتة فورًا، ومن ثم دورتها \(1\). أما الحالات الثلاث الأخرى فتشكّل دورة بطول 3 تحت تأثير \(Q\)، لذا تكون دورتها الدقيقة \(3\). إذن
$$E_2(1)=1,\qquad E_2(3)=3,$$
ومن ثم
$$s(2)=1^2\cdot 1 + 3^2\cdot 3 = 28.$$
وترى صيغة الحالات الثابتة النتيجة نفسها بالضبط: \(A_2(1)=1\) و\(A_2(3)=4\)، ثم تفصل Möbius inversion بين الحالة الثابتة الوحيدة وبين الحالات الثلاث ذات الدورة الحقيقية \(3\).
تبدأ تطبيقات C++ وPython وJava بمنخل لأصغر عامل أولي حتى \(6M\). وهذا يكفي لتحليل جميع قيم الفترات التي ستظهر. بعد ذلك تُحسب \(\pi(p)\) لكل أولي \(p\le M\)، ثم تُرفع هذه القيم إلى كل قوة أولية \(p^e\le M\)، مع تخزين قائمة القواسم الكاملة لكل \(\pi(p^e)\).
لكل قاسم \(d\mid \pi(p^e)\)، يحسب التنفيذ \(A_{p^e}(d)\) انطلاقًا من التقييمات \(p\)-الأدية لعناصر المصفوفة \(Q^d-I\). وبهذا تتحول أثقل خطوة رياضية إلى جدول بحث مفهرس بالقوة الأولية وقاسم الفترة. وفي المرحلة نفسها يُبنى أيضًا العامل \(\prod_{p\mid N}(1-p^2)\) الذي يظهر بعد ضغط مجموع Möbius.
لحساب \(s(m)\)، تحلل الشفرة \(m\) إلى عوامله، وتكوّن \(\pi(m)\) على أنه المضاعف المشترك الأصغر لفترات القوى الأولية الداخلة فيه، ثم تعدد جميع القواسم \(d\mid \pi(m)\)، وتضرب أعداد الحالات الثابتة المناسبة، وتزنها بالعامل \(d^2\prod_{p\mid \pi(m)/d}(1-p^2)\). مجموع هذه الحدود هو \(s(m)\)، ثم يعطي جمع \(s(m)\) على \(1\le m\le M\) القيمة \(S(M)\).
المتطابقات التحققية داخل الشفرة هي \(s(3)=513\)، و\(s(10)=225820\)، و\(S(3)=542\)، و\(S(10)=310897\). والمسار الرياضي نفسه موجود في اللغات الثلاث؛ غير أن نسخة C++ تنفذ كذلك الحساب الكامل للحد الكبير ويمكنها تقسيم المجال الخارجي الخاص بـ \(m\) على عدة خيوط CPU، بينما تتحقق نسختا Python وJava من المنهج على النقاط الصغيرة ثم تطبعان الباقي النهائي المعروف مسبقًا عند أكبر حد منشور.
يُحسب المنخل، وقائمة الأعداد الأولية، وجداول فترات القوى الأولية، وجدول الأوزان الضربية في زمن شبه خطي بالنسبة إلى \(6M\). وبعد ذلك تصبح الكلفة الأساسية لمقدار واحد \(m\) هي تعداد قواسم \(\pi(m)\) وضرب أعداد الحالات الثابتة الموافقة لها. ولا توجد حاجة إلى المرور الصريح على جميع الحالات الابتدائية وعددها \(m^2\).
أما الذاكرة فتستهلكها أساسًا مصفوفة أصغر عامل أولي حتى \(6M\)، وبيانات القوى الأولية، وجداول القواسم والحالات الثابتة الخاصة بفترات القوى الأولية. والمكسب العملي الحاسم هو استبدال محاكاة المدارات بحسابات عددية على قواسم فترات بيسانو.