The implementations compute a modular quantity \(H(N)\) for \(N=10^{15}\), with modulus
$$M=1031^3+2.$$
Instead of scanning all integer pairs \((x,y)\) with \(\max(|x|,|y|)\le N\), the solution starts from the algebraic classification used for the nested-radical identities behind the problem. Once that classification is available, the task becomes:
$$H(N)=\sum_{\substack{(x,y)\ \text{admissible}\\ \max(|x|,|y|)\le N}} (|x|+|y|)\pmod{M}.$$
The key point is that every admissible pair comes from a primitive parameter pair together with a square scaling factor, so the search can be reorganized arithmetically rather than geometrically.
The implementations use two primitive families, depending on the parity of one parameter. To avoid exposing code symbols, let the coprime parameters be \(a,b\in \mathbb{Z}_{>0}\) with \(a\ne 2b\), and let the common scale be \(t\ge 1\).
When \(b\) is odd, the primitive pair is
$$X_{\mathrm{odd}}=b(b+4a)^3,\qquad Y_{\mathrm{odd}}=4a(a-2b)^3.$$
When \(b\) is even, the primitive pair is
$$X_{\mathrm{even}}=2b\left(\frac{b}{2}+2a\right)^3,\qquad Y_{\mathrm{even}}=a(a-2b)^3.$$
These two branches are exactly the two branches used by the C++, Python, and Java implementations. The coprimality condition
$$\gcd(a,b)=1$$
eliminates redundant representations, and the excluded line
$$a=2b$$
is the degenerate case where one cubic factor vanishes.
Each primitive pair produces an infinite family by multiplying both coordinates by a common square:
$$ (x,y)=\left(t^2X_\ast,\ t^2Y_\ast\right), $$
where \((X_\ast,Y_\ast)\) means the appropriate primitive pair from the odd or even branch.
Therefore, for a fixed primitive seed, the admissible values under the box constraint are determined by
$$t^2\max(|X_\ast|,|Y_\ast|)\le N,$$
so the largest allowed scale is
$$t_{\max}=\left\lfloor \sqrt{\frac{N}{\max(|X_\ast|,|Y_\ast|)}} \right\rfloor.$$
This is why the outer arithmetic work is done only on primitive seeds; all non-primitive solutions are absorbed into a single closed-form sum over \(t\).
The implementations precompute the cube-free kernel
$$\operatorname{cf}(n)=\prod_{p} p^{v_p(n)\bmod 3},$$
meaning that every full cube \(p^3\) is stripped away, while exponents \(1\) and \(2\) remain. The primitive pair is kept only if the two relevant kernels differ.
For odd \(b\), the required inequality is
$$\operatorname{cf}(b)\ne \operatorname{cf}(4a).$$
For even \(b\), the required inequality is
$$\operatorname{cf}(2b)\ne \operatorname{cf}(a).$$
This filter removes the cases where the two cube-root components collapse to the same cube-free part and the intended nested-radical structure is lost.
For one primitive seed, every scaled pair contributes
$$|t^2X_\ast|+|t^2Y_\ast|=t^2\bigl(|X_\ast|+|Y_\ast|\bigr).$$
Hence the total contribution of that seed is
$$\bigl(|X_\ast|+|Y_\ast|\bigr)\sum_{t=1}^{t_{\max}} t^2.$$
The implementations use the standard closed form
$$\sum_{t=1}^{u} t^2=\frac{u(u+1)(2u+1)}{6},$$
so the whole scaled family is compressed into one arithmetic term rather than one loop over each individual pair.
The fourth-root bound for the outer parameter comes from the fact that the positive coordinate in each branch already grows like \(b^4\). In particular, the formulas imply
$$b\le \left\lfloor (4N)^{1/4}\right\rfloor.$$
Once \(b\) is fixed, the positive coordinate alone gives an upper bound for \(a\). For odd \(b\),
$$b(b+4a)^3\le N \quad\Longrightarrow\quad a\le \frac{\left\lfloor (N/b)^{1/3}\right\rfloor-b}{4}.$$
For even \(b\),
$$2b\left(\frac{b}{2}+2a\right)^3\le N \quad\Longrightarrow\quad a\le \frac{\left\lfloor (N/(2b))^{1/3}\right\rfloor-\frac{b}{2}}{2}.$$
These cube-root bounds drastically reduce the search space before the remaining filters are applied.
The implementations include a small checkpoint at \(N=1000\). Here
$$\left\lfloor(4N)^{1/4}\right\rfloor=\left\lfloor 4000^{1/4}\right\rfloor=7.$$
Only two primitive seeds survive all conditions and still satisfy the size bound.
First, with \(a=1\) and odd \(b=1\),
$$X_\ast=1(1+4)^3=125,\qquad Y_\ast=4\cdot 1\cdot(1-2)^3=-4.$$
The kernel condition is \( \operatorname{cf}(1)\ne \operatorname{cf}(4)\), so the seed is valid. Its scale limit is
$$t_{\max}=\left\lfloor\sqrt{\frac{1000}{125}}\right\rfloor=2,$$
so its total contribution is
$$ (125+4)(1^2+2^2)=129\cdot 5=645. $$
Second, with \(a=1\) and even \(b=2\),
$$X_\ast=2\cdot 2\left(1+2\right)^3=108,\qquad Y_\ast=1(1-4)^3=-27.$$
Now \( \operatorname{cf}(4)\ne \operatorname{cf}(1)\), and
$$t_{\max}=\left\lfloor\sqrt{\frac{1000}{108}}\right\rfloor=3.$$
This contributes
$$ (108+27)(1^2+2^2+3^2)=135\cdot 14=1890. $$
Adding the two surviving families gives
$$645+1890=2535,$$
which matches the checkpoint in the implementations exactly.
The implementation begins by computing exact integer fourth roots, cube roots, and square roots, so the search bounds are numerically safe even at very large \(N\).
It then builds a smallest-prime-factor sieve up to the limit needed for all cube-free-kernel lookups. From that sieve it derives the table of cube-free kernels \( \operatorname{cf}(n)\).
Next, it scans the outer parameter up to \( \lfloor(4N)^{1/4}\rfloor\). Odd values use the odd primitive formula; even values use the even primitive formula. For each branch it computes the cube-root upper bound for the inner parameter, then applies, in order, the coprimality test, the excluded diagonal test, the cube-free-kernel inequality, and the final max-norm check.
Whenever a primitive seed survives, the implementation computes \(t_{\max}\), evaluates the closed form for \(1^2+\cdots+t_{\max}^2\), multiplies by \(|X_\ast|+|Y_\ast|\), and adds the result modulo \(M\).
The C++ version parallelizes the outer-parameter loop across several threads, while the Python and Java implementations perform the same arithmetic sequentially. In all three languages, the mathematical work is the same.
Let \(B=\lfloor(4N)^{1/4}\rfloor\). The outer loop runs over \(O(B)=O(N^{1/4})\) values.
For a fixed outer parameter \(b\), the inner limit is \(O((N/b)^{1/3})\). Summing this over all \(b\le B\) gives the rough bound
$$\sum_{b\le B} O\!\left((N/b)^{1/3}\right)=O\!\left(N^{1/3}\sum_{b\le B} b^{-1/3}\right)=O(N^{1/2}).$$
So the arithmetic enumeration is about \(O(N^{1/2})\) candidate checks before constant-factor pruning from coprimality and kernel filters. The sieve for cube-free kernels reaches only the largest parameter size, which is \(O(N^{1/3})\), so its cost is \(O(N^{1/3}\log\log N)\) time and \(O(N^{1/3})\) memory.
In practice, the wall-clock time is better than a naive pair scan by many orders of magnitude, and the threaded C++ implementation further reduces elapsed time without changing the asymptotic count.
Die Implementierungen berechnen eine modulare Größe \(H(N)\) für \(N=10^{15}\) mit dem Modulus
$$M=1031^3+2.$$
Statt alle ganzzahligen Paare \((x,y)\) mit \(\max(|x|,|y|)\le N\) direkt zu durchsuchen, beginnen die Programme mit der algebraischen Klassifikation der für die verschachtelten Radikale relevanten Paare. Danach lautet die Aufgabe
$$H(N)=\sum_{\substack{(x,y)\ \text{zulässig}\\ \max(|x|,|y|)\le N}} (|x|+|y|)\pmod{M}.$$
Der entscheidende Punkt ist: Jedes zulässige Paar entsteht aus einem primitiven Parameterpaar und einem gemeinsamen quadratischen Skalierungsfaktor. Dadurch wird aus einer zweidimensionalen Gitter-Suche ein arithmetisches Enumerationsproblem.
Die drei Implementierungen verwenden zwei primitive Familien, je nachdem ob einer der Parameter gerade oder ungerade ist. Zur mathematischen Beschreibung seien \(a,b\in\mathbb{Z}_{>0}\) teilerfremd mit \(a\ne 2b\), und \(t\ge 1\) sei der gemeinsame Skalierungsfaktor.
Für ungerades \(b\) lautet das primitive Paar
$$X_{\mathrm{odd}}=b(b+4a)^3,\qquad Y_{\mathrm{odd}}=4a(a-2b)^3.$$
Für gerades \(b\) lautet es
$$X_{\mathrm{even}}=2b\left(\frac{b}{2}+2a\right)^3,\qquad Y_{\mathrm{even}}=a(a-2b)^3.$$
Genau diese beiden Zweige werden in C++, Python und Java ausgewertet. Die Bedingung
$$\gcd(a,b)=1$$
verhindert Mehrfachdarstellungen, und
$$a=2b$$
wird ausgeschlossen, weil dann ein kubischer Faktor verschwindet und der degenerierte Fall entsteht.
Aus jedem primitiven Paar entsteht die gesamte zugehörige Familie durch Multiplikation beider Koordinaten mit demselben Quadrat:
$$ (x,y)=\left(t^2X_\ast,\ t^2Y_\ast\right). $$
Dabei steht \((X_\ast,Y_\ast)\) für das primitive Paar aus dem passenden Paritätszweig.
Unter der Schranke \(\max(|x|,|y|)\le N\) ist damit
$$t^2\max(|X_\ast|,|Y_\ast|)\le N,$$
also
$$t_{\max}=\left\lfloor \sqrt{\frac{N}{\max(|X_\ast|,|Y_\ast|)}} \right\rfloor.$$
Deshalb muss das Programm nur primitive Startpaare erzeugen; alle skalierten Nachfolger werden anschließend in einem Schritt aufsummiert.
Die Implementierungen berechnen den würfelfreien Kern
$$\operatorname{cf}(n)=\prod_{p} p^{v_p(n)\bmod 3},$$
also den Teil, der nach Entfernung aller vollständigen Kuben übrig bleibt. Ein primitives Paar wird nur dann akzeptiert, wenn sich die beiden relevanten Kerne unterscheiden.
Für ungerades \(b\) gilt die Bedingung
$$\operatorname{cf}(b)\ne \operatorname{cf}(4a),$$
für gerades \(b\) dagegen
$$\operatorname{cf}(2b)\ne \operatorname{cf}(a).$$
Damit werden genau jene Fälle entfernt, in denen beide Kubikwurzeln auf denselben würfelfreien Anteil zurückfallen und die nichttriviale Radikalstruktur verschwindet.
Für einen festen primitiven Startwert trägt jedes skalierte Paar den Betrag
$$|t^2X_\ast|+|t^2Y_\ast|=t^2\bigl(|X_\ast|+|Y_\ast|\bigr)$$
bei. Somit ist der Gesamtbeitrag dieses Startwerts
$$\bigl(|X_\ast|+|Y_\ast|\bigr)\sum_{t=1}^{t_{\max}} t^2.$$
Verwendet wird die Standardformel
$$\sum_{t=1}^{u} t^2=\frac{u(u+1)(2u+1)}{6},$$
wodurch aus einer ganzen Folge zulässiger Skalierungen ein einziger Summand wird.
Der äußere Parameter reicht nur bis zur vierten Wurzel, weil die positive Koordinate in beiden Zweigen bereits wie \(b^4\) wächst. Daraus folgt
$$b\le \left\lfloor (4N)^{1/4}\right\rfloor.$$
Ist \(b\) fest, erhält man aus der positiven Koordinate eine Obergrenze für \(a\). Für ungerades \(b\) gilt
$$b(b+4a)^3\le N \quad\Longrightarrow\quad a\le \frac{\left\lfloor (N/b)^{1/3}\right\rfloor-b}{4},$$
und für gerades \(b\)
$$2b\left(\frac{b}{2}+2a\right)^3\le N \quad\Longrightarrow\quad a\le \frac{\left\lfloor (N/(2b))^{1/3}\right\rfloor-\frac{b}{2}}{2}.$$
Genau diese Kubikwurzelgrenzen verwenden die Programme, bevor sie die restlichen Tests anwenden.
Für \(N=1000\) ist
$$\left\lfloor(4N)^{1/4}\right\rfloor=\left\lfloor 4000^{1/4}\right\rfloor=7.$$
Nach allen Filtern bleiben genau zwei primitive Familien übrig.
Für \(a=1\) und ungerades \(b=1\) erhält man
$$X_\ast=1(1+4)^3=125,\qquad Y_\ast=4\cdot 1\cdot(1-2)^3=-4.$$
Da \( \operatorname{cf}(1)\ne \operatorname{cf}(4)\), ist dieses Paar zulässig. Ferner
$$t_{\max}=\left\lfloor\sqrt{\frac{1000}{125}}\right\rfloor=2,$$
also ist der Beitrag
$$ (125+4)(1^2+2^2)=129\cdot 5=645. $$
Für \(a=1\) und gerades \(b=2\) bekommt man
$$X_\ast=2\cdot 2(1+2)^3=108,\qquad Y_\ast=1(1-4)^3=-27.$$
Hier gilt \( \operatorname{cf}(4)\ne \operatorname{cf}(1)\), und
$$t_{\max}=\left\lfloor\sqrt{\frac{1000}{108}}\right\rfloor=3.$$
Daraus folgt der Beitrag
$$ (108+27)(1^2+2^2+3^2)=135\cdot 14=1890. $$
In Summe ergibt das
$$645+1890=2535,$$
genau der Kontrollwert der Implementierungen.
Zuerst berechnet die Implementierung exakte ganzzahlige vierte Wurzeln, Kubikwurzeln und Quadratwurzeln, damit die Suchgrenzen auch für sehr große \(N\) numerisch zuverlässig bleiben.
Danach wird ein Sieb der kleinsten Primteiler bis zur maximal benötigten Grenze aufgebaut. Aus diesem Sieb entsteht die Tabelle der würfelfreien Kerne \( \operatorname{cf}(n)\).
Anschließend läuft die äußere Schleife über alle Parameter bis \( \lfloor(4N)^{1/4}\rfloor\). Ungerade Werte verwenden den ungeraden Zweig, gerade Werte den geraden Zweig. Für jeden Fall wird erst die Kubikwurzelgrenze für den inneren Parameter bestimmt; dann folgen Teilerfremdheitstest, Ausschluss der Diagonale \(a=2b\), Vergleich der würfelfreien Kerne und schließlich die eigentliche Größenprüfung \(\max(|x|,|y|)\le N\).
Bleibt ein primitives Startpaar übrig, berechnet die Implementierung \(t_{\max}\), setzt die geschlossene Formel für \(1^2+\cdots+t_{\max}^2\) ein, multipliziert mit \(|X_\ast|+|Y_\ast|\) und addiert alles modulo \(M\).
Die C++-Version parallelisiert die äußere Schleife über mehrere Threads; die Python- und Java-Version führen dieselbe Arithmetik sequentiell aus. Mathematisch sind die drei Programme identisch.
Sei \(B=\lfloor(4N)^{1/4}\rfloor\). Dann hat die äußere Schleife \(O(B)=O(N^{1/4})\) Iterationen.
Für festes \(b\) ist die innere Grenze von der Größenordnung \(O((N/b)^{1/3})\). Summiert man dies über alle \(b\le B\), erhält man näherungsweise
$$\sum_{b\le B} O\!\left((N/b)^{1/3}\right)=O\!\left(N^{1/3}\sum_{b\le B} b^{-1/3}\right)=O(N^{1/2}).$$
Die eigentliche Enumeration liegt also ungefähr bei \(O(N^{1/2})\) Kandidatenprüfungen, bevor die konstanten Faktoren durch Teilerfremdheit und Kernfilter reduziert werden. Das Sieb für die würfelfreien Kerne reicht nur bis zur größten Parametergröße, also \(O(N^{1/3})\), und kostet \(O(N^{1/3}\log\log N)\) Zeit sowie \(O(N^{1/3})\) Speicher.
Praktisch ist das um Größenordnungen schneller als eine direkte Suche über alle Paare, und die parallelisierte C++-Version verkürzt zusätzlich die Laufzeit, ohne die asymptotische Struktur zu verändern.
Uygulamalar, \(N=10^{15}\) için şu modül altında bir \(H(N)\) değeri hesaplar:
$$M=1031^3+2.$$
\(\max(|x|,|y|)\le N\) koşulunu sağlayan tüm tamsayı çiftlerini doğrudan taramak yerine, çözüm önce problemdeki iç içe kök yapılarına karşılık gelen uygun çiftlerin cebirsel sınıflandırmasını kullanır. Bu noktadan sonra amaç
$$H(N)=\sum_{\substack{(x,y)\ \text{uygun}\\ \max(|x|,|y|)\le N}} (|x|+|y|)\pmod{M}$$
şeklinde yazılır. Kritik gözlem şudur: Her uygun çift, ilkel bir parametre çifti ile ortak bir kare ölçek çarpanından türetilir. Böylece iki boyutlu kaba tarama yerini aritmetik bir üretim yöntemine bırakır.
Üç uygulama da, parametrelerden birinin tek ya da çift olmasına göre iki ayrı ilkel aile kullanır. Kod içi isimleri tekrar etmeden, aralarında asal olan \(a,b\in\mathbb{Z}_{>0}\) ve \(a\ne 2b\) alalım; ayrıca ortak ölçek çarpanı \(t\ge 1\) olsun.
\(b\) tek ise ilkel çift
$$X_{\mathrm{odd}}=b(b+4a)^3,\qquad Y_{\mathrm{odd}}=4a(a-2b)^3$$
şeklindedir. \(b\) çift ise
$$X_{\mathrm{even}}=2b\left(\frac{b}{2}+2a\right)^3,\qquad Y_{\mathrm{even}}=a(a-2b)^3$$
olur. C++, Python ve Java çözümleri tam olarak bu iki dalı dolaşır. Ayrıca
$$\gcd(a,b)=1$$
koşulu yinelenen gösterimleri atar,
$$a=2b$$
doğrusu ise kübik çarpanlardan birini sıfırladığı için dışarıda bırakılır.
Her ilkel çift, her iki koordinatı da aynı kare ile çarparak sonsuz bir aile üretir:
$$ (x,y)=\left(t^2X_\ast,\ t^2Y_\ast\right). $$
Burada \((X_\ast,Y_\ast)\), tek ya da çift dala ait ilkel çifttir.
Kutu kısıtı altında
$$t^2\max(|X_\ast|,|Y_\ast|)\le N$$
olmalıdır; dolayısıyla izin verilen en büyük ölçek
$$t_{\max}=\left\lfloor \sqrt{\frac{N}{\max(|X_\ast|,|Y_\ast|)}} \right\rfloor$$
olur. Bu nedenle program yalnızca ilkel çekirdekleri üretir; ilkel olmayan tüm çözümler daha sonra tek bir kapalı toplamla hesaba katılır.
Uygulamalar şu küpsüz çekirdeği önceden hesaplar:
$$\operatorname{cf}(n)=\prod_{p} p^{v_p(n)\bmod 3}.$$
Yani her tam küp \(p^3\) çıkarılır; yalnızca üslerin \(1\) veya \(2\) kalıntıları korunur. Bir ilkel çift ancak iki ilgili çekirdek farklıysa kabul edilir.
\(b\) tek olduğunda gerekli eşitsizlik
$$\operatorname{cf}(b)\ne \operatorname{cf}(4a),$$
\(b\) çift olduğunda ise
$$\operatorname{cf}(2b)\ne \operatorname{cf}(a)$$
şeklindedir. Bu filtre, iki küpkök bileşeninin aynı küpsüz parçaya indirgenip özgün iç içe kök yapısını kaybettiği durumları eler.
Sabit bir ilkel çift için ölçeklenmiş her çözümün katkısı
$$|t^2X_\ast|+|t^2Y_\ast|=t^2\bigl(|X_\ast|+|Y_\ast|\bigr)$$
olur. O halde bu ilkel çekirdeğin toplam katkısı
$$\bigl(|X_\ast|+|Y_\ast|\bigr)\sum_{t=1}^{t_{\max}} t^2$$
şeklindedir. Kod burada standart kapalı formu kullanır:
$$\sum_{t=1}^{u} t^2=\frac{u(u+1)(2u+1)}{6}.$$
Böylece tek tek tüm ölçekleri dolaşmak yerine, bir ilkel ailenin toplamı tek aritmetik terime sıkıştırılır.
Dış parametre için dördüncü kök sınırı gelir; çünkü her iki dalda da pozitif koordinat en az \(b^4\) mertebesinde büyür. Buradan
$$b\le \left\lfloor (4N)^{1/4}\right\rfloor$$
elde edilir. \(b\) sabitlendikten sonra, pozitif koordinat tek başına \(a\) için bir üst sınır verir. Tek \(b\) için
$$b(b+4a)^3\le N \quad\Longrightarrow\quad a\le \frac{\left\lfloor (N/b)^{1/3}\right\rfloor-b}{4},$$
çift \(b\) için ise
$$2b\left(\frac{b}{2}+2a\right)^3\le N \quad\Longrightarrow\quad a\le \frac{\left\lfloor (N/(2b))^{1/3}\right\rfloor-\frac{b}{2}}{2}.$$
Uygulamalardaki iç döngü sınırları tam olarak bu küpkök formüllerinden gelir.
\(N=1000\) için
$$\left\lfloor(4N)^{1/4}\right\rfloor=\left\lfloor 4000^{1/4}\right\rfloor=7$$
olur. Bütün koşullar uygulandığında yalnızca iki ilkel aile kalır.
Birincisi \(a=1\), tek \(b=1\) durumudur:
$$X_\ast=1(1+4)^3=125,\qquad Y_\ast=4\cdot 1\cdot(1-2)^3=-4.$$
\(\operatorname{cf}(1)\ne \operatorname{cf}(4)\) olduğu için çift geçerlidir. Ölçek sınırı
$$t_{\max}=\left\lfloor\sqrt{\frac{1000}{125}}\right\rfloor=2$$
ve toplam katkı
$$ (125+4)(1^2+2^2)=129\cdot 5=645 $$
olur. İkincisi \(a=1\), çift \(b=2\) durumudur:
$$X_\ast=2\cdot 2(1+2)^3=108,\qquad Y_\ast=1(1-4)^3=-27.$$
Burada da \(\operatorname{cf}(4)\ne \operatorname{cf}(1)\) ve
$$t_{\max}=\left\lfloor\sqrt{\frac{1000}{108}}\right\rfloor=3.$$
Bu ailenin katkısı
$$ (108+27)(1^2+2^2+3^2)=135\cdot 14=1890 $$
eder. Sonuçta
$$645+1890=2535,$$
ki bu değer uygulamaların küçük kontrol sonucuyla birebir aynıdır.
Uygulama önce tam sayı dördüncü kök, küpkök ve karekök hesaplarını güvenli biçimde yapar; böylece çok büyük \(N\) değerlerinde yuvarlama hataları arama sınırlarını bozmaz.
Daha sonra gerekli üst sınıra kadar en küçük asal bölen süzgeci kurulur. Bu süzgeçten küpsüz çekirdek tablosu \( \operatorname{cf}(n)\) türetilir.
Ardından dış parametre \( \lfloor(4N)^{1/4}\rfloor\) değerine kadar dolaşılır. Tek değerlerde tek dal, çift değerlerde çift dal kullanılır. Her dal için önce iç parametrenin küpkök sınırı hesaplanır; sonra sırasıyla aralarında asallık testi, \(a=2b\) diyagonalinin dışlanması, küpsüz çekirdek eşitsizliği ve son olarak \(\max(|x|,|y|)\le N\) kontrolü uygulanır.
Koşulları geçen her ilkel aile için \(t_{\max}\) bulunur, \(1^2+\cdots+t_{\max}^2\) kapalı formdan hesaplanır, \(|X_\ast|+|Y_\ast|\) ile çarpılır ve sonuç \(M\) modunda toplanır.
C++ sürümü dış döngüyü birden çok iş parçacığına böler; Python ve Java sürümleri aynı aritmetiği tek iş parçacığında yürütür. Matematiksel algoritma üç dilde de aynıdır.
\(B=\lfloor(4N)^{1/4}\rfloor\) olsun. Dış döngü \(O(B)=O(N^{1/4})\) kez çalışır.
Sabit bir \(b\) için iç sınır \(O((N/b)^{1/3})\) mertebesindedir. Bunu tüm \(b\le B\) üzerinde toplarsak
$$\sum_{b\le B} O\!\left((N/b)^{1/3}\right)=O\!\left(N^{1/3}\sum_{b\le B} b^{-1/3}\right)=O(N^{1/2})$$
elde edilir. Yani sabit çarpanlı filtrelerden önce aday kontrol sayısı kabaca \(O(N^{1/2})\) düzeyindedir. Küpsüz çekirdek süzgeci ise yalnızca en büyük parametre ölçeğine kadar gider; bu da \(O(N^{1/3})\) bellek ve \(O(N^{1/3}\log\log N)\) zaman gerektirir.
Pratikte bu yöntem, tüm \((x,y)\) çiftlerini kaba kuvvetle taramaktan çok daha hızlıdır. İş parçacıklı C++ sürümü de asimptotik sınıfı değiştirmeden duvar süresini daha da düşürür.
Las implementaciones calculan una cantidad modular \(H(N)\) para \(N=10^{15}\), con módulo
$$M=1031^3+2.$$
En vez de recorrer directamente todos los pares enteros \((x,y)\) con \(\max(|x|,|y|)\le N\), el método parte de la clasificación algebraica de los pares admisibles asociados a las identidades de radicales anidados del problema. A partir de ahí, la tarea es
$$H(N)=\sum_{\substack{(x,y)\ \text{admisible}\\ \max(|x|,|y|)\le N}} (|x|+|y|)\pmod{M}.$$
La observación decisiva es que cada par admisible nace de un par primitivo de parámetros y de un factor de escala cuadrado común. Eso permite reemplazar una exploración geométrica por una enumeración puramente aritmética.
Las tres implementaciones usan dos familias primitivas, separadas por la paridad de uno de los parámetros. Para describirlas sin reutilizar los nombres del código, tomemos \(a,b\in\mathbb{Z}_{>0}\) coprimos con \(a\ne 2b\), y un factor de escala común \(t\ge 1\).
Si \(b\) es impar, el par primitivo es
$$X_{\mathrm{odd}}=b(b+4a)^3,\qquad Y_{\mathrm{odd}}=4a(a-2b)^3.$$
Si \(b\) es par, el par primitivo es
$$X_{\mathrm{even}}=2b\left(\frac{b}{2}+2a\right)^3,\qquad Y_{\mathrm{even}}=a(a-2b)^3.$$
Esas son exactamente las dos ramas que recorren las versiones en C++, Python y Java. La condición
$$\gcd(a,b)=1$$
elimina representaciones redundantes, mientras que
$$a=2b$$
se descarta porque anula uno de los factores cúbicos y produce el caso degenerado.
Cada semilla primitiva produce una familia infinita multiplicando ambas coordenadas por el mismo cuadrado:
$$ (x,y)=\left(t^2X_\ast,\ t^2Y_\ast\right). $$
Aquí \((X_\ast,Y_\ast)\) denota el par primitivo apropiado, según la rama par o impar.
Bajo la restricción \(\max(|x|,|y|)\le N\), debe cumplirse
$$t^2\max(|X_\ast|,|Y_\ast|)\le N,$$
de donde sale
$$t_{\max}=\left\lfloor \sqrt{\frac{N}{\max(|X_\ast|,|Y_\ast|)}} \right\rfloor.$$
Por eso el algoritmo solo enumera semillas primitivas y luego absorbe todos los múltiplos cuadrados en una única fórmula cerrada.
Las implementaciones precalculan el núcleo libre de cubos
$$\operatorname{cf}(n)=\prod_{p} p^{v_p(n)\bmod 3},$$
es decir, el resultado de quitar todos los cubos completos y conservar solamente los exponentes residuales \(1\) o \(2\). Una semilla se conserva solo si los dos núcleos relevantes son distintos.
Si \(b\) es impar, la condición es
$$\operatorname{cf}(b)\ne \operatorname{cf}(4a).$$
Si \(b\) es par, la condición pasa a ser
$$\operatorname{cf}(2b)\ne \operatorname{cf}(a).$$
Este filtro descarta los casos en los que las dos componentes cúbicas comparten la misma parte libre de cubos y la estructura no trivial del radical anidado se pierde.
Para una semilla fija, cada solución escalada aporta
$$|t^2X_\ast|+|t^2Y_\ast|=t^2\bigl(|X_\ast|+|Y_\ast|\bigr).$$
Así, la contribución total de esa familia primitiva es
$$\bigl(|X_\ast|+|Y_\ast|\bigr)\sum_{t=1}^{t_{\max}} t^2.$$
El código usa la fórmula cerrada clásica
$$\sum_{t=1}^{u} t^2=\frac{u(u+1)(2u+1)}{6},$$
de modo que una familia completa de pares válidos se condensa en un solo término aritmético.
La cota exterior de cuarta raíz aparece porque la coordenada positiva en ambas ramas ya crece como \(b^4\). En consecuencia,
$$b\le \left\lfloor (4N)^{1/4}\right\rfloor.$$
Una vez fijado \(b\), la coordenada positiva también impone una cota para \(a\). Para \(b\) impar,
$$b(b+4a)^3\le N \quad\Longrightarrow\quad a\le \frac{\left\lfloor (N/b)^{1/3}\right\rfloor-b}{4},$$
y para \(b\) par,
$$2b\left(\frac{b}{2}+2a\right)^3\le N \quad\Longrightarrow\quad a\le \frac{\left\lfloor (N/(2b))^{1/3}\right\rfloor-\frac{b}{2}}{2}.$$
Estas son exactamente las cotas de raíz cúbica que usan las implementaciones antes de aplicar los demás filtros.
Para \(N=1000\),
$$\left\lfloor(4N)^{1/4}\right\rfloor=\left\lfloor 4000^{1/4}\right\rfloor=7.$$
Después de todas las condiciones sobreviven solo dos familias primitivas.
La primera aparece con \(a=1\) y \(b=1\) impar:
$$X_\ast=1(1+4)^3=125,\qquad Y_\ast=4\cdot 1\cdot(1-2)^3=-4.$$
Como \( \operatorname{cf}(1)\ne \operatorname{cf}(4)\), la semilla es válida. Además,
$$t_{\max}=\left\lfloor\sqrt{\frac{1000}{125}}\right\rfloor=2,$$
y su contribución es
$$ (125+4)(1^2+2^2)=129\cdot 5=645. $$
La segunda aparece con \(a=1\) y \(b=2\) par:
$$X_\ast=2\cdot 2(1+2)^3=108,\qquad Y_\ast=1(1-4)^3=-27.$$
Aquí también \( \operatorname{cf}(4)\ne \operatorname{cf}(1)\), y
$$t_{\max}=\left\lfloor\sqrt{\frac{1000}{108}}\right\rfloor=3.$$
Su contribución total es
$$ (108+27)(1^2+2^2+3^2)=135\cdot 14=1890. $$
En consecuencia,
$$645+1890=2535,$$
que coincide exactamente con el valor de comprobación de las implementaciones.
La implementación empieza calculando raíces cuarta, cúbica y cuadrada enteras exactas, para que las cotas de búsqueda sean seguras incluso con \(N\) muy grande.
Luego construye una criba de menor factor primo hasta el límite necesario para todas las consultas del núcleo libre de cubos. A partir de esa criba deriva la tabla \( \operatorname{cf}(n)\).
Después recorre el parámetro exterior hasta \( \lfloor(4N)^{1/4}\rfloor\). Los valores impares usan la rama impar y los pares la rama par. En cada caso calcula primero la cota cúbica del parámetro interior; luego aplica, en este orden, la prueba de coprimalidad, la exclusión de la diagonal \(a=2b\), la desigualdad entre núcleos libres de cubos y por último la comprobación de \(\max(|x|,|y|)\le N\).
Cuando una semilla sobrevive, la implementación calcula \(t_{\max}\), evalúa la fórmula cerrada de \(1^2+\cdots+t_{\max}^2\), multiplica por \(|X_\ast|+|Y_\ast|\) y acumula todo módulo \(M\).
La versión en C++ paraleliza el bucle exterior entre varios hilos; las versiones en Python y Java ejecutan la misma aritmética de forma secuencial. El contenido matemático es el mismo en los tres casos.
Sea \(B=\lfloor(4N)^{1/4}\rfloor\). Entonces el bucle exterior tiene \(O(B)=O(N^{1/4})\) iteraciones.
Para un \(b\) fijo, la cota interior es \(O((N/b)^{1/3})\). Sumando sobre todos los \(b\le B\), se obtiene
$$\sum_{b\le B} O\!\left((N/b)^{1/3}\right)=O\!\left(N^{1/3}\sum_{b\le B} b^{-1/3}\right)=O(N^{1/2}).$$
Por tanto, la enumeración aritmética está aproximadamente en \(O(N^{1/2})\) comprobaciones de candidatos antes de la poda por coprimalidad y por núcleos. La criba para los núcleos libres de cubos llega solo hasta la mayor escala de parámetros, es decir, \(O(N^{1/3})\), y cuesta \(O(N^{1/3}\log\log N)\) tiempo y \(O(N^{1/3})\) memoria.
En la práctica esto es muchísimo más rápido que una exploración ingenua de todos los pares \((x,y)\), y la paralelización en C++ reduce aún más el tiempo real sin modificar la complejidad asintótica.
这些实现要求计算一个模量下的函数 \(H(N)\)，其中 \(N=10^{15}\)，模数为
$$M=1031^3+2.$$
程序并不直接枚举所有满足 \(\max(|x|,|y|)\le N\) 的整数对 \((x,y)\)。相反，它从题目对应的嵌套根式恒等式所导出的“可行整数对分类”出发。于是问题被改写成
$$H(N)=\sum_{\substack{(x,y)\ \text{可行}\\ \max(|x|,|y|)\le N}} (|x|+|y|)\pmod{M}.$$
真正关键的事实是：每一个可行整数对都可以唯一地看成“一个原始参数种子”再乘上“一个公共平方尺度”。因此，代码不需要在平面格点上暴力搜索，而是可以通过数论参数化直接生成全部候选。
三种语言的实现都使用两族原始公式，它们只由某个参数的奇偶性区分。为了避免直接重复源码里的局部命名，这里记互素正整数参数为 \(a,b\)，并要求 \(a\ne 2b\)；再记公共缩放因子为 \(t\ge 1\)。
当 \(b\) 为奇数时，原始整数对写成
$$X_{\mathrm{odd}}=b(b+4a)^3,\qquad Y_{\mathrm{odd}}=4a(a-2b)^3.$$
当 \(b\) 为偶数时，原始整数对写成
$$X_{\mathrm{even}}=2b\left(\frac{b}{2}+2a\right)^3,\qquad Y_{\mathrm{even}}=a(a-2b)^3.$$
这正是 C++、Python、Java 三个实现真正遍历的两条分支。条件
$$\gcd(a,b)=1$$
用于去掉重复表示，而
$$a=2b$$
这条直线被排除，因为此时一个三次因子会直接变成零，落入退化情形。
每一个原始种子都会生成一整族解，其形式是同时把两个坐标乘以同一个平方：
$$ (x,y)=\left(t^2X_\ast,\ t^2Y_\ast\right). $$
这里的 \((X_\ast,Y_\ast)\) 表示上面奇分支或偶分支对应的原始整数对。
要满足边界 \(\max(|x|,|y|)\le N\)，就必须有
$$t^2\max(|X_\ast|,|Y_\ast|)\le N,$$
因此最大允许缩放量为
$$t_{\max}=\left\lfloor \sqrt{\frac{N}{\max(|X_\ast|,|Y_\ast|)}} \right\rfloor.$$
这一步非常重要，因为它说明代码只要枚举“原始种子”，再用一个平方和公式把所有非原始放大版本一次性吸收掉即可。
实现会先预处理所谓的立方自由核
$$\operatorname{cf}(n)=\prod_{p} p^{v_p(n)\bmod 3}.$$
含义是：把所有完整的三次幂因子 \(p^3\) 都剥掉，只保留指数对 \(3\) 的余数 \(1\) 或 \(2\)。一个原始种子只有在两边相关的立方自由核不同的时候才会被保留。
若 \(b\) 为奇数，要求
$$\operatorname{cf}(b)\ne \operatorname{cf}(4a).$$
若 \(b\) 为偶数，要求
$$\operatorname{cf}(2b)\ne \operatorname{cf}(a).$$
这个条件的作用，是去掉那种“两项立方根最后落到同一个立方自由部分上，从而失去非平凡嵌套结构”的候选。
对固定的原始种子而言，任意一个缩放后的解对总和的贡献为
$$|t^2X_\ast|+|t^2Y_\ast|=t^2\bigl(|X_\ast|+|Y_\ast|\bigr).$$
所以这一整族的总贡献就是
$$\bigl(|X_\ast|+|Y_\ast|\bigr)\sum_{t=1}^{t_{\max}} t^2.$$
实现直接使用平方和闭式
$$\sum_{t=1}^{u} t^2=\frac{u(u+1)(2u+1)}{6},$$
这样就不必对每一个具体的缩放值分别累加，而是把整族解压缩成一个公式项。
外层参数只有到四次根量级，因为两条分支里的正坐标都至少按 \(b^4\) 级别增长。因此有
$$b\le \left\lfloor (4N)^{1/4}\right\rfloor.$$
在固定 \(b\) 后，正坐标本身又给出了 \(a\) 的上界。对奇数 \(b\)，
$$b(b+4a)^3\le N \quad\Longrightarrow\quad a\le \frac{\left\lfloor (N/b)^{1/3}\right\rfloor-b}{4};$$
对偶数 \(b\)，
$$2b\left(\frac{b}{2}+2a\right)^3\le N \quad\Longrightarrow\quad a\le \frac{\left\lfloor (N/(2b))^{1/3}\right\rfloor-\frac{b}{2}}{2}.$$
这正是实现里内层循环的立方根上界来源。实际程序随后还会继续做互素性、立方自由核以及最终大小约束的筛选。
对于 \(N=1000\)，有
$$\left\lfloor(4N)^{1/4}\right\rfloor=\left\lfloor 4000^{1/4}\right\rfloor=7.$$
把所有条件都套进去之后，只剩下两个原始家族真正有贡献。
第一个来自 \(a=1\)、奇数 \(b=1\)：
$$X_\ast=1(1+4)^3=125,\qquad Y_\ast=4\cdot 1\cdot(1-2)^3=-4.$$
此时 \( \operatorname{cf}(1)\ne \operatorname{cf}(4)\)，所以这一组保留。其最大缩放为
$$t_{\max}=\left\lfloor\sqrt{\frac{1000}{125}}\right\rfloor=2,$$
于是贡献为
$$ (125+4)(1^2+2^2)=129\cdot 5=645. $$
第二个来自 \(a=1\)、偶数 \(b=2\)：
$$X_\ast=2\cdot 2(1+2)^3=108,\qquad Y_\ast=1(1-4)^3=-27.$$
这里同样有 \( \operatorname{cf}(4)\ne \operatorname{cf}(1)\)，并且
$$t_{\max}=\left\lfloor\sqrt{\frac{1000}{108}}\right\rfloor=3.$$
这一族的贡献为
$$ (108+27)(1^2+2^2+3^2)=135\cdot 14=1890. $$
把两族相加得到
$$645+1890=2535,$$
这与实现中使用的小规模校验值完全一致。这个例子也清楚展示了“先找原始种子，再一次性累加全部平方缩放”的工作方式。
实现首先计算精确的整数四次根、立方根和平方根，这样即使 \(N\) 非常大，搜索上界也不会因为浮点舍入而出错。
接着程序建立最小素因子筛，并据此生成立方自由核表 \( \operatorname{cf}(n)\)。这个预处理只需要覆盖所有参数可能达到的最大范围。
之后，外层循环遍历到 \( \lfloor(4N)^{1/4}\rfloor\)。奇数走奇分支，偶数走偶分支。对每一个外层参数，程序先根据正坐标公式算出内层参数的立方根上界；再依次做互素性测试、排除退化直线 \(a=2b\)、比较立方自由核，并最终检查 \(\max(|x|,|y|)\le N\)。
一旦某个原始种子通过所有条件，程序就计算 \(t_{\max}\)，用闭式求出 \(1^2+\cdots+t_{\max}^2\)，再乘上 \(|X_\ast|+|Y_\ast|\)，并对模数 \(M\) 累加。
C++ 版本把外层循环分配到多个线程上；Python 与 Java 版本顺序执行同样的算术逻辑。三者的数学内容完全一致，只是运行方式不同。
设 \(B=\lfloor(4N)^{1/4}\rfloor\)。那么外层循环一共运行 \(O(B)=O(N^{1/4})\) 次。
对固定的 \(b\)，内层候选数量是 \(O((N/b)^{1/3})\)。把它对所有 \(b\le B\) 求和，可得
$$\sum_{b\le B} O\!\left((N/b)^{1/3}\right)=O\!\left(N^{1/3}\sum_{b\le B} b^{-1/3}\right)=O(N^{1/2}).$$
因此，在互素过滤与立方自由核过滤之前，参数枚举的粗略规模约为 \(O(N^{1/2})\) 次候选检查。立方自由核所需的筛表只扩展到最大参数尺度，也就是 \(O(N^{1/3})\)；因此它的预处理代价是 \(O(N^{1/3}\log\log N)\) 时间和 \(O(N^{1/3})\) 空间。
和直接枚举所有 \((x,y)\) 相比，这个复杂度已经低了很多个数量级。多线程 C++ 实现还能进一步缩短实际运行时间，但不会改变渐近复杂度。
Реализации вычисляют величину \(H(N)\) по модулю для \(N=10^{15}\), где
$$M=1031^3+2.$$
Вместо прямого перебора всех целочисленных пар \((x,y)\) с условием \(\max(|x|,|y|)\le N\), программа использует алгебраическую классификацию допустимых пар, связанных с тождествами для вложенных радикалов. После этого задача принимает вид
$$H(N)=\sum_{\substack{(x,y)\ \text{допустима}\\ \max(|x|,|y|)\le N}} (|x|+|y|)\pmod{M}.$$
Главное наблюдение состоит в том, что каждая допустимая пара возникает из примитивной параметризации и общего квадратного масштаба. Поэтому вместо наивного сканирования по всей плоскости можно перечислять только арифметически описанные семейства.
Во всех трех версиях используются две примитивные семьи, различающиеся четностью одного из параметров. Чтобы не повторять локальные имена из исходников, обозначим взаимно простые положительные параметры через \(a\) и \(b\), потребуем \(a\ne 2b\), а общий масштаб обозначим через \(t\ge 1\).
Если \(b\) нечетно, примитивная пара имеет вид
$$X_{\mathrm{odd}}=b(b+4a)^3,\qquad Y_{\mathrm{odd}}=4a(a-2b)^3.$$
Если \(b\) четно, то
$$X_{\mathrm{even}}=2b\left(\frac{b}{2}+2a\right)^3,\qquad Y_{\mathrm{even}}=a(a-2b)^3.$$
Именно эти две формулы проверяют реализации на C++, Python и Java. Условие
$$\gcd(a,b)=1$$
убирает повторяющиеся представления, а прямая
$$a=2b$$
исключается, поскольку тогда один из кубических множителей обращается в нуль и возникает вырожденный случай.
Каждая примитивная пара порождает бесконечное семейство, если умножить обе координаты на один и тот же квадрат:
$$ (x,y)=\left(t^2X_\ast,\ t^2Y_\ast\right). $$
Здесь \((X_\ast,Y_\ast)\) означает нужную примитивную пару из четной или нечетной ветви.
Из ограничения \(\max(|x|,|y|)\le N\) следует
$$t^2\max(|X_\ast|,|Y_\ast|)\le N,$$
поэтому максимальный допустимый масштаб равен
$$t_{\max}=\left\lfloor \sqrt{\frac{N}{\max(|X_\ast|,|Y_\ast|)}} \right\rfloor.$$
Отсюда видно, почему программа перечисляет только примитивные начальные пары: все их квадратные кратные затем суммируются одной формулой.
Реализации заранее вычисляют куб-свободное ядро
$$\operatorname{cf}(n)=\prod_{p} p^{v_p(n)\bmod 3}.$$
Это означает, что из числа удаляются все полные кубы, а показатели \(1\) и \(2\) по модулю \(3\) сохраняются. Примитивная пара принимается только тогда, когда два соответствующих ядра различны.
Для нечетного \(b\) условие таково:
$$\operatorname{cf}(b)\ne \operatorname{cf}(4a).$$
Для четного \(b\) оно меняется на
$$\operatorname{cf}(2b)\ne \operatorname{cf}(a).$$
Такой фильтр убирает случаи, где обе кубические компоненты имеют одну и ту же куб-свободную часть и нужная нетривиальная структура вложенного радикала исчезает.
Для фиксированного примитивного зерна вклад любой масштабированной пары равен
$$|t^2X_\ast|+|t^2Y_\ast|=t^2\bigl(|X_\ast|+|Y_\ast|\bigr).$$
Следовательно, полный вклад этой семьи равен
$$\bigl(|X_\ast|+|Y_\ast|\bigr)\sum_{t=1}^{t_{\max}} t^2.$$
В программах используется стандартная формула
$$\sum_{t=1}^{u} t^2=\frac{u(u+1)(2u+1)}{6},$$
поэтому вместо цикла по каждому масштабу семья целиком сворачивается в один арифметический член.
Внешний параметр ограничен четвертой степенной корнем, поскольку положительная координата в обеих ветвях растет уже как \(b^4\). Значит,
$$b\le \left\lfloor (4N)^{1/4}\right\rfloor.$$
После фиксации \(b\) сама положительная координата ограничивает и \(a\). Для нечетного \(b\)
$$b(b+4a)^3\le N \quad\Longrightarrow\quad a\le \frac{\left\lfloor (N/b)^{1/3}\right\rfloor-b}{4},$$
а для четного \(b\)
$$2b\left(\frac{b}{2}+2a\right)^3\le N \quad\Longrightarrow\quad a\le \frac{\left\lfloor (N/(2b))^{1/3}\right\rfloor-\frac{b}{2}}{2}.$$
Именно эти кубические границы используются в коде до того, как включаются проверки взаимной простоты, ядер и окончательного ограничения по размеру.
При \(N=1000\)
$$\left\lfloor(4N)^{1/4}\right\rfloor=\left\lfloor 4000^{1/4}\right\rfloor=7.$$
После всех условий остаются только две примитивные семьи.
Первая получается при \(a=1\) и нечетном \(b=1\):
$$X_\ast=1(1+4)^3=125,\qquad Y_\ast=4\cdot 1\cdot(1-2)^3=-4.$$
Так как \( \operatorname{cf}(1)\ne \operatorname{cf}(4)\), эта семья допустима. При этом
$$t_{\max}=\left\lfloor\sqrt{\frac{1000}{125}}\right\rfloor=2,$$
и ее вклад равен
$$ (125+4)(1^2+2^2)=129\cdot 5=645. $$
Вторая семья возникает при \(a=1\) и четном \(b=2\):
$$X_\ast=2\cdot 2(1+2)^3=108,\qquad Y_\ast=1(1-4)^3=-27.$$
Здесь снова \( \operatorname{cf}(4)\ne \operatorname{cf}(1)\), а
$$t_{\max}=\left\lfloor\sqrt{\frac{1000}{108}}\right\rfloor=3.$$
Ее вклад равен
$$ (108+27)(1^2+2^2+3^2)=135\cdot 14=1890. $$
Сумма дает
$$645+1890=2535,$$
что точно совпадает с контрольным значением в реализациях.
Сначала программа вычисляет точные целочисленные четвертые корни, кубические корни и квадратные корни, чтобы границы поиска были надежны даже при очень больших \(N\).
Затем строится решето наименьших простых делителей до максимальной нужной границы. По нему формируется таблица куб-свободных ядер \( \operatorname{cf}(n)\).
После этого внешний параметр перебирается до \( \lfloor(4N)^{1/4}\rfloor\). Для нечетных значений используется нечетная формула, для четных - четная. Внутри каждой ветви сначала вычисляется кубическая граница для второго параметра, затем последовательно применяются проверка взаимной простоты, исключение прямой \(a=2b\), сравнение куб-свободных ядер и окончательная проверка \(\max(|x|,|y|)\le N\).
Если примитивное зерно прошло все фильтры, программа находит \(t_{\max}\), вычисляет по закрытой формуле сумму \(1^2+\cdots+t_{\max}^2\), умножает ее на \(|X_\ast|+|Y_\ast|\) и прибавляет результат по модулю \(M\).
В версии на C++ внешний цикл распараллелен между потоками. Версии на Python и Java выполняют ту же арифметику последовательно. С математической точки зрения алгоритм везде одинаков.
Пусть \(B=\lfloor(4N)^{1/4}\rfloor\). Тогда внешний цикл имеет \(O(B)=O(N^{1/4})\) шагов.
Для фиксированного \(b\) число внутренних кандидатов имеет порядок \(O((N/b)^{1/3})\). После суммирования по всем \(b\le B\) получаем
$$\sum_{b\le B} O\!\left((N/b)^{1/3}\right)=O\!\left(N^{1/3}\sum_{b\le B} b^{-1/3}\right)=O(N^{1/2}).$$
Значит, арифметическая часть дает примерно \(O(N^{1/2})\) проверок кандидатов до учета сокращения за счет взаимной простоты и фильтра ядер. Решето для куб-свободных ядер доходит только до максимального размера параметров, то есть до \(O(N^{1/3})\), и требует \(O(N^{1/3}\log\log N)\) времени и \(O(N^{1/3})\) памяти.
На практике это на много порядков быстрее, чем прямой перебор всех \((x,y)\). Многопоточная версия на C++ уменьшает реальное время выполнения, не меняя асимптотического порядка.
تحسب هذه الحلول الكمية \(H(N)\) بترديد معيّن عند \(N=10^{15}\)، حيث
$$M=1031^3+2.$$
وبدلًا من تعداد جميع الأزواج الصحيحة \((x,y)\) التي تحقق \(\max(|x|,|y|)\le N\)، تعتمد الخوارزمية على تصنيف جبري للأزواج المقبولة المرتبطة بهويات الجذور المتداخلة في المسألة. وبعد هذا الاختزال تصبح المهمة
$$H(N)=\sum_{\substack{(x,y)\in \mathcal{A}\\ \max(|x|,|y|)\le N}} (|x|+|y|)\pmod{M}.$$
الفكرة الفاصلة هي أن كل زوج مقبول ينتج من بذرة أولية للمعلمات مع عامل تكبير مربع مشترك. لهذا تتحول المسألة من مسح مباشر لشبكة ثنائية الأبعاد إلى توليد حسابي منظم.
تستعمل التطبيقات الثلاثة عائلتين أوليتين، ويحددهما فقط كون أحد المعاملات زوجيًا أو فرديًا. ومن دون إعادة استعمال أسماء المتغيرات المحلية من الشيفرة، نرمز إلى المعاملين المتوافقين نسبيًا بـ \(a,b\in\mathbb{Z}_{>0}\) مع الشرط \(a\ne 2b\)، ونرمز إلى عامل القياس المشترك بـ \(t\ge 1\).
إذا كان \(b\) فرديًا فإن الزوج الأولي هو
$$X_{\mathrm{odd}}=b(b+4a)^3,\qquad Y_{\mathrm{odd}}=4a(a-2b)^3.$$
وإذا كان \(b\) زوجيًا فإن الزوج الأولي هو
$$X_{\mathrm{even}}=2b\left(\frac{b}{2}+2a\right)^3,\qquad Y_{\mathrm{even}}=a(a-2b)^3.$$
وهذان هما الفرعان اللذان تدور عليهما تطبيقات C++ وPython وJava. كما أن الشرط
$$\gcd(a,b)=1$$
يزيل التمثيلات المكررة، بينما يُستبعد الخط
$$a=2b$$
لأن أحد العوامل التكعيبية يصبح صفرًا فيه، فتظهر الحالة المتدهورة.
كل بذرة أولية تعطي عائلة لا نهائية إذا ضربنا الإحداثيين في المربع نفسه:
$$ (x,y)=\left(t^2X_\ast,\ t^2Y_\ast\right). $$
ويمثل \((X_\ast,Y_\ast)\) الزوج الأولي المناسب من الفرع الفردي أو الزوجي.
ومن شرط الصندوق \(\max(|x|,|y|)\le N\) نحصل على
$$t^2\max(|X_\ast|,|Y_\ast|)\le N,$$
ومن ثم يكون أكبر عامل مسموح هو
$$t_{\max}=\left\lfloor \sqrt{\frac{N}{\max(|X_\ast|,|Y_\ast|)}} \right\rfloor.$$
وهذا يفسر لماذا تكتفي الخوارزمية بتعداد البذور الأولية فقط، ثم تجمع جميع النسخ الموسعة دفعة واحدة باستخدام صيغة مغلقة.
تُحسب مسبقًا النواة الخالية من المكعبات
$$\operatorname{cf}(n)=\prod_{p} p^{v_p(n)\bmod 3}.$$
ومعناها إزالة كل العوامل التكعيبية الكاملة والإبقاء فقط على الباقي من الأسس modulo \(3\)، أي \(1\) أو \(2\). ولا تُقبل البذرة الأولية إلا إذا اختلفت النواتان المرتبطتان بطرفيها.
عندما يكون \(b\) فرديًا يكون الشرط
$$\operatorname{cf}(b)\ne \operatorname{cf}(4a),$$
وعندما يكون \(b\) زوجيًا يصبح
$$\operatorname{cf}(2b)\ne \operatorname{cf}(a).$$
هذا الفلتر يستبعد الحالات التي تسقط فيها المركبتان التكعيبيتان على الجزء نفسه الخالي من المكعبات، وعندها تختفي البنية غير التافهة للجذر المتداخل.
إذا ثبتت بذرة أولية واحدة فإن مساهمة أي نسخة موسعة منها تساوي
$$|t^2X_\ast|+|t^2Y_\ast|=t^2\bigl(|X_\ast|+|Y_\ast|\bigr).$$
إذن المساهمة الكلية لتلك العائلة هي
$$\bigl(|X_\ast|+|Y_\ast|\bigr)\sum_{t=1}^{t_{\max}} t^2.$$
ولهذا تستخدم الحلول الصيغة المعروفة
$$\sum_{t=1}^{u} t^2=\frac{u(u+1)(2u+1)}{6},$$
بحيث تتحول عائلة كاملة من الأزواج إلى حد حسابي واحد بدل المرور على كل قيمة من قيم \(t\) على حدة.
المعامل الخارجي لا يحتاج إلى أكثر من رتبة الجذر الرابع، لأن الإحداثي الموجب في كلا الفرعين ينمو أصلًا مثل \(b^4\). ولذلك
$$b\le \left\lfloor (4N)^{1/4}\right\rfloor.$$
وبعد تثبيت \(b\)، يعطي الإحداثي الموجب نفسه حدًا أعلى للمعامل \(a\). ففي الحالة الفردية
$$b(b+4a)^3\le N \quad\Longrightarrow\quad a\le \frac{\left\lfloor (N/b)^{1/3}\right\rfloor-b}{4},$$
وفي الحالة الزوجية
$$2b\left(\frac{b}{2}+2a\right)^3\le N \quad\Longrightarrow\quad a\le \frac{\left\lfloor (N/(2b))^{1/3}\right\rfloor-\frac{b}{2}}{2}.$$
وهذه هي حدود الجذر التكعيبي نفسها التي تستخدمها التطبيقات قبل بقية الاختبارات.
عند \(N=1000\) يكون
$$\left\lfloor(4N)^{1/4}\right\rfloor=\left\lfloor 4000^{1/4}\right\rfloor=7.$$
وبعد تطبيق جميع الشروط لا تبقى إلا عائلتان أوليتان لهما مساهمة فعلية.
الأولى عند \(a=1\) و\(b=1\) الفردي:
$$X_\ast=1(1+4)^3=125,\qquad Y_\ast=4\cdot 1\cdot(1-2)^3=-4.$$
وبما أن \( \operatorname{cf}(1)\ne \operatorname{cf}(4)\) فهذه البذرة صالحة. كما أن
$$t_{\max}=\left\lfloor\sqrt{\frac{1000}{125}}\right\rfloor=2,$$
فتكون مساهمتها
$$ (125+4)(1^2+2^2)=129\cdot 5=645. $$
والثانية عند \(a=1\) و\(b=2\) الزوجي:
$$X_\ast=2\cdot 2(1+2)^3=108,\qquad Y_\ast=1(1-4)^3=-27.$$
وهنا أيضًا \( \operatorname{cf}(4)\ne \operatorname{cf}(1)\)، ومع
$$t_{\max}=\left\lfloor\sqrt{\frac{1000}{108}}\right\rfloor=3$$
تصبح مساهمة هذه العائلة
$$ (108+27)(1^2+2^2+3^2)=135\cdot 14=1890. $$
ومن ثم نحصل على
$$645+1890=2535,$$
وهو بالضبط مقدار التحقق الصغير الذي تعطيه التطبيقات.
تبدأ الخوارزمية بحساب الجذر الرابع الصحيح والجذر التكعيبي الصحيح والجذر التربيعي الصحيح بدقة، حتى لا تسبب التقريبات العائمة أخطاء في حدود البحث عند القيم الكبيرة جدًا لـ \(N\).
ثم تُبنى غربلة لأصغر عامل أولي حتى الحد اللازم، ومنها تُشتق لاحقًا قيم النواة الخالية من المكعبات \( \operatorname{cf}(n)\).
بعد ذلك تدور الحلقة الخارجية حتى \( \lfloor(4N)^{1/4}\rfloor\). القيم الفردية تستخدم الفرع الفردي، والقيم الزوجية تستخدم الفرع الزوجي. وفي كل فرع يُحسب أولًا الحد التكعيبي للمعامل الداخلي، ثم تُطبق على الترتيب: مفحوصة التوافق النسبي، واستبعاد الخط \(a=2b\)، ومقارنة النواتين الخاليتين من المكعبات، وأخيرًا شرط \(\max(|x|,|y|)\le N\).
إذا نجت البذرة الأولية من كل هذه الاختبارات، تحسب الخوارزمية \(t_{\max}\)، ثم قيمة \(1^2+\cdots+t_{\max}^2\) بصيغة مغلقة، وتضربها في \(|X_\ast|+|Y_\ast|\)، ثم تجمع الناتج بترديد \(M\).
نسخة C++ توزع الحلقة الخارجية على عدة خيوط تنفيذ، بينما تنفذ نسختا Python وJava الحساب نفسه بصورة متسلسلة. أما المحتوى الرياضي فهو واحد في اللغات الثلاث.
لنضع \(B=\lfloor(4N)^{1/4}\rfloor\). عندئذٍ تعمل الحلقة الخارجية \(O(B)=O(N^{1/4})\) مرة.
وبالنسبة إلى قيمة ثابتة من \(b\)، يكون عدد المرشحين في الداخل من رتبة \(O((N/b)^{1/3})\). وعند جمع ذلك على جميع \(b\le B\) نحصل على
$$\sum_{b\le B} O\!\left((N/b)^{1/3}\right)=O\!\left(N^{1/3}\sum_{b\le B} b^{-1/3}\right)=O(N^{1/2}).$$
إذًا فجزء التعداد الحسابي يقارب \(O(N^{1/2})\) فحصًا للمرشحين قبل الاستفادة من عوامل التقليص الناتجة عن شرط التوافق النسبي وفلتر النواة. أما الغربلة اللازمة للنوى الخالية من المكعبات فلا تتجاوز أكبر حجم للمعاملات، أي \(O(N^{1/3})\)، ولذلك تكلفتها \(O(N^{1/3}\log\log N)\) زمنًا و\(O(N^{1/3})\) ذاكرة.
عمليًا هذا أفضل بكثير من تعداد جميع الأزواج \((x,y)\) مباشرةً، كما أن التنفيذ متعدد الخيوط في C++ يخفض الزمن الفعلي من دون أن يغير الرتبة التقاربية.