We seek primitive positive integer triples \((x,y,z)\) satisfying
$$\frac{1}{x^2}+\frac{1}{y^2}=\frac{13}{z^2},\qquad \gcd(x,y,z)=1,$$
together with the bound \(x,y,z\le N\). For every admissible triple we add \(x+y+z\), so the target quantity is
$$S(N)=\sum_{\substack{(x,y,z)\text{ primitive}\\ \frac{1}{x^2}+\frac{1}{y^2}=\frac{13}{z^2}\\ x,y,z\le N}} (x+y+z).$$
A brute-force search over all triples up to \(N=10^{16}\) is impossible, so the solution reorganizes the equation until the search space is generated by primitive Pythagorean parameters and a sharp bound on one auxiliary variable.
Write
$$x=g\,a,\qquad y=g\,b,\qquad \gcd(a,b)=1.$$
Substituting into the reciprocal equation gives
$$z^2(a^2+b^2)=13g^2a^2b^2.$$
Because \(\gcd(a,a^2+b^2)=\gcd(a,b^2)=1\), the factor \(a^2\) must divide \(z^2\), so \(a\mid z\). The same argument shows \(b\mid z\). Since \(\gcd(a,b)=1\), we can write
$$z=abk$$
for some positive integer \(k\). After substitution,
$$k^2(a^2+b^2)=13g^2.$$
The primitive condition \(\gcd(x,y,z)=1\) implies \(\gcd(g,z)=1\), hence \(\gcd(g,k)=1\). Therefore \(k^2\) divides \(13\), so the only possibility is
$$k=1.$$
Thus every primitive solution has the shape
$$x=q\,r,\qquad y=p\,r,\qquad z=pq,$$
where
$$p^2+q^2=13r^2,\qquad \gcd(p,q)=1.$$
The original reciprocal equation has been reduced to finding primitive representations of \(13r^2\) as a sum of two squares.
Any primitive solution of
$$u^2+v^2=r^2$$
is given by Euclid's parametrization
$$u=m^2-n^2,\qquad v=2mn,\qquad r=m^2+n^2,$$
with
$$m>n\ge 0,\qquad m-n\text{ odd},\qquad \gcd(m,n)=1.$$
Allowing \(n=0\) keeps the degenerate base case \(u=1\), \(v=0\), \(r=1\); for larger \(m\) the coprimality test automatically rejects \(n=0\). So the implementation can handle the smallest primitive solution without special-case code.
The identity
$$ (A^2+B^2)(C^2+D^2)=(AC-BD)^2+(AD+BC)^2=(AC+BD)^2+(AD-BC)^2 $$
turns a representation of \(r^2\) into a representation of \(13r^2\). Using \(13=3^2+2^2\) and the pair \((u,v)\) from Step 2 gives two candidate branches:
$$p=\left|3u-2v\right|,\qquad q=\left|2u+3v\right|,$$
or
$$p=\left|3u+2v\right|,\qquad q=\left|2u-3v\right|.$$
In both cases,
$$p^2+q^2=(3^2+2^2)(u^2+v^2)=13r^2.$$
These two branches are the real-coordinate form of multiplying by \(3+2i\) or \(3-2i\) in the Gaussian integers. That viewpoint also explains why the parametrization is complete for primitive solutions.
Once a pair \((p,q)\) with \(p^2+q^2=13r^2\) is known, define
$$x=q\,r,\qquad y=p\,r,\qquad z=pq.$$
Then
$$ \frac{1}{x^2}+\frac{1}{y^2} =\frac{1}{q^2r^2}+\frac{1}{p^2r^2} =\frac{p^2+q^2}{p^2q^2r^2} =\frac{13r^2}{p^2q^2r^2} =\frac{13}{z^2}. $$
The implementations keep only positive pairs and require \(\gcd(p,q)=1\). That already forces the final triple to be primitive. Indeed, if a prime divided both \(r\) and \(p\), then from \(q^2=13r^2-p^2\) it would also divide \(q\), contradicting \(\gcd(p,q)=1\). So \(\gcd(r,p)=\gcd(r,q)=1\), hence \(\gcd(x,y,z)=1\).
To avoid counting the same solution twice with \(x\) and \(y\) swapped, the implementation orders the pair so that \(p\ge q\), which implies \(y\ge x\). It also discards the rare case where the two branches produce the same pair.
The formulas above give
$$x^2+y^2=r^2(p^2+q^2)=13r^4.$$
If \(x\le N\) and \(y\le N\), then
$$x^2+y^2\le 2N^2,$$
so every valid triple must satisfy
$$13r^4\le 2N^2.$$
This yields the sharp search limit
$$r\le r_{\max}=\left\lfloor\left(\frac{2N^2}{13}\right)^{1/4}\right\rfloor.$$
Since \(r=m^2+n^2\), the outer parameter only needs to run up to \(\lfloor\sqrt{r_{\max}}\rfloor\).
For \(N=100\), the bound is
$$13r^4\le 2\cdot 100^2=20000,$$
so \(r_{\max}=6\).
The pair \((m,n)=(1,0)\) gives
$$u=1,\qquad v=0,\qquad r=1.$$
Both branches collapse to the same pair \((p,q)=(3,2)\). Therefore
$$x=2,\qquad y=3,\qquad z=6,$$
and the contribution is \(2+3+6=11\).
The pair \((m,n)=(2,1)\) gives
$$u=3,\qquad v=4,\qquad r=5.$$
The two branches produce \((p,q)=(18,1)\) and \((17,6)\). They map to
$$ (x,y,z)=(5,90,18)\quad \text{and}\quad (30,85,102). $$
Only the first triple satisfies \(x,y,z\le 100\), so its contribution is \(5+90+18=113\).
No other primitive parameter pair has \(r\le 6\). Hence
$$S(100)=11+113=124,$$
which matches the implementation's checkpoint.
The C++, Python, and Java implementations follow the same pipeline. First they compute the largest admissible \(r\) by integer binary search on the inequality \(13r^4\le 2N^2\). That keeps the search exact and avoids floating-point edge cases at the upper boundary.
Next they enumerate all primitive Euclid parameters \((m,n)\) with the parity and coprimality conditions from Step 2. For each such pair, they form the two linear transforms from Step 3, reorder the resulting values so the larger one comes first, discard zero values, suppress a duplicated second branch when both branches coincide, and keep only coprime pairs.
Each surviving pair is converted into \((x,y,z)=(qr,pr,pq)\). If all three coordinates are at most \(N\), the implementation adds \(x+y+z\) to the running total. The arithmetic is done with sufficiently wide integer types so that intermediate products remain exact.
Let
$$R=\left\lfloor\left(\frac{2N^2}{13}\right)^{1/4}\right\rfloor.$$
The parameter search visits primitive pairs \((m,n)\) with \(m^2+n^2\le R\), so the number of tested pairs is \(O(R)\). Because \(R=\Theta(\sqrt{N})\), the total running time is \(O(\sqrt{N})\) with only constant extra memory beyond the integer arithmetic used for safe products and accumulation.
Gesucht sind primitive positive ganzzahlige Tripel \((x,y,z)\) mit
$$\frac{1}{x^2}+\frac{1}{y^2}=\frac{13}{z^2},\qquad \gcd(x,y,z)=1,$$
sowie der Schranke \(x,y,z\le N\). Für jedes zulässige Tripel wird \(x+y+z\) addiert; die Zielgröße ist also
$$S(N)=\sum_{\substack{(x,y,z)\text{ primitiv}\\ \frac{1}{x^2}+\frac{1}{y^2}=\frac{13}{z^2}\\ x,y,z\le N}} (x+y+z).$$
Eine direkte Suche bis \(N=10^{16}\) ist völlig unpraktikabel. Deshalb formt die Lösung die Gleichung so um, dass alle Kandidaten aus primitiven pythagoreischen Parametern und einer scharfen Schranke für eine Hilfsvariable erzeugt werden.
Schreibe
$$x=g\,a,\qquad y=g\,b,\qquad \gcd(a,b)=1.$$
Einsetzen in die Kehrwertgleichung liefert
$$z^2(a^2+b^2)=13g^2a^2b^2.$$
Weil \(\gcd(a,a^2+b^2)=\gcd(a,b^2)=1\), muss \(a^2\) den Ausdruck \(z^2\) teilen; also gilt \(a\mid z\). Genauso folgt \(b\mid z\). Wegen \(\gcd(a,b)=1\) können wir
$$z=abk$$
schreiben. Dann wird die Gleichung zu
$$k^2(a^2+b^2)=13g^2.$$
Aus der Primitivität \(\gcd(x,y,z)=1\) folgt \(\gcd(g,z)=1\) und damit insbesondere \(\gcd(g,k)=1\). Also muss \(k^2\) ein Teiler von \(13\) sein, was nur für
$$k=1$$
möglich ist. Jede primitive Lösung hat somit die Form
$$x=q\,r,\qquad y=p\,r,\qquad z=pq,$$
wobei
$$p^2+q^2=13r^2,\qquad \gcd(p,q)=1.$$
Das ursprüngliche Problem ist damit auf primitive Darstellungen von \(13r^2\) als Summe zweier Quadrate reduziert.
Jede primitive Lösung von
$$u^2+v^2=r^2$$
wird durch Euklids Parametrisierung beschrieben:
$$u=m^2-n^2,\qquad v=2mn,\qquad r=m^2+n^2,$$
mit
$$m>n\ge 0,\qquad m-n\text{ ungerade},\qquad \gcd(m,n)=1.$$
Die Zulassung von \(n=0\) deckt den Grundfall \(u=1\), \(v=0\), \(r=1\) ab; für größere \(m\) verwirft der Teilerfremdheitstest automatisch alle übrigen Fälle mit \(n=0\).
Die Identität
$$ (A^2+B^2)(C^2+D^2)=(AC-BD)^2+(AD+BC)^2=(AC+BD)^2+(AD-BC)^2 $$
macht aus einer Darstellung von \(r^2\) sofort eine Darstellung von \(13r^2\). Mit \(13=3^2+2^2\) und dem Paar \((u,v)\) aus Schritt 2 entstehen zwei Kandidatenzweige:
$$p=\left|3u-2v\right|,\qquad q=\left|2u+3v\right|,$$
oder
$$p=\left|3u+2v\right|,\qquad q=\left|2u-3v\right|.$$
In beiden Fällen gilt
$$p^2+q^2=(3^2+2^2)(u^2+v^2)=13r^2.$$
Diese beiden Zweige entsprechen im Gaußschen Zahlenring der Multiplikation mit \(3+2i\) beziehungsweise \(3-2i\). Dadurch wird zugleich klar, warum die Parametrisierung für primitive Lösungen vollständig ist.
Ist ein Paar \((p,q)\) mit \(p^2+q^2=13r^2\) gefunden, setzen wir
$$x=q\,r,\qquad y=p\,r,\qquad z=pq.$$
Dann folgt sofort
$$ \frac{1}{x^2}+\frac{1}{y^2} =\frac{1}{q^2r^2}+\frac{1}{p^2r^2} =\frac{p^2+q^2}{p^2q^2r^2} =\frac{13r^2}{p^2q^2r^2} =\frac{13}{z^2}. $$
Die Implementierungen behalten nur positive Paare mit \(\gcd(p,q)=1\). Das erzwingt bereits die Primitivität des Endtripels: Teilt eine Primzahl sowohl \(r\) als auch \(p\), dann teilt sie wegen \(q^2=13r^2-p^2\) auch \(q\), im Widerspruch zu \(\gcd(p,q)=1\). Also sind \(r\) und \(p\) sowie \(r\) und \(q\) teilerfremd, daher ist \(\gcd(x,y,z)=1\).
Damit \((x,y,z)\) und \((y,x,z)\) nicht doppelt gezählt werden, wird das Paar so geordnet, dass \(p\ge q\) gilt. Außerdem wird ein seltener doppelter Zweig verworfen, wenn beide Formeln dasselbe Paar erzeugen.
Aus den Formeln folgt
$$x^2+y^2=r^2(p^2+q^2)=13r^4.$$
Wenn \(x\le N\) und \(y\le N\), dann ist
$$x^2+y^2\le 2N^2,$$
also muss jede gültige Lösung die Bedingung
$$13r^4\le 2N^2$$
erfüllen. Damit erhält man die exakte Suchgrenze
$$r\le r_{\max}=\left\lfloor\left(\frac{2N^2}{13}\right)^{1/4}\right\rfloor.$$
Da \(r=m^2+n^2\) ist, genügt es, den äußeren Parameter bis \(\lfloor\sqrt{r_{\max}}\rfloor\) laufen zu lassen.
Für \(N=100\) lautet die Schranke
$$13r^4\le 2\cdot 100^2=20000,$$
also ist \(r_{\max}=6\).
Das Paar \((m,n)=(1,0)\) liefert
$$u=1,\qquad v=0,\qquad r=1.$$
Beide Zweige fallen auf dasselbe Paar \((p,q)=(3,2)\) zusammen. Daraus folgt
$$x=2,\qquad y=3,\qquad z=6,$$
mit Beitrag \(2+3+6=11\).
Das Paar \((m,n)=(2,1)\) ergibt
$$u=3,\qquad v=4,\qquad r=5.$$
Die beiden Zweige erzeugen \((p,q)=(18,1)\) und \((17,6)\). Diese führen zu
$$ (x,y,z)=(5,90,18)\quad \text{und}\quad (30,85,102). $$
Nur das erste Tripel erfüllt \(x,y,z\le 100\), also beträgt sein Beitrag \(5+90+18=113\).
Weitere primitive Parameterpaare mit \(r\le 6\) gibt es nicht. Daher
$$S(100)=11+113=124,$$
genau wie im Prüfwert der Implementierung.
Die C++-, Python- und Java-Implementierungen verwenden dieselbe Ablaufstruktur. Zuerst bestimmen sie per ganzzahliger Binärsuche das größte zulässige \(r\), das noch \(13r^4\le 2N^2\) erfüllt. So bleibt die obere Grenze exakt und es entstehen keine Rundungsprobleme an der Kante.
Danach werden alle primitiven Euklid-Parameter \((m,n)\) mit den Bedingungen aus Schritt 2 durchlaufen. Für jedes Paar werden die beiden linearen Transformationen aus Schritt 3 gebildet, in eine feste Reihenfolge gebracht, bei Null verworfen, im Fall identischer Zweige entdoppelt und nur dann behalten, wenn das resultierende Paar teilerfremd ist.
Jedes verbleibende Paar wird in \((x,y,z)=(qr,pr,pq)\) umgerechnet. Falls alle drei Werte höchstens \(N\) sind, addiert die Implementierung \(x+y+z\) zur Gesamtsumme. Für die Produkte und die Akkumulation werden hinreichend breite Ganzzahltypen benutzt, damit alle Zwischenwerte exakt bleiben.
Setze
$$R=\left\lfloor\left(\frac{2N^2}{13}\right)^{1/4}\right\rfloor.$$
Durchsucht werden primitive Paare \((m,n)\) mit \(m^2+n^2\le R\); deren Anzahl ist \(O(R)\). Da \(R=\Theta(\sqrt{N})\) gilt, beträgt die Gesamtlaufzeit \(O(\sqrt{N})\). Der zusätzliche Speicherbedarf ist abgesehen von den für sichere Ganzzahlarithmetik benötigten Werten konstant.
Aranan şey,
$$\frac{1}{x^2}+\frac{1}{y^2}=\frac{13}{z^2},\qquad \gcd(x,y,z)=1$$
koşullarını sağlayan ilkel pozitif tamsayı üçlüleridir. Ayrıca \(x\le N\), \(y\le N\) ve \(z\le N\) olmalıdır. Her uygun üçlü için \(x+y+z\) toplandığından hedef büyüklük
$$S(N)=\sum_{\substack{(x,y,z)\text{ ilkel}\\ \frac{1}{x^2}+\frac{1}{y^2}=\frac{13}{z^2}\\ x,y,z\le N}} (x+y+z)$$
şeklindedir. \(N=10^{16}\) seviyesinde doğrudan tarama yapılamayacağı için çözüm, denklemi iki kare toplamı parametrizasyonuna indirger.
Şöyle yazalım:
$$x=g\,a,\qquad y=g\,b,\qquad \gcd(a,b)=1.$$
Bunu denkleme koyunca
$$z^2(a^2+b^2)=13g^2a^2b^2$$
elde edilir. \(\gcd(a,a^2+b^2)=\gcd(a,b^2)=1\) olduğundan \(a^2\) mutlaka \(z^2\)'yi bölmelidir; dolayısıyla \(a\mid z\). Aynı mantıkla \(b\mid z\) de çıkar. \(\gcd(a,b)=1\) olduğu için
$$z=abk$$
yazabiliriz. Yerine koyarsak
$$k^2(a^2+b^2)=13g^2$$
olur. İlkel olma koşulu \(\gcd(x,y,z)=1\), özellikle \(\gcd(g,z)=1\) anlamına gelir; buradan \(\gcd(g,k)=1\) gelir. O halde \(k^2\), \(13\)'ü bölmelidir ve tek olasılık
$$k=1$$
olur. Böylece her ilkel çözüm şu biçime iner:
$$x=q\,r,\qquad y=p\,r,\qquad z=pq,$$
ve burada
$$p^2+q^2=13r^2,\qquad \gcd(p,q)=1.$$
Yani asıl problem, \(13r^2\)'yi iki kare toplamı olarak ilkel biçimde yazmaktır.
Şimdi
$$u^2+v^2=r^2$$
eşitliğinin her ilkel çözümü, Euclid parametrizasyonuyla
$$u=m^2-n^2,\qquad v=2mn,\qquad r=m^2+n^2$$
şeklinde yazılır. Koşullar
$$m>n\ge 0,\qquad m-n\text{ tek},\qquad \gcd(m,n)=1$$
olmalıdır. \(n=0\) seçeneği, en küçük temel durum olan \(u=1\), \(v=0\), \(r=1\)'i otomatik olarak kapsar; daha büyük \(m\) değerlerinde ise aralarında asal olma filtresi zaten bunları eler.
Kullanılan temel özdeşlik
$$ (A^2+B^2)(C^2+D^2)=(AC-BD)^2+(AD+BC)^2=(AC+BD)^2+(AD-BC)^2 $$
özdeşliği, \(r^2\)'nin bir temsilini doğrudan \(13r^2\)'nin temsiline dönüştürür. \(13=3^2+2^2\) ve Adım 2'deki \((u,v)\) kullanılırsa iki aday dal elde edilir:
$$p=\left|3u-2v\right|,\qquad q=\left|2u+3v\right|,$$
veya
$$p=\left|3u+2v\right|,\qquad q=\left|2u-3v\right|.$$
Her iki durumda da
$$p^2+q^2=(3^2+2^2)(u^2+v^2)=13r^2$$
sağlanır. Bu iki formül, Gauss tamsayılarında \(3+2i\) ve \(3-2i\) ile çarpmanın gerçek bileşenlere açılmış halidir; bu bakış aynı zamanda parametrizasyonun neden tam olduğunu açıklar.
\(p^2+q^2=13r^2\) sağlayan bir çift bulunduğunda
$$x=q\,r,\qquad y=p\,r,\qquad z=pq$$
tanımlanır. O zaman
$$ \frac{1}{x^2}+\frac{1}{y^2} =\frac{1}{q^2r^2}+\frac{1}{p^2r^2} =\frac{p^2+q^2}{p^2q^2r^2} =\frac{13r^2}{p^2q^2r^2} =\frac{13}{z^2} $$
olur. Uygulama yalnızca pozitif ve \(\gcd(p,q)=1\) olan çiftleri tutar. Bu zaten son üçlünün ilkel olmasını garanti eder: Bir asal hem \(r\)'yi hem \(p\)'yi bölseydi, \(q^2=13r^2-p^2\) nedeniyle \(q\)'yu da bölerdi; bu da \(\gcd(p,q)=1\) ile çelişirdi.
Ayrıca aynı çözümü \(x\) ve \(y\) yer değiştirmiş halde iki kez saymamak için çift \(p\ge q\) olacak şekilde sıralanır. Nadir durumlarda iki aday dal aynı çifti üretirse ikinci kopya atılır.
Yukarıdaki tanımlardan
$$x^2+y^2=r^2(p^2+q^2)=13r^4$$
elde edilir. Eğer \(x\le N\) ve \(y\le N\) ise
$$x^2+y^2\le 2N^2$$
olduğundan her geçerli çözüm için
$$13r^4\le 2N^2$$
zorunludur. Dolayısıyla keskin sınır
$$r\le r_{\max}=\left\lfloor\left(\frac{2N^2}{13}\right)^{1/4}\right\rfloor$$
olur. \(r=m^2+n^2\) olduğu için dış döngüde \(m\) ancak \(\lfloor\sqrt{r_{\max}}\rfloor\) değerine kadar gitmelidir.
\(N=100\) için sınır
$$13r^4\le 2\cdot 100^2=20000$$
olduğundan \(r_{\max}=6\) bulunur.
\((m,n)=(1,0)\) çifti
$$u=1,\qquad v=0,\qquad r=1$$
verir. İki aday dal aynı \((p,q)=(3,2)\) çiftine düşer. Böylece
$$x=2,\qquad y=3,\qquad z=6$$
ve katkı \(2+3+6=11\) olur.
\((m,n)=(2,1)\) çifti için
$$u=3,\qquad v=4,\qquad r=5$$
elde edilir. Dallar \((p,q)=(18,1)\) ve \((17,6)\) üretir. Bunlar
$$ (x,y,z)=(5,90,18)\quad \text{ve}\quad (30,85,102) $$
üçlülerine gider. Bunlardan sadece ilki \(x,y,z\le 100\) koşulunu sağlar; katkısı \(5+90+18=113\)'tür.
\(r\le 6\) için başka ilkel parametre çifti yoktur. Sonuç olarak
$$S(100)=11+113=124$$
bulunur; bu da uygulamadaki kontrol değeriyle aynıdır.
C++, Python ve Java uygulamaları aynı akışı izler. Önce \(13r^4\le 2N^2\) eşitsizliğini tam sayı ikili aramayla çözerek en büyük uygun \(r\) bulunur. Böylece üst sınır kesin olur ve kayan nokta sınır hatalarından kaçınılır.
Ardından Adım 2'deki parity ve aralarında asal olma koşullarını sağlayan tüm Euclid parametreleri \((m,n)\) dolaşılır. Her çift için Adım 3'teki iki lineer dönüşüm hesaplanır, büyük olan başa alınır, sıfır içeren durumlar ve yinelenen ikinci dal atılır, sonra da yalnızca aralarında asal çiftler tutulur.
Kalan her aday \((x,y,z)=(qr,pr,pq)\) biçimine çevrilir. Eğer üç bileşenin hepsi \(N\)'yi aşmıyorsa toplam değere \(x+y+z\) eklenir. Ara çarpımların ve toplamın taşmaması için yeterince geniş tamsayı aritmetiği kullanılır.
Tanım gereği
$$R=\left\lfloor\left(\frac{2N^2}{13}\right)^{1/4}\right\rfloor$$
olsun. Aranan bölge \(m^2+n^2\le R\) koşulunu sağlayan ilkel çiftlerdir; ziyaret edilen çift sayısı \(O(R)\) mertebesindedir. \(R=\Theta(\sqrt{N})\) olduğundan toplam çalışma süresi \(O(\sqrt{N})\), ek bellek kullanımı ise güvenli tamsayı aritmetiği dışında \(O(1)\)'dir.
Buscamos ternas primitivas de enteros positivos \((x,y,z)\) que satisfacen
$$\frac{1}{x^2}+\frac{1}{y^2}=\frac{13}{z^2},\qquad \gcd(x,y,z)=1,$$
con la restricción adicional \(x,y,z\le N\). Para cada terna válida se suma \(x+y+z\), así que la cantidad objetivo es
$$S(N)=\sum_{\substack{(x,y,z)\text{ primitiva}\\ \frac{1}{x^2}+\frac{1}{y^2}=\frac{13}{z^2}\\ x,y,z\le N}} (x+y+z).$$
Recorrer directamente todas las ternas hasta \(N=10^{16}\) es inviable. La solución transforma la ecuación recíproca en una parametrización de suma de dos cuadrados con un límite de búsqueda muy estrecho.
Escribimos
$$x=g\,a,\qquad y=g\,b,\qquad \gcd(a,b)=1.$$
Al sustituir en la ecuación se obtiene
$$z^2(a^2+b^2)=13g^2a^2b^2.$$
Como \(\gcd(a,a^2+b^2)=\gcd(a,b^2)=1\), se deduce que \(a^2\) divide a \(z^2\), luego \(a\mid z\). Lo mismo vale para \(b\). Dado que \(\gcd(a,b)=1\), podemos poner
$$z=abk.$$
Entonces la ecuación pasa a ser
$$k^2(a^2+b^2)=13g^2.$$
La condición de primitividad \(\gcd(x,y,z)=1\) implica \(\gcd(g,z)=1\), así que también \(\gcd(g,k)=1\). Por tanto \(k^2\) debe dividir a \(13\), y la única posibilidad es
$$k=1.$$
Toda solución primitiva queda entonces en la forma
$$x=q\,r,\qquad y=p\,r,\qquad z=pq,$$
con
$$p^2+q^2=13r^2,\qquad \gcd(p,q)=1.$$
El problema original queda reducido a describir representaciones primitivas de \(13r^2\) como suma de dos cuadrados.
Cualquier solución primitiva de
$$u^2+v^2=r^2$$
se escribe con la parametrización de Euclides:
$$u=m^2-n^2,\qquad v=2mn,\qquad r=m^2+n^2,$$
bajo las condiciones
$$m>n\ge 0,\qquad m-n\text{ impar},\qquad \gcd(m,n)=1.$$
Permitir \(n=0\) conserva el caso base degenerado \(u=1\), \(v=0\), \(r=1\); para valores mayores de \(m\), la prueba de coprimalidad elimina automáticamente los demás casos con \(n=0\).
La identidad
$$ (A^2+B^2)(C^2+D^2)=(AC-BD)^2+(AD+BC)^2=(AC+BD)^2+(AD-BC)^2 $$
convierte una representación de \(r^2\) en una representación de \(13r^2\). Usando \(13=3^2+2^2\) y el par \((u,v)\) del paso anterior aparecen dos ramas candidatas:
$$p=\left|3u-2v\right|,\qquad q=\left|2u+3v\right|,$$
o bien
$$p=\left|3u+2v\right|,\qquad q=\left|2u-3v\right|.$$
En ambos casos,
$$p^2+q^2=(3^2+2^2)(u^2+v^2)=13r^2.$$
Estas dos fórmulas son la versión en coordenadas reales de multiplicar por \(3+2i\) o por \(3-2i\) en los enteros gaussianos. Esa interpretación también explica por qué la parametrización cubre todas las soluciones primitivas.
Una vez conocido un par \((p,q)\) con \(p^2+q^2=13r^2\), definimos
$$x=q\,r,\qquad y=p\,r,\qquad z=pq.$$
Entonces
$$ \frac{1}{x^2}+\frac{1}{y^2} =\frac{1}{q^2r^2}+\frac{1}{p^2r^2} =\frac{p^2+q^2}{p^2q^2r^2} =\frac{13r^2}{p^2q^2r^2} =\frac{13}{z^2}. $$
Las implementaciones solo aceptan pares positivos y con \(\gcd(p,q)=1\). Eso ya obliga a que la terna final sea primitiva: si un primo dividiera a la vez \(r\) y \(p\), entonces de \(q^2=13r^2-p^2\) también dividiría a \(q\), contradiciendo \(\gcd(p,q)=1\).
Para no contar dos veces la misma solución con \(x\) e \(y\) intercambiados, el par se ordena de modo que \(p\ge q\). Si las dos ramas producen exactamente el mismo par, la segunda copia se descarta.
De las fórmulas anteriores se deduce
$$x^2+y^2=r^2(p^2+q^2)=13r^4.$$
Si \(x\le N\) y \(y\le N\), entonces
$$x^2+y^2\le 2N^2,$$
y toda solución válida debe cumplir
$$13r^4\le 2N^2.$$
Por tanto, la frontera exacta es
$$r\le r_{\max}=\left\lfloor\left(\frac{2N^2}{13}\right)^{1/4}\right\rfloor.$$
Como \(r=m^2+n^2\), basta con hacer correr el parámetro exterior hasta \(\lfloor\sqrt{r_{\max}}\rfloor\).
Cuando \(N=100\), la cota es
$$13r^4\le 2\cdot 100^2=20000,$$
de modo que \(r_{\max}=6\).
El par \((m,n)=(1,0)\) produce
$$u=1,\qquad v=0,\qquad r=1.$$
Las dos ramas colapsan en el mismo par \((p,q)=(3,2)\). Entonces
$$x=2,\qquad y=3,\qquad z=6,$$
y la contribución es \(2+3+6=11\).
El par \((m,n)=(2,1)\) produce
$$u=3,\qquad v=4,\qquad r=5.$$
Las dos ramas dan \((p,q)=(18,1)\) y \((17,6)\). Eso lleva a
$$ (x,y,z)=(5,90,18)\quad \text{y}\quad (30,85,102). $$
Solo la primera terna cumple \(x,y,z\le 100\), así que aporta \(5+90+18=113\).
No hay más pares primitivos con \(r\le 6\). Por eso
$$S(100)=11+113=124,$$
exactamente el valor de comprobación usado por la implementación.
Las implementaciones en C++, Python y Java siguen la misma secuencia. Primero calculan el mayor \(r\) admisible mediante una búsqueda binaria entera sobre la desigualdad \(13r^4\le 2N^2\). Así el límite superior queda exacto y no depende de redondeos de punto flotante.
Después enumeran todos los pares primitivos \((m,n)\) de Euclides que cumplen las condiciones de paridad y coprimalidad del Paso 2. Para cada par forman las dos transformaciones lineales del Paso 3, ordenan el resultado para que el mayor valor quede primero, descartan ceros, eliminan una segunda rama duplicada cuando coincide con la primera y conservan solo los pares coprimos.
Cada par superviviente se transforma en \((x,y,z)=(qr,pr,pq)\). Si las tres coordenadas no superan \(N\), la implementación suma \(x+y+z\) al acumulado. Las multiplicaciones intermedias y la acumulación se hacen con tipos enteros suficientemente amplios para mantener la exactitud.
Sea
$$R=\left\lfloor\left(\frac{2N^2}{13}\right)^{1/4}\right\rfloor.$$
La búsqueda recorre pares primitivos \((m,n)\) con \(m^2+n^2\le R\), por lo que el número de candidatos visitados es \(O(R)\). Como \(R=\Theta(\sqrt{N})\), el tiempo total es \(O(\sqrt{N})\). El uso de memoria adicional es \(O(1)\), aparte de la aritmética entera necesaria para que los productos sean seguros.
题目要求统计所有本原正整数三元组 \((x,y,z)\)，满足
$$\frac{1}{x^2}+\frac{1}{y^2}=\frac{13}{z^2},\qquad \gcd(x,y,z)=1,$$
并且 \(x\le N\)、\(y\le N\)、\(z\le N\)。对每个合法三元组都要累加 \(x+y+z\)，因此目标函数是
$$S(N)=\sum_{\substack{(x,y,z)\text{ 为本原解}\\ \frac{1}{x^2}+\frac{1}{y^2}=\frac{13}{z^2}\\ x,y,z\le N}} (x+y+z).$$
如果直接在 \([1,N]\) 内枚举三元组，那么在 \(N=10^{16}\) 时完全不可行。实现采用的思路是先把倒数平方方程化成一个“平方和”问题，再用本原勾股参数生成全部候选，并用一个非常紧的上界裁掉绝大多数搜索空间。
设
$$x=g\,a,\qquad y=g\,b,\qquad \gcd(a,b)=1.$$
代回原式可得
$$z^2(a^2+b^2)=13g^2a^2b^2.$$
由于 \(\gcd(a,a^2+b^2)=\gcd(a,b^2)=1\)，所以 \(a^2\) 必须整除 \(z^2\)，于是 \(a\mid z\)。同理也有 \(b\mid z\)。又因为 \(a\) 与 \(b\) 互素，所以可以写成
$$z=abk.$$
再代回去，就得到
$$k^2(a^2+b^2)=13g^2.$$
本原条件 \(\gcd(x,y,z)=1\) 意味着 \(\gcd(g,z)=1\)，因此特别有 \(\gcd(g,k)=1\)。于是 \(k^2\) 必须是 \(13\) 的平方因子，只可能是
$$k=1.$$
所以每个本原解都可以改写为
$$x=q\,r,\qquad y=p\,r,\qquad z=pq,$$
其中
$$p^2+q^2=13r^2,\qquad \gcd(p,q)=1.$$
这一步非常关键，因为它说明原问题等价于寻找 \(13r^2\) 的本原平方和表示。
所有本原解
$$u^2+v^2=r^2$$
都可以由 Euclid 参数化给出：
$$u=m^2-n^2,\qquad v=2mn,\qquad r=m^2+n^2,$$
条件是
$$m>n\ge 0,\qquad m-n\text{ 为奇数},\qquad \gcd(m,n)=1.$$
实现里允许 \(n=0\)，这样最小的退化情形 \(u=1\)、\(v=0\)、\(r=1\) 会自然出现；而当 \(m>1\) 时，互素筛选会自动排掉其余 \(n=0\) 的情况。因此不需要额外写一个“最小解特判”。
平方和有经典恒等式
$$ (A^2+B^2)(C^2+D^2)=(AC-BD)^2+(AD+BC)^2=(AC+BD)^2+(AD-BC)^2 $$
把 \(13=3^2+2^2\) 与步骤 2 中的 \(u^2+v^2=r^2\) 合成，就会得到 \(13r^2\) 的两组候选表示：
$$p=\left|3u-2v\right|,\qquad q=\left|2u+3v\right|,$$
或者
$$p=\left|3u+2v\right|,\qquad q=\left|2u-3v\right|.$$
两种写法都满足
$$p^2+q^2=(3^2+2^2)(u^2+v^2)=13r^2.$$
如果从高层一点的角度看，这正是高斯整数里乘上 \(3+2i\) 或 \(3-2i\) 的实部与虚部。也正因为 \(13\) 在高斯整数中这样分解，所以这两条分支不仅能构造解，而且能覆盖全部本原情形。
只要拿到一个满足 \(p^2+q^2=13r^2\) 的正整数对，就定义
$$x=q\,r,\qquad y=p\,r,\qquad z=pq.$$
于是
$$ \frac{1}{x^2}+\frac{1}{y^2} =\frac{1}{q^2r^2}+\frac{1}{p^2r^2} =\frac{p^2+q^2}{p^2q^2r^2} =\frac{13r^2}{p^2q^2r^2} =\frac{13}{z^2}. $$
实现会强制 \(p,q>0\)，并要求 \(\gcd(p,q)=1\)。这一条已经足以保证最后的 \((x,y,z)\) 也是本原的。原因是：如果某个素数同时整除 \(r\) 和 \(p\)，那么由
$$q^2=13r^2-p^2$$
可知它也会整除 \(q\)，这与 \(\gcd(p,q)=1\) 矛盾。因此 \(r\) 与 \(p\)、\(r\) 与 \(q\) 都互素，最终 \(\gcd(x,y,z)=1\)。
为了避免把 \((x,y,z)\) 和交换了 \(x,y\) 的同一个解重复计算，程序会把候选对整理成 \(p\ge q\) 的顺序。如果两条分支恰好产生同一对值，那么第二份会被直接跳过。
根据上面的定义，
$$x^2+y^2=r^2(p^2+q^2)=13r^4.$$
而对所有合法解都有 \(x\le N\)、\(y\le N\)，因此
$$x^2+y^2\le 2N^2,$$
于是必然满足
$$13r^4\le 2N^2.$$
这就得到一个非常实用的精确上界：
$$r\le r_{\max}=\left\lfloor\left(\frac{2N^2}{13}\right)^{1/4}\right\rfloor.$$
又因为 \(r=m^2+n^2\)，所以外层参数 \(m\) 只需要枚举到 \(\lfloor\sqrt{r_{\max}}\rfloor\) 为止。整个搜索区间因此从三维枚举压缩成了一个很薄的二维区域。
当 \(N=100\) 时，有
$$13r^4\le 2\cdot 100^2=20000,$$
所以 \(r_{\max}=6\)。
先看 \((m,n)=(1,0)\)。这时
$$u=1,\qquad v=0,\qquad r=1.$$
两条分支都给出同一个候选 \((p,q)=(3,2)\)。因此
$$x=2,\qquad y=3,\qquad z=6,$$
贡献是 \(2+3+6=11\)。
再看 \((m,n)=(2,1)\)。这时
$$u=3,\qquad v=4,\qquad r=5.$$
两条分支分别得到 \((p,q)=(18,1)\) 与 \((17,6)\)。映射回原变量后是
\((x,y,z)=(5,90,18)\) 和 \((30,85,102)\)。
第二个三元组因为 \(z=102>100\) 被过滤掉，第一个三元组保留，贡献 \(5+90+18=113\)。
在 \(r\le 6\) 的范围里没有别的本原参数对，所以
$$S(100)=11+113=124,$$
与实现中的检验值完全一致。
C++、Python 和 Java 三个实现的流程是一样的。首先通过整数二分搜索求出满足 \(13r^4\le 2N^2\) 的最大 \(r\)。这样做的好处是边界精确，不会因为浮点近似而漏解或多算。
然后程序遍历所有满足步骤 2 条件的 Euclid 参数 \((m,n)\)。对每一对参数，先算出 \(u\)、\(v\)、\(r\)，再生成步骤 3 中的两条线性变换结果。之后把较大的数放在前面，去掉含零候选，去掉与第一条完全重复的第二条分支，并且只保留互素的 \((p,q)\)。
每个保留下来的候选都按 \((x,y,z)=(qr,pr,pq)\) 转回原问题。如果三项都不超过 \(N\)，就把 \(x+y+z\) 加入总和。实现还使用了足够宽的整数算术，以保证中间乘法和最终累加都保持精确。
记
$$R=\left\lfloor\left(\frac{2N^2}{13}\right)^{1/4}\right\rfloor.$$
搜索只会访问满足 \(m^2+n^2\le R\) 的本原参数对，因此候选数量是 \(O(R)\)。又因为 \(R=\Theta(\sqrt{N})\)，所以总时间复杂度是 \(O(\sqrt{N})\)。除去保证整数乘法安全所需的数值对象之外，额外空间复杂度为 \(O(1)\)。
Нужно найти все примитивные положительные тройки \((x,y,z)\), удовлетворяющие условию
$$\frac{1}{x^2}+\frac{1}{y^2}=\frac{13}{z^2},\qquad \gcd(x,y,z)=1,$$
при ограничении \(x,y,z\le N\). Для каждой допустимой тройки суммируется \(x+y+z\), то есть вычисляется
$$S(N)=\sum_{\substack{(x,y,z)\text{ примитивно}\\ \frac{1}{x^2}+\frac{1}{y^2}=\frac{13}{z^2}\\ x,y,z\le N}} (x+y+z).$$
Прямой перебор троек до \(N=10^{16}\) невозможен. Поэтому решение сначала преобразует уравнение с обратными квадратами в задачу о представлении числа в виде суммы двух квадратов, а затем использует параметризацию примитивных пифагоровых троек.
Пусть
$$x=g\,a,\qquad y=g\,b,\qquad \gcd(a,b)=1.$$
Тогда исходное уравнение превращается в
$$z^2(a^2+b^2)=13g^2a^2b^2.$$
Поскольку \(\gcd(a,a^2+b^2)=\gcd(a,b^2)=1\), число \(a^2\) обязано делить \(z^2\), значит \(a\mid z\). Аналогично \(b\mid z\). Из взаимной простоты \(a\) и \(b\) следует, что
$$z=abk$$
для некоторого целого \(k>0\). После подстановки получаем
$$k^2(a^2+b^2)=13g^2.$$
Примитивность \(\gcd(x,y,z)=1\) влечет \(\gcd(g,z)=1\), а значит и \(\gcd(g,k)=1\). Следовательно, \(k^2\) должно делить \(13\), и единственный возможный вариант:
$$k=1.$$
Тем самым каждая примитивная тройка имеет вид
$$x=q\,r,\qquad y=p\,r,\qquad z=pq,$$
где
$$p^2+q^2=13r^2,\qquad \gcd(p,q)=1.$$
Итак, исходная задача сведена к описанию примитивных представлений \(13r^2\) в виде суммы двух квадратов.
Любое примитивное решение уравнения
$$u^2+v^2=r^2$$
задается формулами Евклида:
$$u=m^2-n^2,\qquad v=2mn,\qquad r=m^2+n^2,$$
при условиях
$$m>n\ge 0,\qquad m-n\text{ нечетно},\qquad \gcd(m,n)=1.$$
Разрешение \(n=0\) удобно тем, что случай \(u=1\), \(v=0\), \(r=1\) появляется автоматически. Для больших \(m\) остальные варианты с \(n=0\) сами отбрасываются проверкой взаимной простоты.
Используется тождество
$$ (A^2+B^2)(C^2+D^2)=(AC-BD)^2+(AD+BC)^2=(AC+BD)^2+(AD-BC)^2. $$
Если взять \(13=3^2+2^2\) и представление \(r^2=u^2+v^2\), получаются две ветви кандидатов:
$$p=\left|3u-2v\right|,\qquad q=\left|2u+3v\right|,$$
или
$$p=\left|3u+2v\right|,\qquad q=\left|2u-3v\right|.$$
В обоих случаях выполняется
$$p^2+q^2=(3^2+2^2)(u^2+v^2)=13r^2.$$
С точки зрения гауссовых целых это просто умножение на \(3+2i\) или на \(3-2i\). Именно эта факторизация объясняет полноту параметризации для примитивных решений.
Как только найдено положительное решение \((p,q)\) для \(p^2+q^2=13r^2\), задаем
$$x=q\,r,\qquad y=p\,r,\qquad z=pq.$$
Тогда
$$ \frac{1}{x^2}+\frac{1}{y^2} =\frac{1}{q^2r^2}+\frac{1}{p^2r^2} =\frac{p^2+q^2}{p^2q^2r^2} =\frac{13r^2}{p^2q^2r^2} =\frac{13}{z^2}. $$
В реализациях сохраняются только положительные взаимно простые пары \((p,q)\). Этого уже достаточно для примитивности итоговой тройки. Если бы простой делитель делил и \(r\), и \(p\), то из равенства \(q^2=13r^2-p^2\) следовало бы, что он делит и \(q\), что невозможно при \(\gcd(p,q)=1\).
Чтобы не считать одну и ту же тройку дважды с переставленными \(x\) и \(y\), пара упорядочивается так, чтобы \(p\ge q\). Если две ветви вдруг дают один и тот же результат, вторая копия просто пропускается.
Из построения следует
$$x^2+y^2=r^2(p^2+q^2)=13r^4.$$
При \(x\le N\) и \(y\le N\) имеем
$$x^2+y^2\le 2N^2,$$
поэтому для любой допустимой тройки обязательно
$$13r^4\le 2N^2.$$
Значит, достаточно рассматривать только
$$r\le r_{\max}=\left\lfloor\left(\frac{2N^2}{13}\right)^{1/4}\right\rfloor.$$
Так как \(r=m^2+n^2\), внешний параметр \(m\) нужно перебирать лишь до \(\lfloor\sqrt{r_{\max}}\rfloor\).
Для \(N=100\) имеем
$$13r^4\le 2\cdot 100^2=20000,$$
откуда \(r_{\max}=6\).
Пара \((m,n)=(1,0)\) дает
$$u=1,\qquad v=0,\qquad r=1.$$
Обе ветви совпадают и дают \((p,q)=(3,2)\). Следовательно,
$$x=2,\qquad y=3,\qquad z=6,$$
а вклад равен \(2+3+6=11\).
Пара \((m,n)=(2,1)\) дает
$$u=3,\qquad v=4,\qquad r=5.$$
Из двух ветвей получаются \((p,q)=(18,1)\) и \((17,6)\). После перехода к исходным переменным это
$$ (x,y,z)=(5,90,18)\quad \text{и}\quad (30,85,102). $$
Вторая тройка отбрасывается, потому что \(z=102>100\). Первая вносит \(5+90+18=113\).
Других примитивных параметров с \(r\le 6\) нет, поэтому
$$S(100)=11+113=124,$$
что совпадает с контрольным значением реализации.
Во всех трех реализациях используется один и тот же конвейер. Сначала целочисленным бинарным поиском находится наибольшее \(r\), для которого верно \(13r^4\le 2N^2\). Это дает точную границу без риска ошибок округления.
Затем перебираются все пары \((m,n)\), удовлетворяющие условиям шага 2. Для каждой пары строятся две линейные формы из шага 3, результат упорядочивается по убыванию, нулевые варианты отбрасываются, совпавшая вторая ветвь подавляется, и остаются только взаимно простые пары.
Каждая оставшаяся пара преобразуется в \((x,y,z)=(qr,pr,pq)\). Если все три значения не превосходят \(N\), к сумме добавляется \(x+y+z\). Для промежуточных произведений и накопления используются достаточно широкие целочисленные типы, чтобы арифметика оставалась точной.
Пусть
$$R=\left\lfloor\left(\frac{2N^2}{13}\right)^{1/4}\right\rfloor.$$
Перебор посещает только примитивные пары \((m,n)\) с \(m^2+n^2\le R\), так что число кандидатов равно \(O(R)\). Поскольку \(R=\Theta(\sqrt{N})\), полная временная сложность составляет \(O(\sqrt{N})\). Дополнительная память постоянна, если не считать объектов, необходимых для безопасных целочисленных произведений.
نبحث عن جميع الثلاثيات الصحيحة الموجبة الأولية \((x,y,z)\) التي تحقق
$$\frac{1}{x^2}+\frac{1}{y^2}=\frac{13}{z^2},\qquad \gcd(x,y,z)=1,$$
مع الشرط \(x,y,z\le N\). ولكل ثلاثية صالحة نجمع \(x+y+z\)، ولذلك تكون الكمية المطلوبة هي
$$S(N)=\sum (x+y+z).$$
ويمتد هذا الجمع على جميع الثلاثيات الأولية التي تحقق \(\frac{1}{x^2}+\frac{1}{y^2}=\frac{13}{z^2}\) مع الشرط \(x,y,z\le N\).
الفحص المباشر لكل الثلاثيات حتى \(N=10^{16}\) غير ممكن عمليا. لذلك تعيد الفكرة صياغة المعادلة ذات المقلوبات إلى مسألة من نوع مجموع مربعين، ثم تولد كل الحلول من بارامترات فيثاغورية أولية مع حد بحث دقيق جدا.
نكتب
$$x=g\,a,\qquad y=g\,b,\qquad \gcd(a,b)=1.$$
وبالتعويض في المعادلة الأصلية نحصل على
$$z^2(a^2+b^2)=13g^2a^2b^2.$$
ولأن \(\gcd(a,a^2+b^2)=\gcd(a,b^2)=1\)، فلا بد أن يقسم \(a^2\) العدد \(z^2\)، ومن ثم \(a\mid z\). وبالطريقة نفسها نحصل على \(b\mid z\). وبما أن \(a\) و\(b\) متباينان أوليا فيمكن كتابة
$$z=abk.$$
بعد ذلك تصبح المعادلة
$$k^2(a^2+b^2)=13g^2.$$
ومن شرط الأولية \(\gcd(x,y,z)=1\) نستنتج \(\gcd(g,z)=1\)، وبالتالي \(\gcd(g,k)=1\). لذا يجب أن يكون \(k^2\) قاسما للعدد \(13\)، وهذا لا يترك إلا الاحتمال
$$k=1.$$
إذن كل حل أولي يمكن كتابته على الصورة
$$x=q\,r,\qquad y=p\,r,\qquad z=pq,$$
حيث
$$p^2+q^2=13r^2,\qquad \gcd(p,q)=1.$$
وهكذا تتحول المسألة من معادلة مقلوبة إلى تمثيل أولي للعدد \(13r^2\) كمجموع مربعين.
كل حل أولي للمعادلة
$$u^2+v^2=r^2$$
يعطى بصيغة إقليدس:
$$u=m^2-n^2,\qquad v=2mn,\qquad r=m^2+n^2,$$
مع الشروط
$$m>n\ge 0,\qquad \gcd(m,n)=1.$$
ويجب أيضا أن يكون \(m-n\) عددا فرديا.
السماح بالحالة \(n=0\) يجعل الحل الأساسي \(u=1\)، \(v=0\)، \(r=1\) يظهر تلقائيا. أما القيم الأخرى ذات \(n=0\) فتستبعدها تلقائيا قاعدة التباين الأولي.
نستخدم الهوية
$$ (A^2+B^2)(C^2+D^2)=(AC-BD)^2+(AD+BC)^2=(AC+BD)^2+(AD-BC)^2. $$
وبما أن \(13=3^2+2^2\)، فإن دمجها مع \(u^2+v^2=r^2\) يعطي فرعين من المرشحين:
$$p=\left|3u-2v\right|,\qquad q=\left|2u+3v\right|,$$
أو
$$p=\left|3u+2v\right|,\qquad q=\left|2u-3v\right|.$$
وفي الحالتين نحصل على
$$p^2+q^2=(3^2+2^2)(u^2+v^2)=13r^2.$$
ومن منظور الأعداد الغاوسية، هذان الفرعان هما فقط ناتج الضرب في \(3+2i\) أو في \(3-2i\). وهذا يفسر أيضا لماذا تكون هذه الصياغة كاملة بالنسبة للحلول الأولية.
بعد الحصول على زوج موجب \((p,q)\) يحقق \(p^2+q^2=13r^2\)، نعرف
$$x=q\,r,\qquad y=p\,r,\qquad z=pq.$$
وعندئذ
$$ \frac{1}{x^2}+\frac{1}{y^2} =\frac{1}{q^2r^2}+\frac{1}{p^2r^2} =\frac{p^2+q^2}{p^2q^2r^2} =\frac{13r^2}{p^2q^2r^2} =\frac{13}{z^2}. $$
التنفيذ يحتفظ فقط بالأزواج الموجبة ذات \(\gcd(p,q)=1\). وهذا يكفي لضمان أن الثلاثية النهائية أولية. فلو كان هناك عدد أولي يقسم كلا من \(r\) و\(p\)، فالعلاقة
$$q^2=13r^2-p^2$$
تجبره على قسمة \(q\) أيضا، وهو ما يناقض \(\gcd(p,q)=1\). لذلك يكون \(r\) متباينا أوليا مع \(p\) ومع \(q\)، وبالتالي \(\gcd(x,y,z)=1\).
ولكي لا نحسب الحل نفسه مرتين بعد تبديل \(x\) و\(y\)، يرتب التنفيذ الزوج بحيث يكون \(p\ge q\). وإذا أعطى الفرعان الزوج نفسه بالضبط، فإن النسخة الثانية تهمل.
من التعريفات السابقة نحصل على
$$x^2+y^2=r^2(p^2+q^2)=13r^4.$$
ومادام \(x\le N\) و\(y\le N\)، فلدينا
$$x^2+y^2\le 2N^2,$$
ومن ثم لا بد لأي حل صالح أن يحقق
$$13r^4\le 2N^2.$$
وعليه يكفي البحث حتى
$$r\le r_{\max}=\left\lfloor\left(\frac{2N^2}{13}\right)^{1/4}\right\rfloor.$$
وبما أن \(r=m^2+n^2\)، فإن البارامتر الخارجي \(m\) لا يحتاج إلى تجاوز \(\lfloor\sqrt{r_{\max}}\rfloor\).
عندما يكون \(N=100\)، نحصل على
$$13r^4\le 2\cdot 100^2=20000,$$
ومن ثم \(r_{\max}=6\).
الزوج \((m,n)=(1,0)\) يعطي
$$u=1,\qquad v=0,\qquad r=1.$$
والفرعان هنا ينهاران إلى الزوج نفسه \((p,q)=(3,2)\). وبالتالي
$$x=2,\qquad y=3,\qquad z=6,$$
ومساهمته \(2+3+6=11\).
أما \((m,n)=(2,1)\) فيعطي
$$u=3,\qquad v=4,\qquad r=5.$$
وينتج عنه الزوجان \((p,q)=(18,1)\) و\((17,6)\). وبعد الرجوع إلى المتغيرات الأصلية نحصل على
\((x,y,z)=(5,90,18)\) و\((30,85,102)\).
الثلاثية الثانية تستبعد لأن \(z=102>100\)، بينما تبقى الأولى وتضيف \(5+90+18=113\).
ولا توجد أزواج أولية أخرى مع \(r\le 6\)، ولذلك
$$S(100)=11+113=124,$$
وهو نفس مقدار الفحص المستخدم في التنفيذ.
تنفذ النسخ المكتوبة بلغة C++ وPython وJava الخطوات نفسها. في البداية تحسب أكبر قيمة ممكنة لـ \(r\) بواسطة بحث ثنائي صحيح على المتراجحة \(13r^4\le 2N^2\). بهذه الطريقة يكون الحد الأعلى دقيقا تماما ولا يعتمد على تقريب الأعداد العائمة.
بعد ذلك تمر على جميع أزواج إقليدس \((m,n)\) التي تحقق شرط الفردية والتباين الأولي في الخطوة الثانية. ولكل زوج تحسب التحويلين الخطيين الواردين في الخطوة الثالثة، ثم ترتب الناتج بحيث تأتي القيمة الأكبر أولا، وتحذف القيم الصفرية، وتتجاهل الفرع الثاني إذا كان مكررا للأول، وتبقي فقط الأزواج المتباينة أوليا.
كل زوج ناجح يحول إلى \((x,y,z)=(qr,pr,pq)\). وإذا كانت القيم الثلاث كلها لا تتجاوز \(N\)، يضاف \(x+y+z\) إلى المجموع الجاري. كما تستخدم الأنواع العددية المناسبة حتى تبقى الضربات الوسيطة وعملية التجميع صحيحة تماما.
لنكتب
$$R=\left\lfloor\left(\frac{2N^2}{13}\right)^{1/4}\right\rfloor.$$
مجال البحث يزور فقط الأزواج الأولية \((m,n)\) التي تحقق \(m^2+n^2\le R\)، ولذلك يكون عدد المرشحين \(O(R)\). وبما أن \(R=\Theta(\sqrt{N})\)، فإن الزمن الكلي هو \(O(\sqrt{N})\). أما الذاكرة الإضافية فتبقى \(O(1)\) باستثناء ما يلزم من حساب صحيح آمن للضرب والتجميع.