Start from the reduced fraction \(1/n\). At each step replace \(x/y\) by \((x+1)/(y-1)\) and reduce again. If the first time the denominator becomes 1 the numerator is \(f(n)\), the problem asks for
$$\sum_{k=1}^{L} f(k^3),\qquad L=2\cdot 10^6.$$
A direct simulation for every cube is far too slow, so the key is to rewrite \(f(n)\) using the prime factors of \(n+1\).
Write any reduced state as \(a/b\) with \(a,b>0\), and define its state sum by
$$s=a+b.$$
Initially we have \(a=1\), \(b=n\), hence \(s=n+1\).
From the reduced fraction \(a/b\), the unreduced next fraction is \((a+1)/(b-1)\). If
$$d=\gcd(a+1,b-1),$$
then after reduction
$$a'=\frac{a+1}{d},\qquad b'=\frac{b-1}{d}.$$
Since \(b=s-a\), we get
$$d=\gcd(a+1,s-a-1)=\gcd(a+1,s).$$
Therefore the next state sum is
$$s'=a'+b'=\frac{s}{\gcd(a+1,s)}.$$
So every reduction replaces the current sum \(s\) by a divisor of \(s\). Also, because \(a/b\) is reduced, \(\gcd(a,b)=1\), and thus \(\gcd(a,s)=1\).
Suppose a phase starts at the special state \(1/(s-1)\). Let \(p\) be the smallest prime factor of \(s\). For every integer \(t\) with \(2 \le t \lt p\), we have \(\gcd(t,s)=1\). Hence no reduction happens while the numerator runs through \(1,2,\dots,p-1\).
When the unreduced numerator becomes \(p\), the reduction factor is
$$d=\gcd(p,s)=p.$$
The fraction at that moment is
$$\frac{p}{s-p}\to \frac{1}{s/p-1}.$$
So one full phase transforms the state sum by
$$s\longmapsto \frac{s}{p},$$
which means: a phase strips off the smallest prime factor of the current sum and resets the numerator to 1.
Starting from \(s_0=n+1\), repeated phases remove the prime factors of \(n+1\) from smallest to largest. After all but the largest prime factor have been removed, the state sum becomes
$$P^+(n+1)=\operatorname{LPF}(n+1),$$
where \(\operatorname{LPF}\) denotes the largest prime factor.
If the current sum \(s\) is prime, then for every \(1\le a\le s-2\) we have \(\gcd(a+1,s)=1\), so no further reductions occur before the denominator reaches 1. Beginning from \(1/(s-1)\), the chain ends at \((s-1)/1\). Therefore
$$\boxed{f(n)=\operatorname{LPF}(n+1)-1.}$$
The iteration is
$$\frac{1}{20}\to\frac{2}{19}\to\frac{3}{18}=\frac{1}{6}\to\frac{2}{5}\to\frac{3}{4}\to\frac{4}{3}\to\frac{5}{2}\to\frac{6}{1}.$$
So \(f(20)=6\). Since \(20+1=21=3\cdot 7\), the closed form gives
$$\operatorname{LPF}(21)-1=7-1=6,$$
which matches the simulation exactly.
For the required sum, set \(n=k^3\):
$$f(k^3)=\operatorname{LPF}(k^3+1)-1.$$
Using the factorization of a sum of cubes,
$$k^3+1=(k+1)(k^2-k+1).$$
If we define
$$Q_k=k^2-k+1,$$
then
$$\operatorname{LPF}(k^3+1)=\max\bigl(\operatorname{LPF}(k+1),\operatorname{LPF}(Q_k)\bigr).$$
This remains true even when the two factors share a factor 3, because the largest prime factor of a product is still the maximum of the largest prime factors of its factors. Hence
$$\sum_{k=1}^{L}f(k^3)=\sum_{k=1}^{L}\left(\max\bigl(\operatorname{LPF}(k+1),\operatorname{LPF}(Q_k)\bigr)-1\right).$$
A standard largest-prime-factor sieve on \([2,L+1]\) fills an array \(\operatorname{LPF}(m)\) for every \(m\le L+1\). This immediately gives all values \(\operatorname{LPF}(k+1)\).
For each \(k\), initialize
$$\mathrm{rem}[k]=Q_k,\qquad \mathrm{lpfQ}[k]=1.$$
Now fix a prime \(p\). We need to know for which \(k\) the congruence
$$Q_k\equiv0\pmod p$$
holds. Multiplying by 4 gives
$$4k^2-4k+4\equiv0\pmod p\qquad\Longleftrightarrow\qquad (2k-1)^2\equiv -3\pmod p.$$
Therefore, for \(p>3\), roots exist exactly when \(-3\) is a quadratic residue modulo \(p\). If \(s^2\equiv -3\pmod p\), then
$$k\equiv \frac{1\pm s}{2}\pmod p.$$
The code finds the square root \(s\) with the Tonelli-Shanks algorithm. The prime \(p=3\) is special: the only class is \(k\equiv2\pmod 3\). The prime \(2\) never divides \(Q_k\), because \(k^2-k=k(k-1)\) is always even, so \(Q_k\) is always odd.
For every matching residue class, the code visits all corresponding \(k\), divides out all powers of \(p\) from \(\mathrm{rem}[k]\), and stores \(p\) as the current largest prime factor. After all primes \(p\le L\) are processed, any leftover \(\mathrm{rem}[k]>1\) must itself be prime and greater than \(L\); otherwise a composite remainder would have a prime divisor \(\le \sqrt{Q_k} \lt L\), which would already have been removed. That leftover is therefore the final large prime factor of \(Q_k\).
Here
$$k^3+1=65=5\cdot 13,$$
so
$$f(64)=\operatorname{LPF}(65)-1=13-1=12.$$
The decomposition gives \(k+1=5\) and \(Q_4=13\), so the larger prime factor is indeed \(13\).
The implementation has two sieve layers. First it fills lpf_small for all integers up to \(L+1\).
Then it stores every \(Q_k\) in the residual array rem and processes primes one by one. For each prime
it computes the admissible residue classes of \(k\), strips that prime from all matching residuals, and updates
lpf_q[k]. Finally it combines lpf_small[k+1] with lpf_q[k], subtracts 1, and
accumulates the sum in a 128-bit total. The C++ version checks \(f(20)=6\),
\(\sum_{k\le100} f(k^3)=118937\), and also matches a brute-force computation for small limits.
The sieve for \(\operatorname{LPF}(k+1)\) is near-linear in practice. The \(Q_k\)-sieve visits about \(L/p\) or \(2L/p\) positions per relevant prime, so its dominant work is also near-linear up to the usual \(\sum_{p\le L}1/p\) factor, i.e. \(O(L\log\log L)\) with extra polylogarithmic cost for modular square roots. Memory usage is \(O(L)\) for the factor tables and residual buffer.
Man startet mit dem vollständig gekürzten Bruch \(1/n\). In jedem Schritt wird \(x/y\) durch \((x+1)/(y-1)\) ersetzt und erneut gekürzt. Wenn beim ersten Auftreten des Nenners 1 der Zähler \(f(n)\) ist, soll berechnet werden
$$\sum_{k=1}^{L} f(k^3),\qquad L=2\cdot 10^6.$$
Eine direkte Simulation für jedes \(k^3\) wäre viel zu langsam; entscheidend ist daher eine geschlossene Formel für \(f(n)\) über die Primfaktoren von \(n+1\).
Schreibe jeden gekürzten Zustand als \(a/b\) mit \(a,b>0\) und definiere die Zustandssumme
$$s=a+b.$$
Am Anfang gilt \(a=1\), \(b=n\), also \(s=n+1\).
Aus dem gekürzten Bruch \(a/b\) wird zunächst der ungekürzte Bruch \((a+1)/(b-1)\). Setze
$$d=\gcd(a+1,b-1).$$
Nach dem Kürzen erhält man
$$a'=\frac{a+1}{d},\qquad b'=\frac{b-1}{d}.$$
Mit \(b=s-a\) folgt
$$d=\gcd(a+1,s-a-1)=\gcd(a+1,s).$$
Die neue Zustandssumme ist also
$$s'=a'+b'=\frac{s}{\gcd(a+1,s)}.$$
Jede Kürzung ersetzt die aktuelle Summe \(s\) somit durch einen Teiler von \(s\). Weil \(a/b\) bereits gekürzt ist, gilt außerdem \(\gcd(a,b)=1\) und damit \(\gcd(a,s)=1\).
Betrachte eine Phase, die im Zustand \(1/(s-1)\) beginnt. Sei \(p\) der kleinste Primfaktor von \(s\). Für alle ganzen Zahlen \(t\) mit \(2\le t\lt p\) gilt \(\gcd(t,s)=1\). Daher findet zunächst keine weitere Kürzung statt, während der Zähler die Werte \(1,2,\dots,p-1\) durchläuft.
Sobald der ungekürzte Zähler den Wert \(p\) erreicht, ist der Kürzungsfaktor
$$d=\gcd(p,s)=p.$$
In diesem Moment wird
$$\frac{p}{s-p}\to \frac{1}{s/p-1}.$$
Eine vollständige Phase bewirkt also
$$s\longmapsto \frac{s}{p},$$
das heißt: Der kleinste Primfaktor der aktuellen Summe wird entfernt und der Zähler springt wieder auf 1 zurück.
Beginnend mit \(s_0=n+1\) entfernen die aufeinanderfolgenden Phasen die Primfaktoren von \(n+1\) in aufsteigender Reihenfolge. Wenn alle bis auf den größten Primfaktor entfernt sind, bleibt
$$P^+(n+1)=\operatorname{LPF}(n+1),$$
wobei \(\operatorname{LPF}\) den größten Primfaktor bezeichnet.
Ist die aktuelle Summe \(s\) prim, so gilt für alle \(1\le a\le s-2\): \(\gcd(a+1,s)=1\). Dann passiert keine weitere Kürzung mehr, bevor der Nenner 1 wird. Ausgehend von \(1/(s-1)\) endet die Folge also bei \((s-1)/1\). Damit folgt
$$\boxed{f(n)=\operatorname{LPF}(n+1)-1.}$$
Die Iteration lautet
$$\frac{1}{20}\to\frac{2}{19}\to\frac{3}{18}=\frac{1}{6}\to\frac{2}{5}\to\frac{3}{4}\to\frac{4}{3}\to\frac{5}{2}\to\frac{6}{1}.$$
Also ist \(f(20)=6\). Da \(20+1=21=3\cdot 7\), ergibt die Formel
$$\operatorname{LPF}(21)-1=7-1=6,$$
genau wie die direkte Simulation.
Für die gesuchte Summe setzt man \(n=k^3\):
$$f(k^3)=\operatorname{LPF}(k^3+1)-1.$$
Mit der Zerlegung einer Kubiksumme
$$k^3+1=(k+1)(k^2-k+1)$$
und
$$Q_k=k^2-k+1$$
erhält man
$$\operatorname{LPF}(k^3+1)=\max\bigl(\operatorname{LPF}(k+1),\operatorname{LPF}(Q_k)\bigr).$$
Auch wenn beide Faktoren gelegentlich den Faktor 3 teilen, bleibt der größte Primfaktor des Produkts einfach das Maximum der größten Primfaktoren der beiden Faktoren. Also
$$\sum_{k=1}^{L}f(k^3)=\sum_{k=1}^{L}\left(\max\bigl(\operatorname{LPF}(k+1),\operatorname{LPF}(Q_k)\bigr)-1\right).$$
Ein gewöhnliches Sieb für größte Primfaktoren auf \([2,L+1]\) füllt die Tabelle \(\operatorname{LPF}(m)\) für alle \(m\le L+1\). Damit sind alle Werte \(\operatorname{LPF}(k+1)\) sofort bekannt.
Für jedes \(k\) initialisiert der Code
$$\mathrm{rem}[k]=Q_k,\qquad \mathrm{lpfQ}[k]=1.$$
Fixiere nun eine Primzahl \(p\). Gesucht sind alle \(k\) mit
$$Q_k\equiv0\pmod p.$$
Nach Multiplikation mit 4 wird dies zu
$$4k^2-4k+4\equiv0\pmod p\qquad\Longleftrightarrow\qquad (2k-1)^2\equiv -3\pmod p.$$
Für \(p>3\) existieren Lösungen also genau dann, wenn \(-3\) ein quadratischer Rest modulo \(p\) ist. Gilt \(s^2\equiv -3\pmod p\), dann sind die Restklassen
$$k\equiv \frac{1\pm s}{2}\pmod p.$$
Die Quadratwurzel \(s\) berechnet der Code mit Tonelli-Shanks. Der Spezialfall \(p=3\) liefert genau die Klasse \(k\equiv2\pmod3\). Die Primzahl 2 teilt \(Q_k\) nie, denn \(k^2-k=k(k-1)\) ist immer gerade, also ist \(Q_k\) stets ungerade.
Für jede passende Restklasse besucht der Code alle entsprechenden \(k\), teilt alle Potenzen von \(p\) aus \(\mathrm{rem}[k]\) heraus und speichert \(p\) als aktuell größten Primfaktor. Nachdem alle Primzahlen \(p\le L\) verarbeitet wurden, muss ein Rest \(\mathrm{rem}[k]>1\) selbst prim und größer als \(L\) sein; wäre er zusammengesetzt, hätte er einen Primteiler \(\le \sqrt{Q_k}\lt L\), der bereits entfernt worden wäre. Dieser Rest ist also der letzte große Primfaktor von \(Q_k\).
Hier gilt
$$k^3+1=65=5\cdot 13,$$
also
$$f(64)=\operatorname{LPF}(65)-1=13-1=12.$$
Die Zerlegung liefert \(k+1=5\) und \(Q_4=13\); der größere Primfaktor ist also tatsächlich \(13\).
Die Implementierung besitzt zwei Siebebenen. Zuerst wird lpf_small für alle Zahlen bis \(L+1\)
gefüllt. Danach werden alle \(Q_k\) im Restfeld rem abgelegt und die Primzahlen der Reihe nach
verarbeitet. Für jede Primzahl bestimmt der Code die zulässigen Restklassen von \(k\), entfernt diesen Primfaktor
aus allen passenden Resten und aktualisiert lpf_q[k]. Am Ende werden
lpf_small[k+1] und lpf_q[k] kombiniert, 1 subtrahiert und die Summe in einem
128-Bit-Akkumulator addiert. Die C++-Version prüft dabei \(f(20)=6\),
\(\sum_{k\le100}f(k^3)=118937\) sowie einen Brute-Force-Abgleich für kleine Grenzen.
Das Sieb für \(\operatorname{LPF}(k+1)\) ist in der Praxis nahezu linear. Das \(Q_k\)-Sieb besucht pro relevanter Primzahl etwa \(L/p\) oder \(2L/p\) Positionen; die dominierende Arbeit ist deshalb ebenfalls fast linear bis auf den üblichen Faktor \(\sum_{p\le L}1/p\), also \(O(L\log\log L)\), plus polylogarithmische Kosten für modulare Quadratwurzeln. Der Speicherbedarf beträgt \(O(L)\) für Faktortabellen und Restpuffer.
Başlangıç kesri sadeleştirilmiş \(1/n\) olsun. Her adımda \(x/y\), \((x+1)/(y-1)\) ile değiştirilip yeniden sadeleştirilir. Paydanın ilk kez 1 olduğu anda pay \(f(n)\) ise, problem şu toplamı ister:
$$\sum_{k=1}^{L} f(k^3),\qquad L=2\cdot 10^6.$$
Her küp için süreci doğrudan simüle etmek çok yavaştır; asıl fikir \(f(n)\)'yi \(n+1\)'in asal çarpanları cinsinden yazmaktır.
Her sadeleştirilmiş durumu \(a/b\) biçiminde yazalım (\(a,b>0\)) ve durum toplamını
$$s=a+b$$
olarak tanımlayalım. Başlangıçta \(a=1\), \(b=n\) olduğundan \(s=n+1\).
Sadeleştirilmiş \(a/b\) kesrinden bir sonraki sadeleştirilmemiş kesir \((a+1)/(b-1)\) olur. Eğer
$$d=\gcd(a+1,b-1)$$
ise sadeleştirmeden sonra
$$a'=\frac{a+1}{d},\qquad b'=\frac{b-1}{d}$$
elde edilir. \(b=s-a\) olduğundan
$$d=\gcd(a+1,s-a-1)=\gcd(a+1,s)$$
yazabiliriz. Böylece yeni durum toplamı
$$s'=a'+b'=\frac{s}{\gcd(a+1,s)}$$
olur. Yani her sadeleştirme, mevcut \(s\) toplamını onun bir bölenine indirir. Ayrıca \(a/b\) zaten sade olduğu için \(\gcd(a,b)=1\) ve dolayısıyla \(\gcd(a,s)=1\) geçerlidir.
Bir fazın \(1/(s-1)\) durumunda başladığını varsayalım. \(p\), \(s\)'nin en küçük asal çarpanı olsun. O zaman \(2\le t\lt p\) için \(\gcd(t,s)=1\) olur. Bu yüzden pay \(1,2,\dots,p-1\) değerlerini alırken herhangi bir sadeleştirme gerçekleşmez.
Sadeleştirilmemiş pay \(p\) olduğunda sadeleştirme katsayısı
$$d=\gcd(p,s)=p$$
olur ve o anda
$$\frac{p}{s-p}\to \frac{1}{s/p-1}$$
dönüşümü meydana gelir. Dolayısıyla tek bir faz
$$s\longmapsto \frac{s}{p}$$
işlemini yapar; yani durum toplamının en küçük asal çarpanını kaldırır ve payı tekrar 1'e döndürür.
Başlangıçta \(s_0=n+1\) olduğundan, ardışık fazlar \(n+1\)'in asal çarpanlarını küçükten büyüğe doğru ayıklar. En sonda yalnızca en büyük asal çarpan kaldığında durum toplamı
$$P^+(n+1)=\operatorname{LPF}(n+1)$$
olur; burada \(\operatorname{LPF}\), en büyük asal çarpan fonksiyonudur.
Eğer mevcut toplam \(s\) asal ise, \(1\le a\le s-2\) için \(\gcd(a+1,s)=1\) olur. Böylece payda 1 olana kadar bir daha sadeleştirme olmaz. \(1/(s-1)\) durumundan başlanırsa zincir \((s-1)/1\) ile biter. Sonuç olarak
$$\boxed{f(n)=\operatorname{LPF}(n+1)-1.}$$
İterasyon şu şekildedir:
$$\frac{1}{20}\to\frac{2}{19}\to\frac{3}{18}=\frac{1}{6}\to\frac{2}{5}\to\frac{3}{4}\to\frac{4}{3}\to\frac{5}{2}\to\frac{6}{1}.$$
Buradan \(f(20)=6\) çıkar. Öte yandan \(20+1=21=3\cdot 7\) olduğu için formül
$$\operatorname{LPF}(21)-1=7-1=6$$
sonucunu verir; simülasyonla tam uyumludur.
İstenen toplam için \(n=k^3\) yazılır:
$$f(k^3)=\operatorname{LPF}(k^3+1)-1.$$
Küpler toplamı çarpanlara ayrılırsa
$$k^3+1=(k+1)(k^2-k+1)$$
ve
$$Q_k=k^2-k+1$$
tanımıyla
$$\operatorname{LPF}(k^3+1)=\max\bigl(\operatorname{LPF}(k+1),\operatorname{LPF}(Q_k)\bigr)$$
elde edilir. İki çarpan bazen ortak 3 çarpanını paylaşsa bile, çarpımın en büyük asal çarpanı yine bu iki terimin en büyük asal çarpanlarının maksimumudur. Bu yüzden
$$\sum_{k=1}^{L}f(k^3)=\sum_{k=1}^{L}\left(\max\bigl(\operatorname{LPF}(k+1),\operatorname{LPF}(Q_k)\bigr)-1\right).$$
\([2,L+1]\) aralığında standart bir en büyük asal çarpan süzgeci çalıştırılarak \(m\le L+1\) için \(\operatorname{LPF}(m)\) tablosu doldurulur. Böylece tüm \(\operatorname{LPF}(k+1)\) değerleri doğrudan elde edilir.
Her \(k\) için başlangıçta
$$\mathrm{rem}[k]=Q_k,\qquad \mathrm{lpfQ}[k]=1$$
atanır. Şimdi bir asal \(p\) sabitleyelim. Amaç
$$Q_k\equiv0\pmod p$$
denkliğini sağlayan \(k\) sınıflarını bulmaktır. 4 ile çarpınca
$$4k^2-4k+4\equiv0\pmod p\qquad\Longleftrightarrow\qquad (2k-1)^2\equiv -3\pmod p$$
olur. Dolayısıyla \(p>3\) için kökler ancak ve ancak \(-3\), \(p\) modunda kuadratik artık ise vardır. Eğer \(s^2\equiv -3\pmod p\) ise
$$k\equiv \frac{1\pm s}{2}\pmod p$$
sınıfları elde edilir. Kod, bu karekökü Tonelli-Shanks algoritması ile bulur. \(p=3\) özel durumunda tek sınıf \(k\equiv2\pmod3\)'tür. \(p=2\) ise hiçbir zaman \(Q_k\)'yi bölmez; çünkü \(k^2-k=k(k-1)\) her zaman çifttir, dolayısıyla \(Q_k\) daima tektir.
Her uygun artık sınıfı için kod ilgili tüm \(k\) değerlerini gezer, \(\mathrm{rem}[k]\) içinden \(p\)'nin bütün kuvvetlerini çıkarır ve \(p\)'yi o ana kadarki en büyük asal çarpan olarak kaydeder. Tüm \(p\le L\) asalları işlendiğinde \(\mathrm{rem}[k]>1\) kalan bir değer, zorunlu olarak asal ve \(L\)'den büyüktür; aksi halde bileşik olsaydı \(\sqrt{Q_k}\lt L\) nedeniyle \(L\)'den küçük bir asal böleni olur ve o daha önce temizlenmiş olurdu. Dolayısıyla bu kalan değer \(Q_k\)'nin son büyük asal çarpanıdır.
Bu durumda
$$k^3+1=65=5\cdot 13$$
olduğundan
$$f(64)=\operatorname{LPF}(65)-1=13-1=12$$
elde edilir. Ayrışımda \(k+1=5\) ve \(Q_4=13\) olduğu için büyük asal çarpan gerçekten \(13\)'tür.
Uygulama iki süzgeç katmanından oluşur. Önce \(L+1\)'e kadar tüm sayılar için lpf_small tablosu
doldurulur. Sonra bütün \(Q_k\) değerleri rem dizisine yerleştirilir ve asallar tek tek işlenir. Her
asal için uygun \(k\) artık sınıfları bulunur, bu asal ilgili kalıntılardan tamamen çıkarılır ve
lpf_q[k] güncellenir. En sonda lpf_small[k+1] ile lpf_q[k] birleştirilir,
1 çıkarılır ve toplam 128 bitlik bir akümülatörde biriktirilir. C++ sürümündeki kontroller \(f(20)=6\),
\(\sum_{k\le100}f(k^3)=118937\) ve küçük limitlerde brute-force eşleşmesini doğrular.
\(\operatorname{LPF}(k+1)\) süzgeci pratikte yaklaşık doğrusal çalışır. \(Q_k\) süzgeci, ilgili her asal için yaklaşık \(L/p\) ya da \(2L/p\) konumu ziyaret eder; bu nedenle baskın maliyet olağan \(\sum_{p\le L}1/p\) faktörü kadar büyür, yani modüler karekök maliyetleri dışında \(O(L\log\log L)\) davranışı gösterir. Bellek kullanımı faktör tabloları ve kalıntı tamponu için \(O(L)\) düzeyindedir.
Se parte de la fracción reducida \(1/n\). En cada paso se reemplaza \(x/y\) por \((x+1)/(y-1)\) y se vuelve a simplificar. Si la primera vez que el denominador llega a 1 el numerador vale \(f(n)\), el problema pide
$$\sum_{k=1}^{L} f(k^3),\qquad L=2\cdot 10^6.$$
Simular el proceso para cada cubo sería demasiado lento; la clave es reescribir \(f(n)\) en función de los factores primos de \(n+1\).
Escribamos cualquier estado reducido como \(a/b\) con \(a,b>0\), y definamos la suma de estado
$$s=a+b.$$
Al inicio \(a=1\), \(b=n\), por lo tanto \(s=n+1\).
Desde la fracción reducida \(a/b\), la siguiente fracción sin simplificar es \((a+1)/(b-1)\). Si
$$d=\gcd(a+1,b-1),$$
entonces tras simplificar obtenemos
$$a'=\frac{a+1}{d},\qquad b'=\frac{b-1}{d}.$$
Como \(b=s-a\), resulta
$$d=\gcd(a+1,s-a-1)=\gcd(a+1,s).$$
Así, la nueva suma de estado vale
$$s'=a'+b'=\frac{s}{\gcd(a+1,s)}.$$
Cada simplificación reemplaza la suma actual \(s\) por un divisor suyo. Además, como \(a/b\) ya está reducida, \(\gcd(a,b)=1\), y por tanto también \(\gcd(a,s)=1\).
Supongamos que una fase comienza en el estado \(1/(s-1)\). Sea \(p\) el menor factor primo de \(s\). Para cualquier entero \(t\) con \(2\le t\lt p\), se cumple \(\gcd(t,s)=1\). Por ello no aparece ninguna simplificación mientras el numerador recorre \(1,2,\dots,p-1\).
Cuando el numerador sin reducir llega a \(p\), el factor de reducción es
$$d=\gcd(p,s)=p.$$
En ese instante
$$\frac{p}{s-p}\to \frac{1}{s/p-1}.$$
Por tanto, una fase completa realiza la transformación
$$s\longmapsto \frac{s}{p},$$
es decir, arranca el menor factor primo de la suma actual y reinicia el numerador en 1.
Comenzando con \(s_0=n+1\), las fases sucesivas eliminan los factores primos de \(n+1\) de menor a mayor. Cuando ya solo queda el mayor, la suma de estado es
$$P^+(n+1)=\operatorname{LPF}(n+1),$$
donde \(\operatorname{LPF}\) significa "mayor factor primo".
Si la suma actual \(s\) es prima, entonces para todo \(1\le a\le s-2\) se tiene \(\gcd(a+1,s)=1\). Así, no vuelve a haber reducciones antes de que el denominador alcance 1. Partiendo de \(1/(s-1)\), la cadena termina en \((s-1)/1\). Luego
$$\boxed{f(n)=\operatorname{LPF}(n+1)-1.}$$
La iteración es
$$\frac{1}{20}\to\frac{2}{19}\to\frac{3}{18}=\frac{1}{6}\to\frac{2}{5}\to\frac{3}{4}\to\frac{4}{3}\to\frac{5}{2}\to\frac{6}{1}.$$
Por tanto \(f(20)=6\). Como \(20+1=21=3\cdot 7\), la fórmula da
$$\operatorname{LPF}(21)-1=7-1=6,$$
exactamente el mismo resultado.
Para la suma pedida tomamos \(n=k^3\):
$$f(k^3)=\operatorname{LPF}(k^3+1)-1.$$
Usando la factorización de suma de cubos,
$$k^3+1=(k+1)(k^2-k+1),$$
y definiendo
$$Q_k=k^2-k+1,$$
obtenemos
$$\operatorname{LPF}(k^3+1)=\max\bigl(\operatorname{LPF}(k+1),\operatorname{LPF}(Q_k)\bigr).$$
Aunque a veces ambos factores comparten un 3, el mayor factor primo del producto sigue siendo el máximo de los mayores factores primos de cada factor. Entonces
$$\sum_{k=1}^{L}f(k^3)=\sum_{k=1}^{L}\left(\max\bigl(\operatorname{LPF}(k+1),\operatorname{LPF}(Q_k)\bigr)-1\right).$$
Una criba estándar de mayor factor primo sobre \([2,L+1]\) llena una tabla \(\operatorname{LPF}(m)\) para todos los \(m\le L+1\). Con ello ya se conocen todos los valores \(\operatorname{LPF}(k+1)\).
Para cada \(k\), el código inicializa
$$\mathrm{rem}[k]=Q_k,\qquad \mathrm{lpfQ}[k]=1.$$
Fijemos ahora un primo \(p\). Queremos resolver
$$Q_k\equiv0\pmod p.$$
Multiplicando por 4:
$$4k^2-4k+4\equiv0\pmod p\qquad\Longleftrightarrow\qquad (2k-1)^2\equiv -3\pmod p.$$
Por tanto, para \(p>3\), existen raíces exactamente cuando \(-3\) es un residuo cuadrático módulo \(p\). Si \(s^2\equiv -3\pmod p\), entonces
$$k\equiv \frac{1\pm s}{2}\pmod p.$$
El código obtiene la raíz cuadrada \(s\) mediante Tonelli-Shanks. El caso especial \(p=3\) produce la única clase \(k\equiv2\pmod3\). El primo 2 nunca divide a \(Q_k\), porque \(k^2-k=k(k-1)\) siempre es par, así que \(Q_k\) siempre es impar.
Para cada clase válida, el código recorre todos los \(k\) correspondientes, divide todas las potencias de \(p\) que aparezcan en \(\mathrm{rem}[k]\) y guarda \(p\) como mayor factor primo encontrado hasta ese momento. Después de procesar todos los primos \(p\le L\), cualquier resto \(\mathrm{rem}[k]>1\) tiene que ser primo y mayor que \(L\); si fuera compuesto, tendría un divisor primo \(\le \sqrt{Q_k}\lt L\), que ya habría sido eliminado. Ese resto es, por tanto, el último gran factor primo de \(Q_k\).
En este caso
$$k^3+1=65=5\cdot 13,$$
así que
$$f(64)=\operatorname{LPF}(65)-1=13-1=12.$$
La descomposición da \(k+1=5\) y \(Q_4=13\); claramente el mayor factor primo es \(13\).
La implementación usa dos capas de cribado. Primero llena lpf_small para todos los enteros hasta
\(L+1\). Después guarda todos los \(Q_k\) en el arreglo residual rem y procesa los primos uno por uno.
Para cada primo calcula las clases de residuos admisibles de \(k\), elimina dicho primo de todos los residuos
compatibles y actualiza lpf_q[k]. Al final combina lpf_small[k+1] con
lpf_q[k], resta 1 y acumula la suma en un total de 128 bits. La versión en C++ comprueba
\(f(20)=6\), \(\sum_{k\le100}f(k^3)=118937\) y también verifica casos pequeños por fuerza bruta.
La criba de \(\operatorname{LPF}(k+1)\) es casi lineal en la práctica. La criba sobre \(Q_k\) visita unas \(L/p\) o \(2L/p\) posiciones por primo relevante; por ello su coste dominante también es casi lineal salvo por el factor habitual \(\sum_{p\le L}1/p\), es decir \(O(L\log\log L)\), más el coste polilogarítmico de las raíces cuadradas modulares. La memoria es \(O(L)\) para tablas de factores y residuos.
从既约分数 \(1/n\) 开始。每一步把 \(x/y\) 变成 \((x+1)/(y-1)\)，然后再次约分。若分母第一次变成 1 时的分子记为 \(f(n)\)，题目要求计算
$$\sum_{k=1}^{L} f(k^3),\qquad L=2\cdot 10^6.$$
如果对每个立方数都直接模拟过程，复杂度会过高；关键是把 \(f(n)\) 改写成 \(n+1\) 的素因子问题。
把任意既约状态写成 \(a/b\)（其中 \(a,b>0\)），并定义它的状态和
$$s=a+b.$$
初始时 \(a=1\)、\(b=n\)，所以 \(s=n+1\)。
从既约分数 \(a/b\) 出发，下一步未约分的分数是 \((a+1)/(b-1)\)。设
$$d=\gcd(a+1,b-1),$$
则约分后得到
$$a'=\frac{a+1}{d},\qquad b'=\frac{b-1}{d}.$$
由于 \(b=s-a\)，于是
$$d=\gcd(a+1,s-a-1)=\gcd(a+1,s).$$
因此新的状态和为
$$s'=a'+b'=\frac{s}{\gcd(a+1,s)}.$$
也就是说，每次约分都会把当前的 \(s\) 替换成它的某个因子。又因为 \(a/b\) 已经最简，所以 \(\gcd(a,b)=1\)，从而 \(\gcd(a,s)=1\)。
设某一阶段从 \(1/(s-1)\) 开始，记 \(p\) 为 \(s\) 的最小素因子。对任意满足 \(2\le t\lt p\) 的整数 \(t\)，都有 \(\gcd(t,s)=1\)。因此当分子依次经过 \(1,2,\dots,p-1\) 时，不会发生新的约分。
当未约分的分子到达 \(p\) 时，约分因子为
$$d=\gcd(p,s)=p.$$
于是此时
$$\frac{p}{s-p}\to \frac{1}{s/p-1}.$$
所以一个完整阶段实现的正是
$$s\longmapsto \frac{s}{p},$$
即把当前状态和的最小素因子剥掉，并把分子重新置为 1。
从 \(s_0=n+1\) 出发，连续的阶段会按从小到大的顺序依次去掉 \(n+1\) 的素因子。当只剩下最大的那个时，状态和变为
$$P^+(n+1)=\operatorname{LPF}(n+1),$$
其中 \(\operatorname{LPF}\) 表示“最大素因子”。
若当前 \(s\) 本身是素数，则对所有 \(1\le a\le s-2\) 都有 \(\gcd(a+1,s)=1\)，因此在分母变成 1 之前不再约分。 从 \(1/(s-1)\) 出发，链最终到达 \((s-1)/1\)。于是得到
$$\boxed{f(n)=\operatorname{LPF}(n+1)-1.}$$
迭代过程为
$$\frac{1}{20}\to\frac{2}{19}\to\frac{3}{18}=\frac{1}{6}\to\frac{2}{5}\to\frac{3}{4}\to\frac{4}{3}\to\frac{5}{2}\to\frac{6}{1}.$$
所以 \(f(20)=6\)。另一方面 \(20+1=21=3\cdot 7\)，因此公式给出
$$\operatorname{LPF}(21)-1=7-1=6,$$
与模拟完全一致。
题目要求的是 \(n=k^3\) 时的值，因此
$$f(k^3)=\operatorname{LPF}(k^3+1)-1.$$
利用立方和分解
$$k^3+1=(k+1)(k^2-k+1),$$
并定义
$$Q_k=k^2-k+1,$$
就有
$$\operatorname{LPF}(k^3+1)=\max\bigl(\operatorname{LPF}(k+1),\operatorname{LPF}(Q_k)\bigr).$$
即使两个因子有时会共同含有因子 3，乘积的最大素因子仍然只是这两个最大素因子的较大者。因此
$$\sum_{k=1}^{L}f(k^3)=\sum_{k=1}^{L}\left(\max\bigl(\operatorname{LPF}(k+1),\operatorname{LPF}(Q_k)\bigr)-1\right).$$
对区间 \([2,L+1]\) 做一次标准的最大素因子筛，就能得到所有 \(m\le L+1\) 的 \(\operatorname{LPF}(m)\)，从而立刻得到全部 \(\operatorname{LPF}(k+1)\)。
对每个 \(k\)，先初始化
$$\mathrm{rem}[k]=Q_k,\qquad \mathrm{lpfQ}[k]=1.$$
现在固定一个素数 \(p\)。我们要找的是满足
$$Q_k\equiv0\pmod p$$
的所有 \(k\)。两边乘以 4，可化为
$$4k^2-4k+4\equiv0\pmod p\qquad\Longleftrightarrow\qquad (2k-1)^2\equiv -3\pmod p.$$
因此对 \(p>3\)，解存在当且仅当 \(-3\) 是模 \(p\) 的二次剩余。若 \(s^2\equiv -3\pmod p\)，则
$$k\equiv \frac{1\pm s}{2}\pmod p.$$
代码用 Tonelli-Shanks 算法求这个平方根 \(s\)。特殊素数 \(p=3\) 只产生一个同余类 \(k\equiv2\pmod3\)。素数 2 永远不会整除 \(Q_k\)，因为 \(k^2-k=k(k-1)\) 总是偶数，所以 \(Q_k\) 总是奇数。
对每个有效同余类，代码遍历所有对应的 \(k\)，把 \(\mathrm{rem}[k]\) 中的 \(p\) 的全部幂次除尽，并把 \(p\) 记作当前找到的最大素因子。处理完所有 \(p\le L\) 后，如果还有 \(\mathrm{rem}[k]>1\)，那么它必然本身就是一个大于 \(L\) 的素数；否则若它是合数，就会有一个 \(\le \sqrt{Q_k}\lt L\) 的素因子，而这类素因子早就会被筛掉。这个剩余值就是 \(Q_k\) 的最后一个大素因子。
此时
$$k^3+1=65=5\cdot 13,$$
所以
$$f(64)=\operatorname{LPF}(65)-1=13-1=12.$$
分解中 \(k+1=5\)、\(Q_4=13\)，最大素因子显然是 \(13\)。
实现中有两层筛。第一层先填好 \(L+1\) 以内所有整数的 lpf_small。第二层把所有 \(Q_k\) 存进残值数组
rem，然后按素数逐个处理。对每个素数，代码先求出允许的 \(k\) 同余类，再把这个素数从所有匹配位置的残值中完全除掉，同时更新
lpf_q[k]。最后把 lpf_small[k+1] 与 lpf_q[k] 取最大值，减 1 后累加到
128 位总和中。C++ 版本还验证了 \(f(20)=6\)、\(\sum_{k\le100}f(k^3)=118937\)，以及小范围下与暴力算法一致。
\(\operatorname{LPF}(k+1)\) 的筛法在实践中接近线性。\(Q_k\) 的筛法对每个相关素数大约访问 \(L/p\) 或 \(2L/p\) 个位置，因此主导成本同样接近线性，只多出通常的 \(\sum_{p\le L}1/p\) 因子，也就是 \(O(L\log\log L)\)，外加模平方根带来的多对数开销。内存复杂度为 \(O(L)\)。
Стартуем с несократимой дроби \(1/n\). На каждом шаге дробь \(x/y\) заменяется на \((x+1)/(y-1)\), после чего снова сокращается. Если в первый момент, когда знаменатель становится равным 1, числитель равен \(f(n)\), то требуется вычислить
$$\sum_{k=1}^{L} f(k^3),\qquad L=2\cdot 10^6.$$
Прямая симуляция для каждого куба слишком медленна, поэтому нужно выразить \(f(n)\) через простые делители числа \(n+1\).
Запишем любое сокращённое состояние как \(a/b\), где \(a,b>0\), и введём сумму состояния
$$s=a+b.$$
Изначально \(a=1\), \(b=n\), следовательно, \(s=n+1\).
Из сокращённой дроби \(a/b\) получаем несокращённую дробь \((a+1)/(b-1)\). Пусть
$$d=\gcd(a+1,b-1).$$
Тогда после сокращения
$$a'=\frac{a+1}{d},\qquad b'=\frac{b-1}{d}.$$
Поскольку \(b=s-a\), имеем
$$d=\gcd(a+1,s-a-1)=\gcd(a+1,s).$$
Значит новая сумма состояния равна
$$s'=a'+b'=\frac{s}{\gcd(a+1,s)}.$$
Итак, каждое сокращение заменяет текущую сумму \(s\) на некоторый её делитель. Так как дробь \(a/b\) уже сокращена, \(\gcd(a,b)=1\), а значит и \(\gcd(a,s)=1\).
Предположим, что фаза начинается в состоянии \(1/(s-1)\). Пусть \(p\) — наименьший простой делитель числа \(s\). Тогда для любого \(t\) с \(2\le t\lt p\) выполняется \(\gcd(t,s)=1\). Следовательно, пока числитель проходит значения \(1,2,\dots,p-1\), новых сокращений не происходит.
Когда несокращённый числитель становится равным \(p\), коэффициент сокращения равен
$$d=\gcd(p,s)=p.$$
В этот момент
$$\frac{p}{s-p}\to \frac{1}{s/p-1}.$$
Таким образом, одна полная фаза реализует преобразование
$$s\longmapsto \frac{s}{p},$$
то есть снимает наименьший простой множитель текущей суммы и возвращает числитель к 1.
Начиная с \(s_0=n+1\), последовательные фазы удаляют простые множители числа \(n+1\) в порядке возрастания. Когда остаётся только наибольший простой множитель, сумма состояния становится равной
$$P^+(n+1)=\operatorname{LPF}(n+1),$$
где \(\operatorname{LPF}\) означает наибольший простой делитель.
Если текущее \(s\) — простое, то для любого \(1\le a\le s-2\) имеем \(\gcd(a+1,s)=1\). Значит до достижения знаменателя 1 больше сокращений не будет. Из состояния \(1/(s-1)\) цепочка заканчивается в \((s-1)/1\). Поэтому
$$\boxed{f(n)=\operatorname{LPF}(n+1)-1.}$$
Итерация имеет вид
$$\frac{1}{20}\to\frac{2}{19}\to\frac{3}{18}=\frac{1}{6}\to\frac{2}{5}\to\frac{3}{4}\to\frac{4}{3}\to\frac{5}{2}\to\frac{6}{1}.$$
Отсюда \(f(20)=6\). Поскольку \(20+1=21=3\cdot 7\), формула даёт
$$\operatorname{LPF}(21)-1=7-1=6,$$
что полностью совпадает с симуляцией.
Для искомой суммы нужно подставить \(n=k^3\):
$$f(k^3)=\operatorname{LPF}(k^3+1)-1.$$
Используя разложение суммы кубов,
$$k^3+1=(k+1)(k^2-k+1),$$
и обозначая
$$Q_k=k^2-k+1,$$
получаем
$$\operatorname{LPF}(k^3+1)=\max\bigl(\operatorname{LPF}(k+1),\operatorname{LPF}(Q_k)\bigr).$$
Даже если оба множителя иногда имеют общий множитель 3, наибольший простой делитель произведения всё равно равен максимуму из наибольших простых делителей множителей. Следовательно,
$$\sum_{k=1}^{L}f(k^3)=\sum_{k=1}^{L}\left(\max\bigl(\operatorname{LPF}(k+1),\operatorname{LPF}(Q_k)\bigr)-1\right).$$
Обычное решето наибольшего простого делителя на отрезке \([2,L+1]\) заполняет таблицу \(\operatorname{LPF}(m)\) для всех \(m\le L+1\). Это сразу даёт все значения \(\operatorname{LPF}(k+1)\).
Для каждого \(k\) код инициализирует
$$\mathrm{rem}[k]=Q_k,\qquad \mathrm{lpfQ}[k]=1.$$
Теперь фиксируем простое \(p\). Нужно найти все \(k\), для которых
$$Q_k\equiv0\pmod p.$$
После умножения на 4 это равносильно
$$4k^2-4k+4\equiv0\pmod p\qquad\Longleftrightarrow\qquad (2k-1)^2\equiv -3\pmod p.$$
Значит, для \(p>3\) решения существуют тогда и только тогда, когда \(-3\) является квадратичным вычетом по модулю \(p\). Если \(s^2\equiv -3\pmod p\), то
$$k\equiv \frac{1\pm s}{2}\pmod p.$$
Квадратный корень \(s\) код находит алгоритмом Тонелли-Шенкса. Особый случай \(p=3\) даёт единственный класс \(k\equiv2\pmod3\). Простое число 2 никогда не делит \(Q_k\), поскольку \(k^2-k=k(k-1)\) всегда чётно, а значит \(Q_k\) всегда нечётно.
Для каждого подходящего класса остатков код проходит по всем соответствующим \(k\), делит из \(\mathrm{rem}[k]\) все степени \(p\) и записывает \(p\) как текущий наибольший простой делитель. После обработки всех простых \(p\le L\) любой остаток \(\mathrm{rem}[k]>1\) обязан быть простым и большим \(L\); если бы он был составным, то имел бы простой делитель \(\le \sqrt{Q_k}\lt L\), который уже был бы удалён. Этот остаток и есть последний большой простой делитель числа \(Q_k\).
Здесь
$$k^3+1=65=5\cdot 13,$$
поэтому
$$f(64)=\operatorname{LPF}(65)-1=13-1=12.$$
Разложение даёт \(k+1=5\) и \(Q_4=13\), так что максимальный простой делитель действительно равен \(13\).
Реализация использует два слоя решета. Сначала заполняется массив lpf_small для всех чисел до
\(L+1\). Затем все значения \(Q_k\) помещаются в массив остатков rem, и простые числа обрабатываются
по одному. Для каждого простого кода вычисляет допустимые классы \(k\), полностью вычищает этот простой из всех
подходящих остатков и обновляет lpf_q[k]. В конце берётся максимум из
lpf_small[k+1] и lpf_q[k], вычитается 1, и результат суммируется в 128-битном
аккумуляторе. Версия на C++ проверяет \(f(20)=6\), \(\sum_{k\le100}f(k^3)=118937\), а также совпадение с brute
force на малых пределах.
Решето для \(\operatorname{LPF}(k+1)\) на практике близко к линейному. Решето по \(Q_k\) посещает примерно \(L/p\) или \(2L/p\) позиций на каждый подходящий простой, поэтому основной объём работы также близок к линейному, с обычным множителем \(\sum_{p\le L}1/p\), то есть \(O(L\log\log L)\), плюс полилогарифмическая стоимость вычисления квадратных корней по модулю. Память: \(O(L)\).
نبدأ بالكسر المختزل \(1/n\). في كل خطوة نستبدل \(x/y\) بـ \((x+1)/(y-1)\) ثم نختزل من جديد. إذا كان البسط عند أول لحظة يصبح فيها المقام مساويًا لـ 1 هو \(f(n)\)، فالمطلوب حساب
$$\sum_{k=1}^{L} f(k^3),\qquad L=2\cdot 10^6.$$
المحاكاة المباشرة لكل مكعب بطيئة جدًا، ولذلك الفكرة الأساسية هي إعادة كتابة \(f(n)\) بدلالة العوامل الأولية للعدد \(n+1\).
لنكتب أي حالة مختزلة على الصورة \(a/b\) حيث \(a,b>0\)، ولنعرّف مجموع الحالة بأنه
$$s=a+b.$$
في البداية لدينا \(a=1\) و\(b=n\)، ومن ثم \(s=n+1\).
انطلاقًا من الكسر المختزل \(a/b\)، يكون الكسر التالي قبل الاختزال هو \((a+1)/(b-1)\). إذا وضعنا
$$d=\gcd(a+1,b-1),$$
فعند الاختزال نحصل على
$$a'=\frac{a+1}{d},\qquad b'=\frac{b-1}{d}.$$
وبما أن \(b=s-a\)، فإن
$$d=\gcd(a+1,s-a-1)=\gcd(a+1,s).$$
إذًا يصبح مجموع الحالة الجديد
$$s'=a'+b'=\frac{s}{\gcd(a+1,s)}.$$
أي أن كل عملية اختزال تستبدل المجموع الحالي \(s\) بقاسم من قواسمه. وبما أن \(a/b\) مختزل أصلًا، فإن \(\gcd(a,b)=1\)، وبالتالي \(\gcd(a,s)=1\).
افترض أن مرحلة ما تبدأ من الحالة \(1/(s-1)\). ولْيكن \(p\) أصغر عامل أولي للعدد \(s\). عندئذ لكل عدد صحيح \(t\) يحقق \(2\le t\lt p\) يكون \(\gcd(t,s)=1\). لذلك لا يحدث أي اختزال جديد بينما يمر البسط على القيم \(1,2,\dots,p-1\).
عندما يصل البسط غير المختزل إلى القيمة \(p\)، يكون عامل الاختزال
$$d=\gcd(p,s)=p.$$
وفي هذه اللحظة يتحول
$$\frac{p}{s-p}\to \frac{1}{s/p-1}.$$
إذن مرحلة كاملة تنفذ التحويل
$$s\longmapsto \frac{s}{p},$$
أي إنها تنزع أصغر عامل أولي من مجموع الحالة الحالي وتعيد البسط إلى 1.
بدءًا من \(s_0=n+1\)، تزيل المراحل المتعاقبة العوامل الأولية للعدد \(n+1\) من الأصغر إلى الأكبر. وعندما لا يبقى إلا أكبر عامل أولي، يصبح مجموع الحالة
$$P^+(n+1)=\operatorname{LPF}(n+1),$$
حيث ترمز \(\operatorname{LPF}\) إلى أكبر عامل أولي.
إذا كان \(s\) الحالي عددًا أوليًا، فإن \(\gcd(a+1,s)=1\) لكل \(1\le a\le s-2\)، ولذلك لا يحدث أي اختزال آخر قبل أن يصل المقام إلى 1. وإذا بدأنا من \(1/(s-1)\)، تنتهي السلسلة عند \((s-1)/1\). وعليه
$$\boxed{f(n)=\operatorname{LPF}(n+1)-1.}$$
تكون السلسلة
$$\frac{1}{20}\to\frac{2}{19}\to\frac{3}{18}=\frac{1}{6}\to\frac{2}{5}\to\frac{3}{4}\to\frac{4}{3}\to\frac{5}{2}\to\frac{6}{1}.$$
إذن \(f(20)=6\). ومن جهة أخرى \(20+1=21=3\cdot 7\)، لذا تعطي الصيغة
$$\operatorname{LPF}(21)-1=7-1=6,$$
وهو نفس الناتج تمامًا.
للمجموع المطلوب نضع \(n=k^3\)، فنحصل على
$$f(k^3)=\operatorname{LPF}(k^3+1)-1.$$
وباستخدام تحليل مجموع المكعبين
$$k^3+1=(k+1)(k^2-k+1),$$
ومع التعريف
$$Q_k=k^2-k+1,$$
يصبح لدينا
$$\operatorname{LPF}(k^3+1)=\max\bigl(\operatorname{LPF}(k+1),\operatorname{LPF}(Q_k)\bigr).$$
حتى لو اشترك العاملان أحيانًا في العامل 3، فإن أكبر عامل أولي للجداء يظل ببساطة هو أكبر القيمتين السابقتين. ومن ثم
$$\sum_{k=1}^{L}f(k^3)=\sum_{k=1}^{L}\left(\max\bigl(\operatorname{LPF}(k+1),\operatorname{LPF}(Q_k)\bigr)-1\right).$$
تُشغّل غربلة قياسية لأكبر عامل أولي على المجال \([2,L+1]\)، فتملأ جدول \(\operatorname{LPF}(m)\) لكل \(m\le L+1\). وهكذا نحصل مباشرة على جميع قيم \(\operatorname{LPF}(k+1)\).
لكل \(k\) يبدأ البرنامج بالقيم
$$\mathrm{rem}[k]=Q_k,\qquad \mathrm{lpfQ}[k]=1.$$
والآن نثبت عددًا أوليًا \(p\). نريد تحديد جميع \(k\) التي تحقق
$$Q_k\equiv0\pmod p.$$
وبعد الضرب في 4 تصبح المعادلة
$$4k^2-4k+4\equiv0\pmod p\qquad\Longleftrightarrow\qquad (2k-1)^2\equiv -3\pmod p.$$
لذلك عندما \(p>3\)، توجد جذور إذا وفقط إذا كانت \(-3\) بقايا تربيعية modulo \(p\). وإذا كان \(s^2\equiv -3\pmod p\)، فإن أصناف البواقي هي
$$k\equiv \frac{1\pm s}{2}\pmod p.$$
ويحسب البرنامج الجذر \(s\) بخوارزمية Tonelli-Shanks. أما الحالة الخاصة \(p=3\) فتعطي الصنف الوحيد \(k\equiv2\pmod3\). والعدد الأولي 2 لا يقسم \(Q_k\) أبدًا، لأن \(k^2-k=k(k-1)\) عدد زوجي دائمًا، وبالتالي \(Q_k\) عدد فردي دائمًا.
لكل صنف متوافق، يمر البرنامج على كل قيم \(k\) الموافقة له، ويقسم جميع قوى \(p\) من \(\mathrm{rem}[k]\)، ثم يسجل \(p\) كأكبر عامل أولي معروف حتى تلك اللحظة. وبعد معالجة كل الأعداد الأولية \(p\le L\)، فإن أي قيمة متبقية \(\mathrm{rem}[k]>1\) لا بد أن تكون عددًا أوليًا أكبر من \(L\)؛ لأنها لو كانت مركبة لامتلكت عاملًا أوليًا \(\le \sqrt{Q_k}\lt L\)، وكان هذا العامل سيُزال سابقًا. وهذه البقية هي إذن العامل الأولي الكبير الأخير لـ \(Q_k\).
لدينا هنا
$$k^3+1=65=5\cdot 13,$$
ومن ثم
$$f(64)=\operatorname{LPF}(65)-1=13-1=12.$$
وفي هذا التفكيك نجد \(k+1=5\) و\(Q_4=13\)، لذا فإن أكبر عامل أولي هو فعلًا \(13\).
التنفيذ يعتمد على طبقتين من الغربلة. أولًا يملأ المصفوفة lpf_small لجميع الأعداد حتى \(L+1\). ثم
يخزن كل قيم \(Q_k\) في مصفوفة البواقي rem ويعالج الأعداد الأولية واحدًا بعد آخر. لكل عدد أولي يحسب
أصناف البواقي المسموح بها لـ \(k\)، ويزيل هذا العامل من جميع البواقي المطابقة، ويحدّث
lpf_q[k]. في النهاية يأخذ القيمة الكبرى بين lpf_small[k+1] وlpf_q[k]،
ثم يطرح 1 ويجمع النتيجة في مُراكم 128-بت. وتتحقق نسخة ++C من القيم \(f(20)=6\) و
\(\sum_{k\le100}f(k^3)=118937\)، إضافة إلى مطابقة brute force عند الحدود الصغيرة.
غربلة \(\operatorname{LPF}(k+1)\) شبه خطية عمليًا. أما غربلة \(Q_k\) فتزور تقريبًا \(L/p\) أو \(2L/p\) موضعًا لكل عدد أولي ذي صلة، ولذلك يبقى العمل المسيطر شبه خطي أيضًا مع عامل \(\sum_{p\le L}1/p\) المعتاد، أي \(O(L\log\log L)\)، إضافة إلى كلفة كثيرة الحدود اللوغاريتمية للجذور التربيعية المعيارية. الذاكرة المطلوبة هي \(O(L)\).