Define
$$M=123456789,\qquad x_1=1,\qquad x_{n+1}\equiv (x_n-1)^3+2 \pmod{M},\qquad 0\le x_n<M.$$
The target quantity is the gcd-weighted double sum
$$G(N)=\sum_{1\le i,j\le N} w_{\gcd(i,j)} \pmod{M},\qquad w_n=x_{x_n},$$
with the Project Euler input \(N=10^8\). A direct evaluation is impossible: the outer expression ranges over \(N^2\) ordered pairs, while each weight \(w_n\) is itself a nested lookup inside the same cubic recurrence. The solution works only because the gcd part can be counted arithmetically and the recurrence part becomes periodic.
The derivation has two main pieces. First, the gcd-weighted double sum is reorganized by possible gcd values and reduced to the summatory totient function. Second, the cubic orbit is shown to have a short preperiod and a manageable cycle, which turns the nested weight \(w_n=x_{x_n}\) into a table lookup inside one periodic tail.
For a fixed ordered pair \((i,j)\), write
$$i=d\,a,\qquad j=d\,b,\qquad d=\gcd(i,j),\qquad \gcd(a,b)=1.$$
Once \(d\) is fixed, the reduced pair \((a,b)\) must satisfy
$$1\le a,b\le \left\lfloor \frac{N}{d}\right\rfloor.$$
If we define
$$R(q)=\#\{(a,b)\in [1,q]^2:\gcd(a,b)=1\},$$
then every ordered pair whose gcd is exactly \(d\) contributes the same weight \(w_d\). Therefore
$$G(N)=\sum_{d=1}^{N} w_d\,R\!\left(\left\lfloor \frac{N}{d}\right\rfloor\right).$$
This is the first decisive simplification: the original two-dimensional sum is replaced by a one-dimensional sum over gcd classes.
The quantity \(R(q)\) counts ordered coprime pairs inside a \(q\times q\) grid. For each denominator \(m\ge 2\), there are exactly \(\varphi(m)\) reduced numerators \(1\le a<m\). Because the pairs are ordered, both \((a,m)\) and \((m,a)\) appear, and the diagonal contributes the extra point \((1,1)\). Hence
$$R(q)=1+2\sum_{m=2}^{q}\varphi(m)=2\Phi(q)-1,$$
where
$$\Phi(q)=\sum_{m=1}^{q}\varphi(m).$$
Substituting into the gcd decomposition gives the exact formula used by the implementations:
$$\boxed{G(N)=\sum_{d=1}^{N} w_d\left(2\Phi\!\left(\left\lfloor \frac{N}{d}\right\rfloor\right)-1\right)\pmod{M}.}$$
The remaining arithmetic bottleneck is the summatory totient \(\Phi(n)\). The standard identity behind the code is
$$\sum_{k=1}^{n}\Phi\!\left(\left\lfloor \frac{n}{k}\right\rfloor\right)=\sum_{m=1}^{n}\varphi(m)\left\lfloor \frac{n}{m}\right\rfloor=\sum_{t=1}^{n}\sum_{m\mid t}\varphi(m)=\sum_{t=1}^{n} t=\frac{n(n+1)}{2}.$$
Isolating the \(k=1\) term yields
$$\Phi(n)=\frac{n(n+1)}{2}-\sum_{k=2}^{n}\Phi\!\left(\left\lfloor \frac{n}{k}\right\rfloor\right).$$
The quotient \(\left\lfloor n/k\right\rfloor\) stays constant over intervals. If \(q=\left\lfloor n/\ell\right\rfloor\), then the last index with the same quotient is
$$r=\left\lfloor \frac{n}{q}\right\rfloor=\left\lfloor \frac{n}{\lfloor n/\ell\rfloor}\right\rfloor.$$
Grouping equal quotients gives the block form
$$\Phi(n)=\frac{n(n+1)}{2}-\sum_{\ell=2}^{n}(r-\ell+1)\,\Phi\!\left(\left\lfloor \frac{n}{\ell}\right\rfloor\right).$$
This is what makes large values of \(\Phi\) practical: the recursion touches only the distinct quotient values, not all integers one by one.
The recurrence
$$f(t)\equiv (t-1)^3+2 \pmod{M}$$
acts on a finite state space, so the orbit starting at \(x_1=1\) must eventually repeat. Floyd's cycle-finding method, applied to the actual orbit, gives a preperiod of length \(53\) and a cycle of length \(33705\). Therefore the periodic tail begins at \(x_{54}\), and for every \(r\ge 54\),
$$x_r=x_{\,54+((r-54)\bmod 33705)}.$$
The implementations therefore precompute the orbit through \(x_{33759}\), which is enough to cover the whole non-periodic prefix and every periodic representative that is later queried.
The nested index is the subtle part of the problem. Since the same cubic polynomial defines the sequence, reducing the state modulo the cycle length \(33705\) is compatible with iteration. Define the companion orbit
$$y_1=1,\qquad y_{n+1}\equiv (y_n-1)^3+2 \pmod{33705}.$$
Then \(y_n\equiv x_n\pmod{33705}\) for every \(n\). On the actual orbit, the only inner indices that need the non-periodic prefix are the first four:
$$x_1=1,\qquad x_2=2,\qquad x_3=3,\qquad x_4=10.$$
After these exceptional cases, the inner index is lifted back into the periodic tail by
$$I_n=54+((y_n-54)\bmod 33705),\qquad w_n=x_{I_n}.$$
So the huge-looking object \(x_{x_n}\) is never evaluated by expanding the sequence out to position \(x_n\); it is recovered from one residue modulo the cycle length and one lookup inside the stored tail.
The outer factor depends on \(d\) only through
$$q=\left\lfloor \frac{N}{d}\right\rfloor.$$
If this quotient is constant on an interval \(d\in [L,R]\), then the whole block contributes
$$\left(\sum_{d=L}^{R} w_d\right)\bigl(2\Phi(q)-1\bigr).$$
This is the second major compression. Instead of multiplying each \(w_d\) by its totient factor separately, the implementations maintain a running sum of the weights over each quotient block and apply the factor once. Since \(\left\lfloor N/d\right\rfloor\) has only \(O(\sqrt{N})\) distinct values, the expensive arithmetic attached to \(\Phi\) is heavily reused.
The first terms of the orbit are
$$x_1=1,\qquad x_2=2,\qquad x_3=3,\qquad x_4=10,\qquad x_5=731,\qquad x_{10}=24{,}881{,}905.$$
Hence the first gcd-weights are
$$w_1=1,\qquad w_2=2,\qquad w_3=3,\qquad w_4=x_{10}=24{,}881{,}905.$$
For \(N=4\), the relevant summatory-totient values are
$$\Phi(1)=1,\qquad \Phi(2)=2,\qquad \Phi(4)=6,$$
so
$$2\Phi(4)-1=11,\qquad 2\Phi(2)-1=3,\qquad 2\Phi(1)-1=1.$$
The divisor formula gives
$$G(4)=w_1\cdot 11+w_2\cdot 3+w_3\cdot 1+w_4\cdot 1.$$
Numerically,
$$G(4)=1\cdot 11+2\cdot 3+3\cdot 1+24{,}881{,}905\cdot 1=24{,}881{,}925.$$
This tiny example already shows the block structure: \(\lfloor 4/3\rfloor=\lfloor 4/4\rfloor=1\), so the last two divisor classes share the same factor.
The C++, Python, and Java implementations all follow the same mathematical plan. They never iterate over all \(N^2\) pairs. Instead, they combine orbit preprocessing, a fast summatory-totient routine, and a quotient-block sweep over the divisor index \(d\).
The implementation first applies Floyd's tortoise-and-hare method to the cubic recurrence and recovers the exact preperiod length \(53\) and cycle length \(33705\). It then stores the actual orbit values through the end of the useful periodic window, so every later access to \(x_k\) inside the tail becomes an array lookup. In parallel with the main divisor sweep, it also advances the same recurrence modulo \(33705\) so that the residue class of the inner index is always available.
For smaller arguments, the implementation builds a prefix table of \(\varphi(n)\) up to \(10^6\) with a linear sieve. For larger arguments, it uses the quotient-block recurrence for \(\Phi(n)\) together with memoization. This hybrid is important: values below the cutoff are answered in constant time from the prefix table, while larger values reuse the same cached quotient results again and again.
The divisor range \(1\le d\le N\) is partitioned into maximal intervals on which \(\left\lfloor N/d\right\rfloor\) is constant. While sweeping those intervals from left to right, the implementation generates the next weights \(w_d\), updates a running prefix sum of weights modulo \(M\), extracts the segment sum for the current block, computes the factor \(2\Phi(q)-1\), and adds the block contribution. This is why the final loop is linear in \(N\) rather than quadratic in \(N\).
Let \(L=\min(N,10^6)\). Cycle detection and orbit precomputation take \(O(53+33705)\) time and \(O(33705)\) memory, which are negligible compared with the main sweep. Building the totient prefix table costs \(O(L)\) time and \(O(L)\) memory because the sieve is linear. The memoized summatory-totient recurrence visits only \(O(\sqrt{N})\) distinct quotient values, and the quotient partition itself also has \(O(\sqrt{N})\) blocks.
However, the algorithm still generates every weight \(w_d\) for \(1\le d\le N\), so the dominant running time is \(O(N)\). Overall, the implementations run in \(O(N+L+\sqrt{N})\) time and use \(O(L+\sqrt{N}+33705)\) memory.
Definiere
$$M=123456789,\qquad x_1=1,\qquad x_{n+1}\equiv (x_n-1)^3+2 \pmod{M},\qquad 0\le x_n<M.$$
Gesucht ist die ggT-gewichtete Doppelsumme
$$G(N)=\sum_{1\le i,j\le N} w_{\gcd(i,j)} \pmod{M},\qquad w_n=x_{x_n},$$
wobei in der eigentlichen Aufgabe \(N=10^8\) ist. Eine direkte Auswertung ist ausgeschlossen: Die Außensumme läuft über \(N^2\) geordnete Paare, und jedes Gewicht \(w_n\) ist selbst wieder ein verschachtelter Zugriff auf dieselbe kubische Folge. Die Lösung funktioniert nur, weil der ggT-Teil arithmetisch gezählt werden kann und die Rekursion schließlich periodisch wird.
Die Herleitung besteht aus zwei Bausteinen. Zuerst wird die ggT-gewichtete Doppelsumme nach möglichen ggT-Werten umgeordnet und auf die summatorische Phi-Funktion zurückgeführt. Danach wird gezeigt, dass die kubische Bahn eine kurze Vorperiode und einen gut handhabbaren Zyklus besitzt, sodass das verschachtelte Gewicht \(w_n=x_{x_n}\) auf ein Nachschlagen innerhalb dieses periodischen Abschnitts reduziert werden kann.
Für ein festes geordnetes Paar \((i,j)\) schreiben wir
$$i=d\,a,\qquad j=d\,b,\qquad d=\gcd(i,j),\qquad \gcd(a,b)=1.$$
Ist \(d\) fest, dann muss das reduzierte Paar \((a,b)\) die Schranken
$$1\le a,b\le \left\lfloor \frac{N}{d}\right\rfloor$$
erfüllen. Definiert man
$$R(q)=\#\{(a,b)\in [1,q]^2:\gcd(a,b)=1\},$$
dann tragen alle geordneten Paare mit ggT genau \(d\) dasselbe Gewicht \(w_d\) bei. Also gilt
$$G(N)=\sum_{d=1}^{N} w_d\,R\!\left(\left\lfloor \frac{N}{d}\right\rfloor\right).$$
Das ist die erste entscheidende Vereinfachung: Aus der zweidimensionalen Summe wird eine eindimensionale Summe über ggT-Klassen.
\(R(q)\) zählt geordnete teilerfremde Paare im Quadrat \(q\times q\). Für jeden Nenner \(m\ge 2\) gibt es genau \(\varphi(m)\) reduzierte Zähler \(1\le a<m\). Da die Paare geordnet sind, erscheinen sowohl \((a,m)\) als auch \((m,a)\), und auf der Diagonalen kommt noch \((1,1)\) hinzu. Deshalb
$$R(q)=1+2\sum_{m=2}^{q}\varphi(m)=2\Phi(q)-1,$$
wobei
$$\Phi(q)=\sum_{m=1}^{q}\varphi(m).$$
Einsetzen in die ggT-Zerlegung liefert genau die Formel, die in den Implementierungen verwendet wird:
$$\boxed{G(N)=\sum_{d=1}^{N} w_d\left(2\Phi\!\left(\left\lfloor \frac{N}{d}\right\rfloor\right)-1\right)\pmod{M}.}$$
Der verbleibende arithmetische Engpass ist die summatorische Phi-Funktion \(\Phi(n)\). Die grundlegende Identität lautet
$$\sum_{k=1}^{n}\Phi\!\left(\left\lfloor \frac{n}{k}\right\rfloor\right)=\sum_{m=1}^{n}\varphi(m)\left\lfloor \frac{n}{m}\right\rfloor=\sum_{t=1}^{n}\sum_{m\mid t}\varphi(m)=\sum_{t=1}^{n} t=\frac{n(n+1)}{2}.$$
Isoliert man den Term \(k=1\), erhält man
$$\Phi(n)=\frac{n(n+1)}{2}-\sum_{k=2}^{n}\Phi\!\left(\left\lfloor \frac{n}{k}\right\rfloor\right).$$
Der Quotient \(\left\lfloor n/k\right\rfloor\) bleibt auf ganzen Intervallen konstant. Ist \(q=\left\lfloor n/\ell\right\rfloor\), dann ist der letzte Index mit demselben Quotienten
$$r=\left\lfloor \frac{n}{q}\right\rfloor=\left\lfloor \frac{n}{\lfloor n/\ell\rfloor}\right\rfloor.$$
Damit folgt die blockweise Form
$$\Phi(n)=\frac{n(n+1)}{2}-\sum_{\ell=2}^{n}(r-\ell+1)\,\Phi\!\left(\left\lfloor \frac{n}{\ell}\right\rfloor\right).$$
Genau dadurch werden große Werte von \(\Phi\) praktikabel: Die Rekursion besucht nur die verschiedenen Quotientenwerte, nicht alle ganzen Zahlen einzeln.
Die Rekursion
$$f(t)\equiv (t-1)^3+2 \pmod{M}$$
wirkt auf einem endlichen Zustandsraum, daher muss die Bahn ab \(x_1=1\) irgendwann zyklisch werden. Floyds Zyklensuche liefert für die tatsächliche Bahn eine Vorperiode der Länge \(53\) und einen Zyklus der Länge \(33705\). Der periodische Abschnitt beginnt also bei \(x_{54}\), und für jedes \(r\ge 54\) gilt
$$x_r=x_{\,54+((r-54)\bmod 33705)}.$$
Deshalb berechnen die Implementierungen die Bahn bis \(x_{33759}\) vor. Das deckt den ganzen nichtperiodischen Anfang und alle später benötigten periodischen Repräsentanten ab.
Der verschachtelte Index ist der heikle Teil. Weil dieselbe kubische Polynomabbildung die Folge definiert, verträgt sich die Reduktion modulo der Zykluslänge \(33705\) mit der Iteration. Man definiert die Begleitbahn
$$y_1=1,\qquad y_{n+1}\equiv (y_n-1)^3+2 \pmod{33705}.$$
Dann gilt für alle \(n\): \(y_n\equiv x_n\pmod{33705}\). Auf der tatsächlichen Bahn müssen nur die ersten vier inneren Indizes noch im Vorperiodenbereich direkt behandelt werden:
$$x_1=1,\qquad x_2=2,\qquad x_3=3,\qquad x_4=10.$$
Nach diesen Ausnahmefällen wird der innere Index durch
$$I_n=54+((y_n-54)\bmod 33705),\qquad w_n=x_{I_n}$$
zurück in den periodischen Abschnitt gehoben. Damit wird das scheinbar riesige Objekt \(x_{x_n}\) nie durch explizites Ausrollen bis zur Position \(x_n\) berechnet, sondern aus einem Rest modulo der Zykluslänge und einem Tabellenzugriff rekonstruiert.
Der äußere Faktor hängt von \(d\) nur über
$$q=\left\lfloor \frac{N}{d}\right\rfloor$$
ab. Ist dieser Quotient auf einem Intervall \(d\in [L,R]\) konstant, dann ist der gesamte Beitrag des Blocks
$$\left(\sum_{d=L}^{R} w_d\right)\bigl(2\Phi(q)-1\bigr).$$
Das ist die zweite große Kompression. Anstatt jedes \(w_d\) einzeln mit seinem Totientenfaktor zu multiplizieren, halten die Implementierungen eine laufende Summe der Gewichte innerhalb des Quotientenblocks und wenden den Faktor nur einmal an. Weil \(\left\lfloor N/d\right\rfloor\) nur \(O(\sqrt{N})\) verschiedene Werte annimmt, wird die teure \(\Phi\)-Arithmetik stark wiederverwendet.
Die ersten Folgenglieder sind
$$x_1=1,\qquad x_2=2,\qquad x_3=3,\qquad x_4=10,\qquad x_5=731,\qquad x_{10}=24{,}881{,}905.$$
Daraus folgen die ersten ggT-Gewichte
$$w_1=1,\qquad w_2=2,\qquad w_3=3,\qquad w_4=x_{10}=24{,}881{,}905.$$
Für \(N=4\) braucht man
$$\Phi(1)=1,\qquad \Phi(2)=2,\qquad \Phi(4)=6,$$
also
$$2\Phi(4)-1=11,\qquad 2\Phi(2)-1=3,\qquad 2\Phi(1)-1=1.$$
Die Divisorformel liefert
$$G(4)=w_1\cdot 11+w_2\cdot 3+w_3\cdot 1+w_4\cdot 1.$$
Numerisch ergibt das
$$G(4)=1\cdot 11+2\cdot 3+3\cdot 1+24{,}881{,}905\cdot 1=24{,}881{,}925.$$
Schon dieses kleine Beispiel zeigt die Blockstruktur: \(\lfloor 4/3\rfloor=\lfloor 4/4\rfloor=1\), also teilen die letzten beiden Divisorklassen denselben Faktor.
Die Implementierungen in C++, Python und Java folgen alle demselben mathematischen Plan. Sie iterieren niemals über alle \(N^2\) Paare. Stattdessen kombinieren sie eine Vorverarbeitung der Bahn, eine schnelle Routine für die summatorische Phi-Funktion und einen Quotientenblock-Durchlauf über den Divisorindex \(d\).
Zuerst bestimmt die Implementierung mit Floyds Schildkröte-und-Hase-Methode die exakte Vorperiodenlänge \(53\) und die Zykluslänge \(33705\) der kubischen Rekursion. Danach werden die tatsächlichen Folgenglieder bis zum Ende des nützlichen periodischen Fensters gespeichert, sodass jeder spätere Zugriff auf \(x_k\) im periodischen Abschnitt nur noch ein Array-Lookup ist. Parallel dazu wird während des Hauptdurchlaufs dieselbe Rekursion modulo \(33705\) fortgesetzt, damit die Restklasse des inneren Indexes stets verfügbar ist.
Für kleinere Argumente baut die Implementierung mit einem linearen Sieb eine Präfixtabelle von \(\varphi(n)\) bis \(10^6\) auf. Für größere Argumente benutzt sie die quotientenblockierte Rekursion für \(\Phi(n)\) zusammen mit Memoisierung. Diese Kombination ist entscheidend: Werte unterhalb der Grenze werden in konstanter Zeit aus der Präfixtabelle beantwortet, und größere Werte profitieren davon, dass dieselben Quotientenergebnisse immer wieder aus dem Cache gelesen werden.
Der Bereich \(1\le d\le N\) wird in maximale Intervalle zerlegt, auf denen \(\left\lfloor N/d\right\rfloor\) konstant bleibt. Beim Durchlauf dieser Intervalle erzeugt die Implementierung die nächsten Gewichte \(w_d\), aktualisiert eine laufende Präfixsumme der Gewichte modulo \(M\), gewinnt daraus die Segmentsumme des aktuellen Blocks, berechnet den Faktor \(2\Phi(q)-1\) und addiert den Blockbeitrag. Genau deshalb ist die Endschleife linear in \(N\) und nicht quadratisch.
Sei \(L=\min(N,10^6)\). Zyklensuche und Bahnvorberechnung benötigen \(O(53+33705)\) Zeit und \(O(33705)\) Speicher, was gegenüber dem Hauptdurchlauf klein ist. Das Erzeugen der Totienten-Präfixtabelle kostet wegen des linearen Siebs \(O(L)\) Zeit und \(O(L)\) Speicher. Die memoisierten Werte der summatorischen Phi-Funktion betreffen nur \(O(\sqrt{N})\) verschiedene Quotienten, und auch die Quotientenzerlegung selbst hat nur \(O(\sqrt{N})\) Blöcke.
Trotzdem wird jedes Gewicht \(w_d\) für \(1\le d\le N\) explizit erzeugt, also dominiert insgesamt \(O(N)\) Laufzeit. Damit ergibt sich \(O(N+L+\sqrt{N})\) Zeit und \(O(L+\sqrt{N}+33705)\) Speicher.
Şöyle tanımlayalım:
$$M=123456789,\qquad x_1=1,\qquad x_{n+1}\equiv (x_n-1)^3+2 \pmod{M},\qquad 0\le x_n<M.$$
Aranan büyüklük, gcd ile ağırlıklandırılmış şu çift toplamdır:
$$G(N)=\sum_{1\le i,j\le N} w_{\gcd(i,j)} \pmod{M},\qquad w_n=x_{x_n},$$
ve gerçek Euler girdisinde \(N=10^8\) kullanılır. Doğrudan hesaplama imkansızdır: dış ifade \(N^2\) sıralı çift içerir, her \(w_n\) de aynı kübik dizinin içine ikinci bir kez bakmayı gerektirir. Çözümün çalışmasının nedeni, gcd tarafının aritmetik olarak sayılabilmesi ve yinelemenin sonunda periyodik hâle gelmesidir.
Türetim iki ana parçadan oluşur. Önce gcd-ağırlıklı çift toplam, olası gcd değerlerine göre yeniden düzenlenip toplam totient fonksiyonuna indirgenir. Ardından kübik yörüngenin kısa bir ön-kuyruğa ve yönetilebilir bir döneme sahip olduğu gösterilir; böylece iç içe tanımlı \(w_n=x_{x_n}\) ifadesi, tek bir periyodik kuyruk içindeki tablo erişimine dönüşür.
Sabit bir \((i,j)\) çifti için
$$i=d\,a,\qquad j=d\,b,\qquad d=\gcd(i,j),\qquad \gcd(a,b)=1$$
yazalım. \(d\) sabitlendiğinde indirgenmiş \((a,b)\) çifti
$$1\le a,b\le \left\lfloor \frac{N}{d}\right\rfloor$$
koşulunu sağlamalıdır. Eğer
$$R(q)=\#\{(a,b)\in [1,q]^2:\gcd(a,b)=1\}$$
dersek, gcd'si tam olarak \(d\) olan her sıralı çift aynı \(w_d\) ağırlığını getirir. Dolayısıyla
$$G(N)=\sum_{d=1}^{N} w_d\,R\!\left(\left\lfloor \frac{N}{d}\right\rfloor\right)$$
elde edilir. Bu, ilk büyük sadeleştirmedir: iki boyutlu toplam, gcd sınıfları üzerindeki tek boyutlu toplama iner.
\(R(q)\), \(q\times q\) kare içindeki sıralı aralarında asal çiftleri sayar. Her \(m\ge 2\) için \(1\le a<m\) ve \(\gcd(a,m)=1\) olan tam \(\varphi(m)\) tane pay vardır. Çiftler sıralı olduğu için hem \((a,m)\) hem de \((m,a)\) sayılır; ayrıca köşegenden \((1,1)\) gelir. Bu yüzden
$$R(q)=1+2\sum_{m=2}^{q}\varphi(m)=2\Phi(q)-1,$$
burada
$$\Phi(q)=\sum_{m=1}^{q}\varphi(m).$$
Bunu gcd ayrışımına yerleştirince uygulamaların kullandığı tam formül ortaya çıkar:
$$\boxed{G(N)=\sum_{d=1}^{N} w_d\left(2\Phi\!\left(\left\lfloor \frac{N}{d}\right\rfloor\right)-1\right)\pmod{M}.}$$
Kalan aritmetik darboğaz toplam totient \(\Phi(n)\) değerleridir. Kodun dayandığı temel özdeşlik şudur:
$$\sum_{k=1}^{n}\Phi\!\left(\left\lfloor \frac{n}{k}\right\rfloor\right)=\sum_{m=1}^{n}\varphi(m)\left\lfloor \frac{n}{m}\right\rfloor=\sum_{t=1}^{n}\sum_{m\mid t}\varphi(m)=\sum_{t=1}^{n} t=\frac{n(n+1)}{2}.$$
Buradan \(k=1\) terimini ayırırsak
$$\Phi(n)=\frac{n(n+1)}{2}-\sum_{k=2}^{n}\Phi\!\left(\left\lfloor \frac{n}{k}\right\rfloor\right)$$
elde edilir. \(\left\lfloor n/k\right\rfloor\) değeri belli aralıklarda sabit kalır. Eğer \(q=\left\lfloor n/\ell\right\rfloor\) ise aynı bölümü veren son indis
$$r=\left\lfloor \frac{n}{q}\right\rfloor=\left\lfloor \frac{n}{\lfloor n/\ell\rfloor}\right\rfloor$$
olur. Böylece bloklu biçim
$$\Phi(n)=\frac{n(n+1)}{2}-\sum_{\ell=2}^{n}(r-\ell+1)\,\Phi\!\left(\left\lfloor \frac{n}{\ell}\right\rfloor\right)$$
çıkar. Büyük \(\Phi\) değerlerinin pratikleşmesini sağlayan nokta budur; yineleme bütün sayıları değil, yalnızca farklı bölüm değerlerini dolaşır.
Şu dönüşüm
$$f(t)\equiv (t-1)^3+2 \pmod{M}$$
sonlu bir durum uzayında çalıştığı için \(x_1=1\) ile başlayan yörünge eninde sonunda döngüye girmek zorundadır. Floyd döngü tespiti, gerçek yörünge için ön-kuyruk uzunluğunu \(53\), dönem uzunluğunu ise \(33705\) olarak verir. Dolayısıyla periyodik kuyruk \(x_{54}\) ile başlar ve her \(r\ge 54\) için
$$x_r=x_{\,54+((r-54)\bmod 33705)}$$
geçerlidir. Bu yüzden uygulamalar \(x_{33759}\)'a kadar ön hesap yapar; bu aralık hem periyodik olmayan başlangıcı hem de daha sonra sorgulanan periyodik temsilcileri kapsar.
İç içe indeks, problemin hassas kısmıdır. Aynı kübik polinom diziyi tanımladığı için, durumun \(33705\) dönem uzunluğuna göre indirgenmesi yineleme ile uyumludur. Bu amaçla yardımcı yörünge
$$y_1=1,\qquad y_{n+1}\equiv (y_n-1)^3+2 \pmod{33705}$$
tanımlanır. Böylece her \(n\) için \(y_n\equiv x_n\pmod{33705}\) olur. Gerçek yörüngede periyodik olmayan ön-kuyruğa doğrudan bakılması gereken tek iç indisler ilk dörttür:
$$x_1=1,\qquad x_2=2,\qquad x_3=3,\qquad x_4=10.$$
Bu istisnalardan sonra iç indis şu kaldırma ile periyodik kuyruğa taşınır:
$$I_n=54+((y_n-54)\bmod 33705),\qquad w_n=x_{I_n}.$$
Böylece devasa görünen \(x_{x_n}\) ifadesi, \(x_n\) kadar uzağa kadar diziyi açmadan, dönem uzunluğu modundaki bir artık ve saklanan kuyruk içindeki bir tablo erişimiyle elde edilir.
Dış katsayı, \(d\)'ye yalnızca
$$q=\left\lfloor \frac{N}{d}\right\rfloor$$
üzerinden bağlıdır. Eğer bu bölüm \(d\in [L,R]\) aralığında sabitse, o bloğun toplam katkısı
$$\left(\sum_{d=L}^{R} w_d\right)\bigl(2\Phi(q)-1\bigr)$$
olur. Bu ikinci büyük sıkıştırmadır. Her \(w_d\) için totient çarpanı ayrı ayrı hesaplanmaz; uygulamalar bölüm bloğu içindeki ağırlıkların kısmi toplamını tutar ve katsayıyı bir kez uygular. \(\left\lfloor N/d\right\rfloor\) yalnızca \(O(\sqrt{N})\) farklı değer aldığı için, pahalı \(\Phi\) aritmetiği yoğun biçimde yeniden kullanılır.
Yörüngenin ilk terimleri
$$x_1=1,\qquad x_2=2,\qquad x_3=3,\qquad x_4=10,\qquad x_5=731,\qquad x_{10}=24{,}881{,}905$$
şeklindedir. Buna karşılık ilk gcd-ağırlıkları
$$w_1=1,\qquad w_2=2,\qquad w_3=3,\qquad w_4=x_{10}=24{,}881{,}905$$
olur. \(N=4\) için gerekli toplam totient değerleri
$$\Phi(1)=1,\qquad \Phi(2)=2,\qquad \Phi(4)=6$$
ve dolayısıyla
$$2\Phi(4)-1=11,\qquad 2\Phi(2)-1=3,\qquad 2\Phi(1)-1=1.$$
Asıl formül
$$G(4)=w_1\cdot 11+w_2\cdot 3+w_3\cdot 1+w_4\cdot 1$$
şeklini alır. Sayısal olarak
$$G(4)=1\cdot 11+2\cdot 3+3\cdot 1+24{,}881{,}905\cdot 1=24{,}881{,}925.$$
Bu küçük örnek bile blok yapısını gösterir: \(\lfloor 4/3\rfloor=\lfloor 4/4\rfloor=1\) olduğundan son iki bölen sınıfı aynı katsayıyı paylaşır.
C++, Python ve Java uygulamaları aynı matematiksel planı izler. Hiçbiri bütün \(N^2\) çifti gezmez. Bunun yerine yörünge ön işlemesi, hızlı bir toplam totient yordamı ve bölen indisi \(d\) üzerinde bölüm-bloklu bir tarama birleştirilir.
Önce Floyd'un kaplumbağa-tavşan yöntemiyle kübik yinelemenin tam ön-kuyruk uzunluğu \(53\) ve dönem uzunluğu \(33705\) bulunur. Sonra gerçek dizi değerleri, gerekli periyodik pencerenin sonuna kadar saklanır; böylece kuyruğun içindeki herhangi bir \(x_k\) erişimi doğrudan dizi okumasına dönüşür. Ana bölen taramasıyla birlikte aynı yineleme bir kez de \(33705\) modunda ilerletilir; böylece iç indisin artık sınıfı her an elde tutulur.
Daha küçük argümanlar için uygulama \(10^6\)'ya kadar lineer sieve ile \(\varphi(n)\) önek tablosu kurar. Daha büyük argümanlarda ise bloklu bölüm yinelemesi ve memoization ile \(\Phi(n)\) hesaplanır. Bu hibrit yaklaşım önemlidir: eşik altındaki değerler önek tablodan sabit sürede gelir, büyük değerler ise aynı bölüm sonuçlarını tekrar tekrar önbellekten kullanır.
\(1\le d\le N\) aralığı, \(\left\lfloor N/d\right\rfloor\) sabit kalan en büyük bloklara ayrılır. Kod bu blokları soldan sağa gezerken yeni \(w_d\) değerlerini üretir, ağırlıkların mod \(M\) altındaki kısmi toplamını günceller, o bloğun ağırlık toplamını fark alarak çıkarır, \(2\Phi(q)-1\) katsayısını hesaplar ve blok katkısını sonuca ekler. Son döngünün kuadratik değil lineer olmasının nedeni tam olarak budur.
\(L=\min(N,10^6)\) olsun. Döngü tespiti ve yörünge önhesabı \(O(53+33705)\) zaman ve \(O(33705)\) bellek kullanır; bunlar ana taramaya göre küçüktür. Totient önek tablosu lineer sieve sayesinde \(O(L)\) zamanda ve \(O(L)\) bellekte kurulur. Belleklemeli toplam totient yinelemesi yalnızca \(O(\sqrt{N})\) farklı bölüm değerine dokunur; bölüm bloklarının sayısı da \(O(\sqrt{N})\)'dir.
Buna rağmen \(1\le d\le N\) için her \(w_d\) değeri tek tek üretilir, dolayısıyla baskın çalışma süresi \(O(N)\)'dir. Toplamda yöntem \(O(N+L+\sqrt{N})\) zaman ve \(O(L+\sqrt{N}+33705)\) bellek kullanır.
Definimos
$$M=123456789,\qquad x_1=1,\qquad x_{n+1}\equiv (x_n-1)^3+2 \pmod{M},\qquad 0\le x_n<M.$$
La cantidad buscada es la suma doble ponderada por el gcd
$$G(N)=\sum_{1\le i,j\le N} w_{\gcd(i,j)} \pmod{M},\qquad w_n=x_{x_n},$$
con \(N=10^8\) en la entrada real del problema. Un cálculo directo es inviable: la parte exterior recorre \(N^2\) pares ordenados y, además, cada peso \(w_n\) exige una consulta anidada dentro de la misma recurrencia cúbica. La solución solo funciona porque la parte del gcd se puede contar aritméticamente y la órbita de la recurrencia termina siendo periódica.
La derivación tiene dos ideas centrales. Primero, la suma doble ponderada por el gcd se reorganiza por valores posibles del gcd y se reduce a la función sumatoria del totiente. Segundo, la órbita cúbica tiene un prefijo corto y un ciclo manejable, de modo que el peso anidado \(w_n=x_{x_n}\) se convierte en una consulta dentro de una única cola periódica.
Para un par ordenado fijo \((i,j)\), escribimos
$$i=d\,a,\qquad j=d\,b,\qquad d=\gcd(i,j),\qquad \gcd(a,b)=1.$$
Una vez fijado \(d\), el par reducido \((a,b)\) debe satisfacer
$$1\le a,b\le \left\lfloor \frac{N}{d}\right\rfloor.$$
Si definimos
$$R(q)=\#\{(a,b)\in [1,q]^2:\gcd(a,b)=1\},$$
entonces todos los pares ordenados cuyo gcd es exactamente \(d\) aportan el mismo peso \(w_d\). Por tanto,
$$G(N)=\sum_{d=1}^{N} w_d\,R\!\left(\left\lfloor \frac{N}{d}\right\rfloor\right).$$
Esta es la primera simplificación decisiva: la suma bidimensional se reemplaza por una suma unidimensional sobre clases de gcd.
\(R(q)\) cuenta pares ordenados coprimos dentro del cuadrado \(q\times q\). Para cada denominador \(m\ge 2\), existen exactamente \(\varphi(m)\) numeradores reducidos \(1\le a<m\). Como los pares están ordenados, aparecen tanto \((a,m)\) como \((m,a)\), y además la diagonal aporta \((1,1)\). Así,
$$R(q)=1+2\sum_{m=2}^{q}\varphi(m)=2\Phi(q)-1,$$
donde
$$\Phi(q)=\sum_{m=1}^{q}\varphi(m).$$
Sustituyendo en la descomposición por gcd, obtenemos la fórmula exacta usada por las implementaciones:
$$\boxed{G(N)=\sum_{d=1}^{N} w_d\left(2\Phi\!\left(\left\lfloor \frac{N}{d}\right\rfloor\right)-1\right)\pmod{M}.}$$
El cuello de botella aritmético restante es la función sumatoria del totiente \(\Phi(n)\). La identidad clave es
$$\sum_{k=1}^{n}\Phi\!\left(\left\lfloor \frac{n}{k}\right\rfloor\right)=\sum_{m=1}^{n}\varphi(m)\left\lfloor \frac{n}{m}\right\rfloor=\sum_{t=1}^{n}\sum_{m\mid t}\varphi(m)=\sum_{t=1}^{n} t=\frac{n(n+1)}{2}.$$
Separando el término \(k=1\), queda
$$\Phi(n)=\frac{n(n+1)}{2}-\sum_{k=2}^{n}\Phi\!\left(\left\lfloor \frac{n}{k}\right\rfloor\right).$$
El cociente \(\left\lfloor n/k\right\rfloor\) permanece constante en intervalos. Si \(q=\left\lfloor n/\ell\right\rfloor\), entonces el último índice con ese mismo cociente es
$$r=\left\lfloor \frac{n}{q}\right\rfloor=\left\lfloor \frac{n}{\lfloor n/\ell\rfloor}\right\rfloor.$$
Al agrupar cocientes iguales aparece la forma por bloques
$$\Phi(n)=\frac{n(n+1)}{2}-\sum_{\ell=2}^{n}(r-\ell+1)\,\Phi\!\left(\left\lfloor \frac{n}{\ell}\right\rfloor\right).$$
Eso es lo que vuelve manejables los valores grandes de \(\Phi\): la recurrencia visita solo los distintos cocientes, no todos los enteros uno por uno.
La transformación
$$f(t)\equiv (t-1)^3+2 \pmod{M}$$
actúa sobre un espacio finito de estados, así que la órbita que parte de \(x_1=1\) debe repetir finalmente. El algoritmo de Floyd, aplicado a la órbita real, encuentra un preperíodo de longitud \(53\) y un período de longitud \(33705\). Por tanto, la cola periódica empieza en \(x_{54}\), y para todo \(r\ge 54\) se cumple
$$x_r=x_{\,54+((r-54)\bmod 33705)}.$$
Por eso las implementaciones precalculan la órbita hasta \(x_{33759}\), suficiente para cubrir todo el prefijo no periódico y todos los representantes periódicos que luego se consultan.
El índice anidado es la parte delicada. Como la misma transformación cúbica define la sucesión, reducir el estado módulo la longitud del ciclo \(33705\) es compatible con iterar. Se define entonces la órbita auxiliar
$$y_1=1,\qquad y_{n+1}\equiv (y_n-1)^3+2 \pmod{33705}.$$
Así, \(y_n\equiv x_n\pmod{33705}\) para todo \(n\). En la órbita real, los únicos índices internos que todavía requieren el prefijo no periódico son los cuatro primeros:
$$x_1=1,\qquad x_2=2,\qquad x_3=3,\qquad x_4=10.$$
Después de esos casos excepcionales, el índice interno se eleva a la cola periódica mediante
$$I_n=54+((y_n-54)\bmod 33705),\qquad w_n=x_{I_n}.$$
En otras palabras, el objeto aparentemente enorme \(x_{x_n}\) nunca se obtiene expandiendo la sucesión hasta la posición \(x_n\); se recupera a partir de un residuo módulo la longitud del ciclo y una lectura dentro de la cola almacenada.
El factor exterior depende de \(d\) solo a través de
$$q=\left\lfloor \frac{N}{d}\right\rfloor.$$
Si este cociente es constante en un intervalo \(d\in [L,R]\), entonces el bloque completo aporta
$$\left(\sum_{d=L}^{R} w_d\right)\bigl(2\Phi(q)-1\bigr).$$
Esta es la segunda gran compresión. En lugar de multiplicar cada \(w_d\) por su factor de totiente por separado, las implementaciones mantienen una suma acumulada de pesos dentro de cada bloque y aplican el factor una sola vez. Como \(\left\lfloor N/d\right\rfloor\) solo toma \(O(\sqrt{N})\) valores distintos, la parte cara asociada a \(\Phi\) se reutiliza intensamente.
Los primeros términos de la órbita son
$$x_1=1,\qquad x_2=2,\qquad x_3=3,\qquad x_4=10,\qquad x_5=731,\qquad x_{10}=24{,}881{,}905.$$
Por tanto, los primeros pesos por gcd son
$$w_1=1,\qquad w_2=2,\qquad w_3=3,\qquad w_4=x_{10}=24{,}881{,}905.$$
Para \(N=4\), se necesitan
$$\Phi(1)=1,\qquad \Phi(2)=2,\qquad \Phi(4)=6,$$
de modo que
$$2\Phi(4)-1=11,\qquad 2\Phi(2)-1=3,\qquad 2\Phi(1)-1=1.$$
La fórmula por divisores da
$$G(4)=w_1\cdot 11+w_2\cdot 3+w_3\cdot 1+w_4\cdot 1.$$
Numéricamente,
$$G(4)=1\cdot 11+2\cdot 3+3\cdot 1+24{,}881{,}905\cdot 1=24{,}881{,}925.$$
Incluso este ejemplo pequeño ya muestra la estructura por bloques: \(\lfloor 4/3\rfloor=\lfloor 4/4\rfloor=1\), así que las dos últimas clases de divisores comparten el mismo factor.
Las implementaciones en C++, Python y Java siguen el mismo plan matemático. Ninguna recorre los \(N^2\) pares. En su lugar combinan un preprocesamiento de la órbita, una rutina rápida para la función sumatoria del totiente y un barrido por bloques de cociente sobre el índice divisor \(d\).
Primero, la implementación aplica el método de la tortuga y la liebre de Floyd a la recurrencia cúbica y obtiene exactamente el preperíodo \(53\) y el período \(33705\). Después almacena los valores reales de la órbita hasta el final de la ventana periódica útil, de modo que cualquier acceso posterior a \(x_k\) dentro de la cola se convierte en una lectura de arreglo. En paralelo con el barrido principal por divisores, la misma recurrencia avanza módulo \(33705\), dejando siempre disponible la clase residual del índice interno.
Para argumentos pequeños, la implementación construye una tabla prefijo de \(\varphi(n)\) hasta \(10^6\) mediante un tamiz lineal. Para argumentos grandes, usa la recurrencia por bloques de cociente de \(\Phi(n)\) junto con memoización. Este enfoque híbrido es importante: los valores por debajo del umbral se responden en tiempo constante desde la tabla prefijo, mientras que los grandes reutilizan una y otra vez los mismos resultados almacenados en caché.
El rango \(1\le d\le N\) se divide en intervalos máximos donde \(\left\lfloor N/d\right\rfloor\) permanece constante. Mientras recorre esos intervalos de izquierda a derecha, la implementación genera los siguientes pesos \(w_d\), actualiza una suma prefijo de pesos módulo \(M\), extrae la suma del segmento del bloque actual, calcula el factor \(2\Phi(q)-1\) y añade la contribución del bloque. Por eso el bucle final es lineal en \(N\) y no cuadrático.
Sea \(L=\min(N,10^6)\). La detección del ciclo y el precálculo de la órbita cuestan \(O(53+33705)\) tiempo y \(O(33705)\) memoria, cantidades pequeñas frente al barrido principal. Construir la tabla prefijo del totiente cuesta \(O(L)\) tiempo y \(O(L)\) memoria porque el tamiz es lineal. La recurrencia memoizada de la función sumatoria del totiente visita solo \(O(\sqrt{N})\) cocientes distintos, y la propia partición en bloques de cociente también tiene \(O(\sqrt{N})\) bloques.
Aun así, el algoritmo sigue generando todos los pesos \(w_d\) para \(1\le d\le N\), de modo que el tiempo dominante es \(O(N)\). En conjunto, las implementaciones usan \(O(N+L+\sqrt{N})\) tiempo y \(O(L+\sqrt{N}+33705)\) memoria.
定义
$$M=123456789,\qquad x_1=1,\qquad x_{n+1}\equiv (x_n-1)^3+2 \pmod{M},\qquad 0\le x_n<M.$$
目标量是按 gcd 加权的二重和
$$G(N)=\sum_{1\le i,j\le N} w_{\gcd(i,j)} \pmod{M},\qquad w_n=x_{x_n},$$
在正式输入中 \(N=10^8\)。直接计算完全不可行：外层要处理 \(N^2\) 个有序数对，而每个权重 \(w_n\) 又不是普通的数组读取，而是对同一条三次递推序列进行一次嵌套索引。这个问题之所以能做，是因为 gcd 这一层可以用数论计数重写，而递推轨道本身最终会进入周期。
整个推导由两个核心简化组成。第一步，把 gcd 加权的双重求和按 gcd 的取值分层，化为与欧拉函数前缀和有关的一维求和。第二步，证明这条三次递推轨道只有很短的前导段和一个可控的周期，从而把嵌套权重 \(w_n=x_{x_n}\) 变成周期尾部中的一次查表。
对任意固定的有序对 \((i,j)\)，写成
$$i=d\,a,\qquad j=d\,b,\qquad d=\gcd(i,j),\qquad \gcd(a,b)=1.$$
一旦 \(d\) 固定，约去公因子后的 \((a,b)\) 必须满足
$$1\le a,b\le \left\lfloor \frac{N}{d}\right\rfloor.$$
记
$$R(q)=\#\{(a,b)\in [1,q]^2:\gcd(a,b)=1\},$$
即边长为 \(q\) 的方格里有多少个有序互素对。那么所有 gcd 恰好等于 \(d\) 的项都带来同一个权重 \(w_d\)，于是
$$G(N)=\sum_{d=1}^{N} w_d\,R\!\left(\left\lfloor \frac{N}{d}\right\rfloor\right).$$
这是第一个关键化简：原本二维的求和，被压缩成按照 gcd 类别进行的一维求和。
\(R(q)\) 统计的是 \(q\times q\) 方格中的有序互素对。对每个 \(m\ge 2\)，满足 \(1\le a<m\) 且 \(\gcd(a,m)=1\) 的 \(a\) 恰有 \(\varphi(m)\) 个。由于这里统计的是有序对，所以 \((a,m)\) 与 \((m,a)\) 都要算上，再加上对角线上的 \((1,1)\)，得到
$$R(q)=1+2\sum_{m=2}^{q}\varphi(m)=2\Phi(q)-1,$$
其中
$$\Phi(q)=\sum_{m=1}^{q}\varphi(m).$$
代回 gcd 分层公式后，就得到实现真正计算的主公式：
$$\boxed{G(N)=\sum_{d=1}^{N} w_d\left(2\Phi\!\left(\left\lfloor \frac{N}{d}\right\rfloor\right)-1\right)\pmod{M}.}$$
剩下的数论瓶颈是如何高效计算欧拉函数前缀和 \(\Phi(n)\)。程序使用的关键恒等式是
$$\sum_{k=1}^{n}\Phi\!\left(\left\lfloor \frac{n}{k}\right\rfloor\right)=\sum_{m=1}^{n}\varphi(m)\left\lfloor \frac{n}{m}\right\rfloor=\sum_{t=1}^{n}\sum_{m\mid t}\varphi(m)=\sum_{t=1}^{n} t=\frac{n(n+1)}{2}.$$
把 \(k=1\) 这一项单独拿出来，就有
$$\Phi(n)=\frac{n(n+1)}{2}-\sum_{k=2}^{n}\Phi\!\left(\left\lfloor \frac{n}{k}\right\rfloor\right).$$
\(\left\lfloor n/k\right\rfloor\) 会在一整段区间上保持不变。若 \(q=\left\lfloor n/\ell\right\rfloor\)，那么具有同样商值的最后一个位置是
$$r=\left\lfloor \frac{n}{q}\right\rfloor=\left\lfloor \frac{n}{\lfloor n/\ell\rfloor}\right\rfloor.$$
于是递推可以写成分块形式
$$\Phi(n)=\frac{n(n+1)}{2}-\sum_{\ell=2}^{n}(r-\ell+1)\,\Phi\!\left(\left\lfloor \frac{n}{\ell}\right\rfloor\right).$$
这一步非常重要，因为它说明程序不需要逐个递归所有整数，而只需要处理那些真正出现过的不同商值。
递推映射
$$f(t)\equiv (t-1)^3+2 \pmod{M}$$
作用在有限状态空间上，所以从 \(x_1=1\) 出发的轨道必然最终进入循环。对真实轨道应用 Floyd 判圈，可得到前导长度为 \(53\)，周期长度为 \(33705\)。因此从 \(x_{54}\) 开始进入周期尾部，并且对所有 \(r\ge 54\)，都有
$$x_r=x_{\,54+((r-54)\bmod 33705)}.$$
这就是为什么实现只需要预先保存到 \(x_{33759}\) 为止：这个范围足以覆盖非周期前缀以及后续所有会被查询到的周期代表位置。
嵌套索引是本题最微妙的部分。由于定义序列的还是同一个整数系数三次多项式，所以先对状态取模 \(33705\) 再迭代，与先迭代再取模 \(33705\) 是相容的。于是定义辅助轨道
$$y_1=1,\qquad y_{n+1}\equiv (y_n-1)^3+2 \pmod{33705}.$$
这样对每个 \(n\) 都有 \(y_n\equiv x_n\pmod{33705}\)。在真实轨道上，只有最开始四个内层索引还落在周期开始之前，需要直接处理：
$$x_1=1,\qquad x_2=2,\qquad x_3=3,\qquad x_4=10.$$
过了这些例外之后，内层索引就可以通过
$$I_n=54+((y_n-54)\bmod 33705),\qquad w_n=x_{I_n}$$
抬升回周期尾部。也就是说，程序并不会真的把序列展开到第 \(x_n\) 项；它只需要知道 \(x_n\) 在周期长度模意义下的位置，然后在保存好的尾部中查一次表即可。
外层因子对 \(d\) 的依赖只通过
$$q=\left\lfloor \frac{N}{d}\right\rfloor$$
体现。如果这个商在区间 \(d\in [L,R]\) 上保持不变，那么整个区间的贡献就是
$$\left(\sum_{d=L}^{R} w_d\right)\bigl(2\Phi(q)-1\bigr).$$
这是第二个关键压缩。程序不再对每个 \(d\) 单独计算 totient 因子，而是维护当前商块中所有权重的和，然后统一乘一次 \(2\Phi(q)-1\)。由于 \(\left\lfloor N/d\right\rfloor\) 只有 \(O(\sqrt{N})\) 个不同取值，和 \(\Phi\) 相关的昂贵计算就能被大量复用。
轨道的前几项为
$$x_1=1,\qquad x_2=2,\qquad x_3=3,\qquad x_4=10,\qquad x_5=731,\qquad x_{10}=24{,}881{,}905.$$
因此最开始几个 gcd 权重是
$$w_1=1,\qquad w_2=2,\qquad w_3=3,\qquad w_4=x_{10}=24{,}881{,}905.$$
当 \(N=4\) 时，需要的欧拉函数前缀和为
$$\Phi(1)=1,\qquad \Phi(2)=2,\qquad \Phi(4)=6,$$
于是
$$2\Phi(4)-1=11,\qquad 2\Phi(2)-1=3,\qquad 2\Phi(1)-1=1.$$
主公式变成
$$G(4)=w_1\cdot 11+w_2\cdot 3+w_3\cdot 1+w_4\cdot 1.$$
代入数值后得到
$$G(4)=1\cdot 11+2\cdot 3+3\cdot 1+24{,}881{,}905\cdot 1=24{,}881{,}925.$$
这个小例子已经能看出分块思想：\(\lfloor 4/3\rfloor=\lfloor 4/4\rfloor=1\)，所以最后两个约数类共享同一个因子。
C++、Python 和 Java 三个实现都遵循同一套数学方案。它们完全不会去枚举全部 \(N^2\) 个数对，而是把轨道预处理、快速的 totient 前缀和求法，以及按商值分块的 \(d\) 扫描组合起来。
实现首先用 Floyd 判圈精确求出三次递推的前导长度 \(53\) 和周期长度 \(33705\)。然后把真实轨道的值预计算到所需的周期窗口末端，这样以后任何落在周期尾部中的 \(x_k\) 都能直接数组读取。同时，在主除数扫描过程中，程序还会同步推进一次模 \(33705\) 的递推，以便随时拿到内层索引的剩余类。
对较小的参数，实现用线性筛构造到 \(10^6\) 为止的 \(\varphi(n)\) 前缀和表。对更大的参数，则使用前面推导出的按商分块递推，并配合记忆化缓存。这个混合结构很关键：阈值以下的结果可以直接常数时间返回，阈值以上的结果则不断复用已经算过的商值。
区间 \(1\le d\le N\) 被拆成若干个 \(\left\lfloor N/d\right\rfloor\) 不变的极大区间。实现按顺序扫描这些区间，逐步生成新的 \(w_d\)，更新模 \(M\) 下的权重前缀和，取出当前块对应的权重和，再乘上 \(2\Phi(q)-1\) 并加入答案。这就是为什么最终主循环是线性的，而不是二次的。
令 \(L=\min(N,10^6)\)。判圈和轨道预计算需要 \(O(53+33705)\) 时间与 \(O(33705)\) 空间，相对于主扫描来说很小。由于采用线性筛，totient 前缀表的构建成本是 \(O(L)\) 时间和 \(O(L)\) 空间。记忆化后的 \(\Phi\) 递推只会访问 \(O(\sqrt{N})\) 个不同商值，而商值分块本身也只有 \(O(\sqrt{N})\) 段。
不过，算法仍然必须为每个 \(1\le d\le N\) 生成一次 \(w_d\)，因此主导时间复杂度仍是 \(O(N)\)。整体复杂度为 \(O(N+L+\sqrt{N})\) 时间，空间复杂度为 \(O(L+\sqrt{N}+33705)\)。
Положим
$$M=123456789,\qquad x_1=1,\qquad x_{n+1}\equiv (x_n-1)^3+2 \pmod{M},\qquad 0\le x_n<M.$$
Искомая величина — двойная сумма с весом по gcd:
$$G(N)=\sum_{1\le i,j\le N} w_{\gcd(i,j)} \pmod{M},\qquad w_n=x_{x_n},$$
где в основной задаче \(N=10^8\). Прямое вычисление нереально: внешняя часть содержит \(N^2\) упорядоченных пар, а каждый вес \(w_n\) сам является вложенным обращением к той же кубической последовательности. Решение возможно только потому, что gcd-структуру можно пересчитать средствами арифметики, а сама рекурсия в конце концов становится периодической.
Вывод держится на двух упрощениях. Сначала двойная сумма с весами по gcd перераскладывается по возможным значениям gcd и сводится к сумматорной функции Эйлера. Затем для кубической орбиты находится короткий предпериод и достаточно небольшой цикл, так что вложенный вес \(w_n=x_{x_n}\) превращается в обращение к одной сохранённой периодической части.
Для фиксированной упорядоченной пары \((i,j)\) запишем
$$i=d\,a,\qquad j=d\,b,\qquad d=\gcd(i,j),\qquad \gcd(a,b)=1.$$
После фиксации \(d\) сокращённая пара \((a,b)\) должна удовлетворять
$$1\le a,b\le \left\lfloor \frac{N}{d}\right\rfloor.$$
Обозначим
$$R(q)=\#\{(a,b)\in [1,q]^2:\gcd(a,b)=1\}.$$
Тогда все упорядоченные пары, у которых gcd равен ровно \(d\), дают один и тот же вес \(w_d\), и потому
$$G(N)=\sum_{d=1}^{N} w_d\,R\!\left(\left\lfloor \frac{N}{d}\right\rfloor\right).$$
Это первое принципиальное упрощение: двумерная сумма заменяется одномерной суммой по gcd-классам.
\(R(q)\) считает упорядоченные взаимно простые пары внутри квадрата \(q\times q\). Для каждого знаменателя \(m\ge 2\) существует ровно \(\varphi(m)\) сокращённых числителей \(1\le a<m\). Так как пары упорядочены, учитываются и \((a,m)\), и \((m,a)\), а диагональ даёт ещё точку \((1,1)\). Следовательно,
$$R(q)=1+2\sum_{m=2}^{q}\varphi(m)=2\Phi(q)-1,$$
где
$$\Phi(q)=\sum_{m=1}^{q}\varphi(m).$$
Подстановка в разложение по gcd даёт точную формулу, используемую в реализациях:
$$\boxed{G(N)=\sum_{d=1}^{N} w_d\left(2\Phi\!\left(\left\lfloor \frac{N}{d}\right\rfloor\right)-1\right)\pmod{M}.}$$
Оставшееся арифметическое узкое место — сумматорная функция Эйлера \(\Phi(n)\). Ключевое тождество имеет вид
$$\sum_{k=1}^{n}\Phi\!\left(\left\lfloor \frac{n}{k}\right\rfloor\right)=\sum_{m=1}^{n}\varphi(m)\left\lfloor \frac{n}{m}\right\rfloor=\sum_{t=1}^{n}\sum_{m\mid t}\varphi(m)=\sum_{t=1}^{n} t=\frac{n(n+1)}{2}.$$
Выделяя член при \(k=1\), получаем
$$\Phi(n)=\frac{n(n+1)}{2}-\sum_{k=2}^{n}\Phi\!\left(\left\lfloor \frac{n}{k}\right\rfloor\right).$$
Значение \(\left\lfloor n/k\right\rfloor\) постоянно на целых интервалах. Если \(q=\left\lfloor n/\ell\right\rfloor\), то последний индекс с тем же частным равен
$$r=\left\lfloor \frac{n}{q}\right\rfloor=\left\lfloor \frac{n}{\lfloor n/\ell\rfloor}\right\rfloor.$$
После группировки одинаковых частных получается блочная форма
$$\Phi(n)=\frac{n(n+1)}{2}-\sum_{\ell=2}^{n}(r-\ell+1)\,\Phi\!\left(\left\lfloor \frac{n}{\ell}\right\rfloor\right).$$
Именно это делает большие значения \(\Phi\) вычислимыми на практике: рекурсия проходит только по различным значениям частного, а не по всем целым числам подряд.
Преобразование
$$f(t)\equiv (t-1)^3+2 \pmod{M}$$
действует на конечном пространстве состояний, поэтому орбита, начинающаяся в \(x_1=1\), неизбежно войдёт в цикл. Алгоритм Флойда, применённый к реальной орбите, даёт длину предпериода \(53\) и длину периода \(33705\). Значит, периодический хвост начинается с \(x_{54}\), и для любого \(r\ge 54\) выполняется
$$x_r=x_{\,54+((r-54)\bmod 33705)}.$$
Поэтому реализации предвычисляют орбиту до \(x_{33759}\): этого достаточно, чтобы покрыть непериодический префикс и все периодические представители, которые потом понадобятся.
Вложенный индекс — самая тонкая часть задачи. Поскольку ту же последовательность задаёт одно и то же кубическое полиномиальное отображение, редукция состояния по модулю длины цикла \(33705\) совместима с итерацией. Вводится вспомогательная орбита
$$y_1=1,\qquad y_{n+1}\equiv (y_n-1)^3+2 \pmod{33705}.$$
Тогда для каждого \(n\) справедливо \(y_n\equiv x_n\pmod{33705}\). На реальной орбите только первые четыре внутренних индекса ещё требуют прямого обращения к непериодическому префиксу:
$$x_1=1,\qquad x_2=2,\qquad x_3=3,\qquad x_4=10.$$
После этих исключений внутренний индекс переносится в хвост цикла формулой
$$I_n=54+((y_n-54)\bmod 33705),\qquad w_n=x_{I_n}.$$
То есть объект \(x_{x_n}\), который выглядит как обращение к очень далёкому месту последовательности, на самом деле восстанавливается по одному остатку по модулю длины цикла и одному чтению из сохранённого периодического хвоста.
Внешний множитель зависит от \(d\) только через
$$q=\left\lfloor \frac{N}{d}\right\rfloor.$$
Если этот частный постоянен на интервале \(d\in [L,R]\), то вклад всего блока равен
$$\left(\sum_{d=L}^{R} w_d\right)\bigl(2\Phi(q)-1\bigr).$$
Это второе большое сжатие. Вместо того чтобы по отдельности умножать каждое \(w_d\) на соответствующий totient-фактор, реализации поддерживают накопленную сумму весов внутри блока одинакового частного и применяют множитель один раз. Так как \(\left\lfloor N/d\right\rfloor\) принимает лишь \(O(\sqrt{N})\) различных значений, дорогая часть, связанная с \(\Phi\), используется многократно.
Первые члены орбиты таковы:
$$x_1=1,\qquad x_2=2,\qquad x_3=3,\qquad x_4=10,\qquad x_5=731,\qquad x_{10}=24{,}881{,}905.$$
Значит, первые gcd-веса равны
$$w_1=1,\qquad w_2=2,\qquad w_3=3,\qquad w_4=x_{10}=24{,}881{,}905.$$
Для \(N=4\) нужны значения
$$\Phi(1)=1,\qquad \Phi(2)=2,\qquad \Phi(4)=6,$$
следовательно,
$$2\Phi(4)-1=11,\qquad 2\Phi(2)-1=3,\qquad 2\Phi(1)-1=1.$$
Формула по делителям даёт
$$G(4)=w_1\cdot 11+w_2\cdot 3+w_3\cdot 1+w_4\cdot 1.$$
После подстановки получаем
$$G(4)=1\cdot 11+2\cdot 3+3\cdot 1+24{,}881{,}905\cdot 1=24{,}881{,}925.$$
Уже этот маленький пример показывает блочную структуру: \(\lfloor 4/3\rfloor=\lfloor 4/4\rfloor=1\), значит, два последних класса делителей получают один и тот же множитель.
Реализации на C++, Python и Java придерживаются одного и того же математического плана. Они вообще не перебирают все \(N^2\) пар. Вместо этого сочетаются предобработка орбиты, быстрая процедура для сумматорной функции Эйлера и блочный проход по индексам делителя \(d\).
Сначала реализация применяет метод черепахи и зайца Флойда к кубической рекурсии и точно находит длину предпериода \(53\) и длину цикла \(33705\). Затем реальные значения орбиты сохраняются до конца полезного периодического окна, так что любое позднее обращение к \(x_k\) внутри хвоста становится чтением из массива. Одновременно в основном проходе по делителям та же рекурсия продвигается по модулю \(33705\), чтобы всегда была известна остаточная позиция внутреннего индекса.
Для малых аргументов реализация строит префиксную таблицу \(\varphi(n)\) до \(10^6\) с помощью линейного решета. Для больших аргументов используется описанная выше рекурсия по блокам одинакового частного и мемоизация. Такой гибрид важен: значения ниже порога отвечаются за константное время из таблицы, а большие значения многократно используют одни и те же закэшированные частные.
Диапазон \(1\le d\le N\) разбивается на максимальные интервалы, на которых \(\left\lfloor N/d\right\rfloor\) постоянно. Во время прохода по этим интервалам реализация порождает очередные веса \(w_d\), обновляет префиксную сумму весов по модулю \(M\), извлекает сумму на текущем блоке, вычисляет множитель \(2\Phi(q)-1\) и добавляет вклад блока. Именно поэтому финальный цикл линейный по \(N\), а не квадратичный.
Пусть \(L=\min(N,10^6)\). Поиск цикла и предвычисление орбиты требуют \(O(53+33705)\) времени и \(O(33705)\) памяти, что мало по сравнению с основным проходом. Построение префиксной таблицы тотиентов стоит \(O(L)\) времени и \(O(L)\) памяти благодаря линейному решету. Мемоизированная рекурсия для сумматорной функции Эйлера посещает лишь \(O(\sqrt{N})\) различных частных, и сама блочная разбивка по частным тоже имеет размер \(O(\sqrt{N})\).
Тем не менее каждый вес \(w_d\) для \(1\le d\le N\) всё равно порождается явно, поэтому доминирующее время работы равно \(O(N)\). В итоге получаем \(O(N+L+\sqrt{N})\) по времени и \(O(L+\sqrt{N}+33705)\) по памяти.
لنعرّف
$$M=123456789,\qquad x_1=1,\qquad x_{n+1}\equiv (x_n-1)^3+2 \pmod{M},\qquad 0\le x_n<M.$$
والكمية المطلوبة هي المجموع المزدوج الموزون بحسب القاسم المشترك الأكبر:
$$G(N)=\sum_{1\le i,j\le N} w_{\gcd(i,j)} \pmod{M},\qquad w_n=x_{x_n},$$
وفي مسألة Project Euler الحقيقية يكون \(N=10^8\). الحساب المباشر مستحيل عملياً، لأن الحد الخارجي يتعامل مع \(N^2\) من الأزواج المرتبة، ولكل وزن \(w_n\) فهرسة متداخلة داخل المتتالية التكعيبية نفسها. سبب إمكان الحل هو أن جزء gcd يمكن إعادة صياغته بعدٍّ عددي، وأن مدار العلاقة التكرارية يصبح دورياً في النهاية.
يعتمد الاشتقاق على تبسيطين أساسيين. أولاً نعيد ترتيب المجموع المزدوج بحسب قيم gcd الممكنة ونختزله إلى الدالة التراكمية لتوتينت أويلر. ثانياً نثبت أن المدار التكعيبي يملك مقدمة قصيرة ودورة بطول يمكن التحكم فيه، وبذلك يتحول الوزن المتداخل \(w_n=x_{x_n}\) إلى قراءة واحدة من ذيل دوري محفوظ.
لكل زوج مرتب ثابت \((i,j)\) نكتب
$$i=d\,a,\qquad j=d\,b,\qquad d=\gcd(i,j),\qquad \gcd(a,b)=1.$$
بعد تثبيت \(d\)، يجب أن يحقق الزوج المختزل \((a,b)\)
$$1\le a,b\le \left\lfloor \frac{N}{d}\right\rfloor.$$
إذا عرّفنا
$$R(q)=\#\{(a,b)\in [1,q]^2:\gcd(a,b)=1\},$$
فإن كل الأزواج المرتبة التي يكون gcd لها مساوياً تماماً لـ \(d\) تعطي الوزن نفسه \(w_d\)، ولذلك
$$G(N)=\sum_{d=1}^{N} w_d\,R\!\left(\left\lfloor \frac{N}{d}\right\rfloor\right).$$
هذا هو التبسيط الحاسم الأول: المجموع ثنائي الأبعاد يتحول إلى مجموع أحادي الأبعاد على طبقات gcd.
الكمية \(R(q)\) تعد الأزواج المرتبة المتباينة أولياً داخل مربع \(q\times q\). لكل مقام \(m\ge 2\) يوجد بالضبط \(\varphi(m)\) من البسوط المختزلة \(1\le a<m\). وبما أن الأزواج مرتبة فإن \((a,m)\) و\((m,a)\) كلاهما يُحسب، ثم تضيف القطرية النقطة \((1,1)\). ومن ثم
$$R(q)=1+2\sum_{m=2}^{q}\varphi(m)=2\Phi(q)-1,$$
حيث
$$\Phi(q)=\sum_{m=1}^{q}\varphi(m).$$
بالتعويض في تفكيك gcd نحصل على الصيغة الدقيقة التي تعتمدها التطبيقات:
$$\boxed{G(N)=\sum_{d=1}^{N} w_d\left(2\Phi\!\left(\left\lfloor \frac{N}{d}\right\rfloor\right)-1\right)\pmod{M}.}$$
العنق الزجاجي العددي المتبقي هو الدالة التراكمية \(\Phi(n)\). والهوية الأساسية هي
$$\sum_{k=1}^{n}\Phi\!\left(\left\lfloor \frac{n}{k}\right\rfloor\right)=\sum_{m=1}^{n}\varphi(m)\left\lfloor \frac{n}{m}\right\rfloor=\sum_{t=1}^{n}\sum_{m\mid t}\varphi(m)=\sum_{t=1}^{n} t=\frac{n(n+1)}{2}.$$
وبفصل الحد الموافق لـ \(k=1\) نحصل على
$$\Phi(n)=\frac{n(n+1)}{2}-\sum_{k=2}^{n}\Phi\!\left(\left\lfloor \frac{n}{k}\right\rfloor\right).$$
القيمة \(\left\lfloor n/k\right\rfloor\) تبقى ثابتة على فترات كاملة. إذا كان \(q=\left\lfloor n/\ell\right\rfloor\)، فإن آخر فهرس يملك الخارج نفسه هو
$$r=\left\lfloor \frac{n}{q}\right\rfloor=\left\lfloor \frac{n}{\lfloor n/\ell\rfloor}\right\rfloor.$$
وبجمع الحدود ذات الخارج المتساوي نحصل على الصيغة الكتلية
$$\Phi(n)=\frac{n(n+1)}{2}-\sum_{\ell=2}^{n}(r-\ell+1)\,\Phi\!\left(\left\lfloor \frac{n}{\ell}\right\rfloor\right).$$
هذه هي الخطوة التي تجعل القيم الكبيرة لـ \(\Phi\) قابلة للحساب عملياً، لأن العلاقة لا تزور كل الأعداد على حدة بل تمر فقط على قيم الخارج المختلفة فعلاً.
التحويل
$$f(t)\equiv (t-1)^3+2 \pmod{M}$$
يعمل على فضاء حالات منتهٍ، لذا فإن المدار الذي يبدأ من \(x_1=1\) لا بد أن يدخل في دورة في نهاية المطاف. تطبيق خوارزمية Floyd على المدار الحقيقي يعطي طول مقدمة يساوي \(53\) وطول دورة يساوي \(33705\). لذلك يبدأ الذيل الدوري عند \(x_{54}\)، ولكل \(r\ge 54\) يكون
$$x_r=x_{\,54+((r-54)\bmod 33705)}.$$
ولهذا السبب تُحسب القيم مسبقاً حتى \(x_{33759}\)، وهو مدى يكفي لتغطية المقدمة غير الدورية وكل الممثلين الدوريين الذين ستحتاج إليهم العمليات اللاحقة.
الفهرسة المتداخلة هي الجزء الأدق في المسألة. وبما أن المتتالية يحددها التحويل التكعيبي نفسه، فإن اختزال الحالة بترديد طول الدورة \(33705\) ينسجم مع التكرار. لذلك نعرّف المدار المساعد
$$y_1=1,\qquad y_{n+1}\equiv (y_n-1)^3+2 \pmod{33705}.$$
وعندئذ يتحقق \(y_n\equiv x_n\pmod{33705}\) لكل \(n\). وفي المدار الحقيقي لا تحتاج المقدمة غير الدورية إلا إلى أول أربعة فهارس داخلية:
$$x_1=1,\qquad x_2=2,\qquad x_3=3,\qquad x_4=10.$$
بعد هذه الحالات الاستثنائية يرفع الفهرس الداخلي إلى الذيل الدوري بالصيغة
$$I_n=54+((y_n-54)\bmod 33705),\qquad w_n=x_{I_n}.$$
وهكذا فإن المقدار الكبير ظاهرياً \(x_{x_n}\) لا يُحسب بتمديد المتتالية حتى الموضع \(x_n\)، بل يستعاد من باقي واحد بترديد طول الدورة ومن قراءة واحدة داخل الذيل المخزن.
العامل الخارجي يعتمد على \(d\) فقط عبر
$$q=\left\lfloor \frac{N}{d}\right\rfloor.$$
فإذا كان هذا الخارج ثابتاً على الفترة \(d\in [L,R]\)، فإن مساهمة الكتلة كلها تساوي
$$\left(\sum_{d=L}^{R} w_d\right)\bigl(2\Phi(q)-1\bigr).$$
وهذا هو التبسيط الكبير الثاني. بدلاً من ضرب كل \(w_d\) في عامله المرتبط بالتوتينت على حدة، تحتفظ التطبيقات بمجموع جارٍ للأوزان داخل كل كتلة ثم تطبق العامل مرة واحدة فقط. وبما أن \(\left\lfloor N/d\right\rfloor\) يأخذ \(O(\sqrt{N})\) قيمة مختلفة فقط، فإن الحساب المكلف المرتبط بـ \(\Phi\) يعاد استخدامه كثيراً.
أوائل حدود المدار هي
$$x_1=1,\qquad x_2=2,\qquad x_3=3,\qquad x_4=10,\qquad x_5=731,\qquad x_{10}=24{,}881{,}905.$$
ومن ثم تكون أول أوزان gcd هي
$$w_1=1,\qquad w_2=2,\qquad w_3=3,\qquad w_4=x_{10}=24{,}881{,}905.$$
عند \(N=4\) نحتاج إلى
$$\Phi(1)=1,\qquad \Phi(2)=2,\qquad \Phi(4)=6,$$
ومن ثم
$$2\Phi(4)-1=11,\qquad 2\Phi(2)-1=3,\qquad 2\Phi(1)-1=1.$$
فتصبح الصيغة
$$G(4)=w_1\cdot 11+w_2\cdot 3+w_3\cdot 1+w_4\cdot 1.$$
وبالتعويض العددي نحصل على
$$G(4)=1\cdot 11+2\cdot 3+3\cdot 1+24{,}881{,}905\cdot 1=24{,}881{,}925.$$
حتى هذا المثال الصغير يكشف فكرة الكتل: إذ إن \(\lfloor 4/3\rfloor=\lfloor 4/4\rfloor=1\)، ولذلك تشترك آخر فئتين من القواسم في العامل نفسه.
تتبع تطبيقات C++ وPython وJava الخطة الرياضية نفسها. فهي لا تعد جميع الأزواج البالغ عددها \(N^2\) إطلاقاً، بل تجمع بين تمهيد المدار، وروتين سريع للدالة التراكمية لتوتينت، ومسح كتلي على فهرس القاسم \(d\).
تبدأ الشيفرة بتطبيق طريقة السلحفاة والأرنب لفلويد على العلاقة التكعيبية لاستخراج طول المقدمة \(53\) وطول الدورة \(33705\) بدقة. بعد ذلك تُخزن قيم المدار الحقيقي حتى نهاية النافذة الدورية المفيدة، بحيث تصبح أي قراءة لاحقة لـ \(x_k\) داخل الذيل مجرد وصول إلى مصفوفة. وبالتوازي مع المسح الرئيسي على القواسم، تُحدَّث العلاقة نفسها أيضاً بترديد \(33705\) حتى تبقى فئة باقي الفهرس الداخلي متاحة دائماً.
للوسائط الصغيرة تبني الشيفرة جدولاً سابقاً لقيم \(\varphi(n)\) حتى \(10^6\) باستعمال غربال خطي. وللوسائط الأكبر تستخدم العلاقة الكتلية السابقة مع الحفظ في الذاكرة. هذا الدمج مهم جداً: القيم الأصغر من الحد تعاد فوراً من الجدول، والقيم الأكبر تعيد استخدام النتائج نفسها المخزنة لمرات عديدة.
يُقسَّم المجال \(1\le d\le N\) إلى فترات عظمى يبقى فيها \(\left\lfloor N/d\right\rfloor\) ثابتاً. وأثناء المرور على هذه الفترات من اليسار إلى اليمين، تولد الشيفرة قيم \(w_d\) التالية، وتحدث مجموعاً سابقاً للأوزان modulo \(M\)، وتستخرج مجموع الكتلة الحالية، ثم تحسب العامل \(2\Phi(q)-1\) وتضيف مساهمة الكتلة. ولهذا يكون الحلقـة النهائية خطية في \(N\) وليست تربيعية.
لنضع \(L=\min(N,10^6)\). إن كشف الدورة وتمهيد المدار يكلفان \(O(53+33705)\) زمناً و\(O(33705)\) ذاكرة، وهو مقدار صغير مقارنة بالمسح الرئيسي. وبفضل الغربال الخطي، فإن بناء جدول التوتينت السابق يكلف \(O(L)\) زمناً و\(O(L)\) ذاكرة. أما العلاقة المحفوظة للدالة التراكمية لتوتينت فلا تزور إلا \(O(\sqrt{N})\) من قيم الخارج المختلفة، كما أن عدد كتل الخارج نفسه يساوي \(O(\sqrt{N})\).
ومع ذلك لا بد من توليد كل وزن \(w_d\) لكل \(1\le d\le N\)، ولذلك يبقى الزمن الغالب هو \(O(N)\). وبصورة كلية نحصل على \(O(N+L+\sqrt{N})\) زمناً و\(O(L+\sqrt{N}+33705)\) ذاكرة.