We must count the number \(F(L)\) of positive integer triples \((x,y,n)\) satisfying
$$x \lt y \le L,\qquad \frac1x+\frac1y=\frac1n.$$
The symmetry between \(x\) and \(y\) means that the strict inequality \(x \lt y\) removes the mirror-image duplicate of every solution. The published checkpoints are \(F(15)=4\) and \(F(1000)=1069\), while the real target is \(F(10^{12})\).
Let
$$g=\gcd(x,y),\qquad x=g\,u,\qquad y=g\,v,$$
with \(u \lt v\) and \(\gcd(u,v)=1\). Substituting into the reciprocal equation gives
$$n(u+v)=g\,u\,v.$$
Because \(\gcd(u,u+v)=1\) and \(\gcd(v,u+v)=1\), we can read divisibility information directly from this identity. Since \(u\mid n(u+v)\) and \(u\) is coprime to \(u+v\), it follows that \(u\mid n\). The same argument gives \(v\mid n\). As \(\gcd(u,v)=1\), we conclude that \(uv\mid n\), so we may write
$$n=tuv,\qquad t\in\mathbb Z_{>0}.$$
Then the previous equation simplifies to \(g=t(u+v)\). Therefore every admissible solution has the form
$$x=t\,u(u+v),\qquad y=t\,v(u+v),\qquad n=t\,u\,v,$$
with \(u,v,t \in \mathbb Z_{>0}\), \(u \lt v\), and \(\gcd(u,v)=1\). Conversely, every triple of this form satisfies \(\frac1x+\frac1y=\frac1n\), so this is a bijection and not just a one-way construction.
Introduce the new variables
$$a=v,\qquad s=u+v.$$
Then \(u=s-a\). The condition \(1\le u \lt v\) becomes
$$a+1\le s\le 2a-1.$$
The coprimality condition is preserved because
$$\gcd(a,s)=\gcd(v,u+v)=\gcd(u,v)=1.$$
The upper bound on \(y\) is now
$$y=t\,a\,s\le L.$$
So for a fixed coprime pair \((a,s)\), the number of admissible scale factors \(t\) is exactly
$$1\le t\le \left\lfloor\frac{L}{as}\right\rfloor,$$
which contributes \(\left\lfloor L/(as)\right\rfloor\) solutions. Hence
$$F(L)=\sum_{a=1}^{A_{\max}}\ \sum_{s=a+1}^{U(a)} \mathbf{1}_{\gcd(a,s)=1}\left\lfloor\frac{L}{as}\right\rfloor,$$
where
$$U(a)=\min\left(2a-1,\left\lfloor\frac{L}{a}\right\rfloor\right).$$
The outer loop can stop once no valid \(s\) exists. Since the smallest allowed \(s\) is \(a+1\), we need \(a(a+1)\le L\). Therefore
$$A_{\max}=\left\lfloor\frac{\sqrt{1+4L}-1}{2}\right\rfloor.$$
For a fixed outer parameter \(a\), define
$$C_a(\alpha,\beta)=\#\left\{s\in\mathbb Z:\alpha\le s\le \beta,\ \gcd(a,s)=1\right\}.$$
The standard Möbius identity is
$$\mathbf{1}_{\gcd(a,s)=1}=\sum_{d\mid \gcd(a,s)}\mu(d).$$
Summing this over an interval gives
$$C_a(\alpha,\beta)=\sum_{d\mid a}\mu(d)\left(\left\lfloor\frac{\beta}{d}\right\rfloor-\left\lfloor\frac{\alpha-1}{d}\right\rfloor\right).$$
Only squarefree divisors contribute, because \(\mu(d)=0\) whenever \(d\) contains a repeated prime factor. The implementation therefore precomputes, for each outer value, exactly those divisors with nonzero Möbius sign and reuses them for every interval query.
Fix \(a\) and set
$$M=\left\lfloor\frac{L}{a}\right\rfloor,\qquad q(s)=\left\lfloor\frac{M}{s}\right\rfloor.$$
As \(s\) increases, \(q(s)\) remains constant across contiguous ranges. If the current block starts at \(\ell\) and
$$q=\left\lfloor\frac{M}{\ell}\right\rfloor,$$
then the same quotient remains valid up to
$$r=\min\left(U(a),\left\lfloor\frac{M}{q}\right\rfloor\right).$$
Every \(s\in[\ell,r]\) contributes the same multiplier \(q\), so the whole block contributes
$$q\cdot C_a(\ell,r).$$
This is the main speedup: instead of scanning each inner value one by one, the algorithm jumps directly between maximal ranges on which the floor quotient is constant.
Here
$$A_{\max}=\left\lfloor\frac{\sqrt{61}-1}{2}\right\rfloor=3.$$
For \(a=1\), the interval \(a+1\le s\le 2a-1\) is empty, so there is no contribution.
For \(a=2\), the only possible value is \(s=3\). It is coprime to \(2\), and it contributes
$$\left\lfloor\frac{15}{2\cdot 3}\right\rfloor=2.$$
For \(a=3\), the admissible values are \(s=4\) and \(s=5\). Both are coprime to \(3\), giving
$$\left\lfloor\frac{15}{3\cdot 4}\right\rfloor=1,\qquad \left\lfloor\frac{15}{3\cdot 5}\right\rfloor=1.$$
Therefore
$$F(15)=2+1+1=4,$$
which matches the checkpoint exactly.
The C++, Python, and Java implementations all follow this derivation. They first compute \(A_{\max}\) using an integer square root so the outer range is exact. Then they build the Möbius function up to that limit with a sieve, and from it they assemble, for every outer value, the signed list of squarefree divisors needed by the interval formula above.
During the main summation, the implementation determines the valid inner interval, advances through maximal quotient blocks, and for each block evaluates the coprime count using the precomputed signed divisors. The C++ and Java versions store this data in compact flattened form to reduce allocation overhead, while the Python version keeps direct per-value lists. The Java implementation also parallelizes independent outer ranges across available CPU cores.
Let \(A=A_{\max}=O(\sqrt L)\). Computing the Möbius function is sieve-like and near-linear in \(A\). Building the stored squarefree-divisor data costs
$$\sum_{\substack{d\le A\\ \mu(d)\ne 0}}\left\lfloor\frac{A}{d}\right\rfloor=O(A\log A),$$
and the memory needed for that data has the same order.
The summation stage is much smaller than a naive search over all \(x\) and \(y\). It processes quotient blocks rather than individual inner values, and each block performs one Möbius interval count over the stored divisors of the current outer value. In practice this behaves close to \(O(\sqrt L\log L)\), which is why the method remains feasible for \(L=10^{12}\), whereas direct enumeration would be hopelessly large.
Gesucht ist die Anzahl \(F(L)\) der positiven ganzzahligen Tripel \((x,y,n)\) mit
$$x \lt y \le L,\qquad \frac1x+\frac1y=\frac1n.$$
Da \(x\) und \(y\) symmetrisch auftreten, entfernt die Bedingung \(x \lt y\) genau die gespiegelte Doppelzählung jeder Lösung. Die Kontrollwerte sind \(F(15)=4\) und \(F(1000)=1069\); die eigentliche Aufgabe fragt nach \(F(10^{12})\).
Setze
$$g=\gcd(x,y),\qquad x=g\,u,\qquad y=g\,v,$$
mit \(u \lt v\) und \(\gcd(u,v)=1\). Einsetzen in die Reziprokgleichung liefert
$$n(u+v)=g\,u\,v.$$
Wegen \(\gcd(u,u+v)=1\) und \(\gcd(v,u+v)=1\) kann man sofort Teilbarkeiten ablesen. Aus \(u\mid n(u+v)\) und \(\gcd(u,u+v)=1\) folgt \(u\mid n\). Analog erhält man \(v\mid n\). Weil \(u\) und \(v\) teilerfremd sind, gilt also \(uv\mid n\). Daher schreiben wir
$$n=tuv,\qquad t\in\mathbb Z_{>0}.$$
Dann vereinfacht sich die Gleichung zu \(g=t(u+v)\). Jede zulässige Lösung besitzt somit die Form
$$x=t\,u(u+v),\qquad y=t\,v(u+v),\qquad n=t\,u\,v,$$
mit \(u,v,t\in\mathbb Z_{>0}\), \(u \lt v\) und \(\gcd(u,v)=1\). Umgekehrt erfüllt jedes Tripel dieser Form die Gleichung \(\frac1x+\frac1y=\frac1n\). Die Parametrisierung ist also bijektiv.
Führe die Variablen
$$a=v,\qquad s=u+v$$
ein. Dann ist \(u=s-a\). Die Bedingung \(1\le u \lt v\) wird zu
$$a+1\le s\le 2a-1.$$
Die Teilerfremdheit bleibt erhalten, denn
$$\gcd(a,s)=\gcd(v,u+v)=\gcd(u,v)=1.$$
Die Schranke für \(y\) lautet jetzt
$$y=t\,a\,s\le L.$$
Für ein festes teilerfremdes Paar \((a,s)\) gibt es also genau
$$1\le t\le \left\lfloor\frac{L}{as}\right\rfloor$$
zulässige Skalierungsfaktoren. Daher gilt
$$F(L)=\sum_{a=1}^{A_{\max}}\ \sum_{s=a+1}^{U(a)} \mathbf{1}_{\gcd(a,s)=1}\left\lfloor\frac{L}{as}\right\rfloor,$$
wobei
$$U(a)=\min\left(2a-1,\left\lfloor\frac{L}{a}\right\rfloor\right).$$
Der äußere Bereich endet, sobald es kein zulässiges \(s\) mehr gibt. Da das kleinste mögliche \(s\) gleich \(a+1\) ist, brauchen wir \(a(a+1)\le L\). Also
$$A_{\max}=\left\lfloor\frac{\sqrt{1+4L}-1}{2}\right\rfloor.$$
Für festes \(a\) definieren wir
$$C_a(\alpha,\beta)=\#\left\{s\in\mathbb Z:\alpha\le s\le \beta,\ \gcd(a,s)=1\right\}.$$
Die klassische Möbius-Identität lautet
$$\mathbf{1}_{\gcd(a,s)=1}=\sum_{d\mid \gcd(a,s)}\mu(d).$$
Durch Summation über ein Intervall erhält man
$$C_a(\alpha,\beta)=\sum_{d\mid a}\mu(d)\left(\left\lfloor\frac{\beta}{d}\right\rfloor-\left\lfloor\frac{\alpha-1}{d}\right\rfloor\right).$$
Nur quadratfreie Teiler tragen bei, denn \(\mu(d)=0\), sobald \(d\) einen mehrfachen Primfaktor besitzt. Genau deshalb speichert die Implementierung für jeden äußeren Wert nur die Teiler mit nichtverschwindendem Möbius-Vorzeichen und verwendet sie für alle Intervallabfragen wieder.
Fixiere \(a\) und setze
$$M=\left\lfloor\frac{L}{a}\right\rfloor,\qquad q(s)=\left\lfloor\frac{M}{s}\right\rfloor.$$
Wenn \(s\) wächst, bleibt \(q(s)\) auf ganzen Bereichen konstant. Beginnt der aktuelle Block bei \(\ell\) und gilt
$$q=\left\lfloor\frac{M}{\ell}\right\rfloor,$$
dann reicht derselbe Quotient bis
$$r=\min\left(U(a),\left\lfloor\frac{M}{q}\right\rfloor\right).$$
Der gesamte Block \([\ell,r]\) liefert somit den Beitrag
$$q\cdot C_a(\ell,r).$$
Das ist der zentrale Beschleunigungsschritt: Statt jedes innere \(s\) einzeln zu untersuchen, springt der Algorithmus direkt zwischen maximalen Bereichen mit konstantem Floor-Quotienten.
Hier ist
$$A_{\max}=\left\lfloor\frac{\sqrt{61}-1}{2}\right\rfloor=3.$$
Für \(a=1\) ist das Intervall \(a+1\le s\le 2a-1\) leer, also gibt es keinen Beitrag.
Für \(a=2\) ist nur \(s=3\) möglich. Dieser Wert ist teilerfremd zu \(2\) und trägt bei mit
$$\left\lfloor\frac{15}{2\cdot 3}\right\rfloor=2.$$
Für \(a=3\) sind \(s=4\) und \(s=5\) zulässig. Beide sind teilerfremd zu \(3\), also erhält man
$$\left\lfloor\frac{15}{3\cdot 4}\right\rfloor=1,\qquad \left\lfloor\frac{15}{3\cdot 5}\right\rfloor=1.$$
Folglich
$$F(15)=2+1+1=4,$$
genau wie im Kontrollwert.
Die C++-, Python- und Java-Implementierungen folgen genau dieser Ableitung. Zuerst bestimmen sie \(A_{\max}\) mit einer ganzzahligen Quadratwurzel, damit der äußere Bereich exakt endet. Danach wird die Möbius-Funktion bis zu dieser Grenze per Sieb berechnet, und daraus werden für jeden äußeren Wert die vorzeichenbehafteten quadratfreien Teiler aufgebaut, die in der Intervallformel benötigt werden.
Während der Hauptsumme ermittelt die Implementierung das zulässige Intervall für den inneren Parameter, springt blockweise über identische Quotienten und berechnet pro Block die Zahl teilerfremder Werte mit den vorberechneten signierten Teilern. Die C++- und Java-Versionen speichern diese Daten in kompakter abgeflachter Form, während Python direkte Listen pro äußerem Wert verwendet. Die Java-Implementierung verteilt den äußeren Bereich zusätzlich parallel auf mehrere Kerne.
Sei \(A=A_{\max}=O(\sqrt L)\). Die Berechnung der Möbius-Funktion ist siebartig und nahezu linear in \(A\). Das Aufbauen der gespeicherten quadratfreien Teiler kostet
$$\sum_{\substack{d\le A\\ \mu(d)\ne 0}}\left\lfloor\frac{A}{d}\right\rfloor=O(A\log A),$$
und dieselbe Größenordnung gilt für den Speicherbedarf dieser Daten.
Die Summationsphase ist wesentlich kleiner als eine naive Suche über alle \(x\) und \(y\). Sie verarbeitet Quotientenblöcke statt einzelner innerer Werte, und jeder Block benötigt nur eine Möbius-Intervallauswertung über die gespeicherten Teiler des aktuellen äußeren Werts. Praktisch verhält sich das Verfahren nahe bei \(O(\sqrt L\log L)\). Genau deshalb bleibt \(L=10^{12}\) handhabbar, während direkte Enumeration völlig unbrauchbar wäre.
Aranan şey,
$$x \lt y \le L,\qquad \frac1x+\frac1y=\frac1n$$
koşullarını sağlayan pozitif tamsayı üçlülerinin \((x,y,n)\) sayısıdır. \(x\) ile \(y\) simetrik olduğu için \(x \lt y\) şartı her çözümün ayna görüntüsünü ikinci kez saymayı engeller. Kontrol değerleri \(F(15)=4\) ve \(F(1000)=1069\) olup asıl hedef \(F(10^{12})\) değeridir.
Şöyle yazalım:
$$g=\gcd(x,y),\qquad x=g\,u,\qquad y=g\,v,$$
burada \(u \lt v\) ve \(\gcd(u,v)=1\). Bunu denklemde yerine koyunca
$$n(u+v)=g\,u\,v$$
elde edilir. \(\gcd(u,u+v)=1\) ve \(\gcd(v,u+v)=1\) olduğundan bölüm özellikleri doğrudan çıkar. \(u\mid n(u+v)\) ve \(u\) ile \(u+v\) aralarında asal olduğu için \(u\mid n\). Aynı şekilde \(v\mid n\). Ayrıca \(\gcd(u,v)=1\) olduğundan \(uv\mid n\) sonucuna varırız. Bu nedenle
$$n=tuv,\qquad t\in\mathbb Z_{>0}$$
yazabiliriz. Böylece denklem \(g=t(u+v)\) olur ve her çözüm
$$x=t\,u(u+v),\qquad y=t\,v(u+v),\qquad n=t\,u\,v$$
biçimindedir; burada \(u,v,t\in\mathbb Z_{>0}\), \(u \lt v\) ve \(\gcd(u,v)=1\). Tersi de doğrudur: bu formdaki her üçlü \(\frac1x+\frac1y=\frac1n\) eşitliğini sağlar. Yani bu parametrizasyon bire bir eşlemedir.
Yeni değişkenleri
$$a=v,\qquad s=u+v$$
olarak tanımlayalım. O zaman \(u=s-a\) olur. \(1\le u \lt v\) koşulu
$$a+1\le s\le 2a-1$$
şeklini alır. Aralarında asallık da korunur, çünkü
$$\gcd(a,s)=\gcd(v,u+v)=\gcd(u,v)=1.$$
\(y\) için üst sınır artık
$$y=t\,a\,s\le L$$
olur. Sabit bir aralarında asal \((a,s)\) çifti için izin verilen ölçekleme katsayısı sayısı tam olarak
$$1\le t\le \left\lfloor\frac{L}{as}\right\rfloor$$
kadardır. Bu yüzden
$$F(L)=\sum_{a=1}^{A_{\max}}\ \sum_{s=a+1}^{U(a)} \mathbf{1}_{\gcd(a,s)=1}\left\lfloor\frac{L}{as}\right\rfloor,$$
burada
$$U(a)=\min\left(2a-1,\left\lfloor\frac{L}{a}\right\rfloor\right).$$
Dış döngü, artık geçerli bir \(s\) kalmadığında durabilir. En küçük olası \(s\) değeri \(a+1\) olduğundan \(a(a+1)\le L\) gereklidir. Dolayısıyla
$$A_{\max}=\left\lfloor\frac{\sqrt{1+4L}-1}{2}\right\rfloor.$$
Sabit \(a\) için
$$C_a(\alpha,\beta)=\#\left\{s\in\mathbb Z:\alpha\le s\le \beta,\ \gcd(a,s)=1\right\}$$
tanımını yapalım. Standart Möbius özdeşliği
$$\mathbf{1}_{\gcd(a,s)=1}=\sum_{d\mid \gcd(a,s)}\mu(d)$$
şeklindedir. Bunu aralık üzerinde toplarsak
$$C_a(\alpha,\beta)=\sum_{d\mid a}\mu(d)\left(\left\lfloor\frac{\beta}{d}\right\rfloor-\left\lfloor\frac{\alpha-1}{d}\right\rfloor\right)$$
elde edilir. Yalnızca kareden arınmış bölenler katkı verir; çünkü bir asalın karesini içeren \(d\) için \(\mu(d)=0\) olur. Bu yüzden implementasyon, her dış değer için yalnızca Möbius işareti sıfır olmayan bölenleri önceden hazırlar ve bütün aralık sorgularında tekrar kullanır.
Sabit \(a\) için
$$M=\left\lfloor\frac{L}{a}\right\rfloor,\qquad q(s)=\left\lfloor\frac{M}{s}\right\rfloor$$
olsun. \(s\) arttıkça \(q(s)\) uzun bitişik aralıklarda sabit kalır. Geçerli bloğun sol ucu \(\ell\) ise ve
$$q=\left\lfloor\frac{M}{\ell}\right\rfloor$$
ise, aynı bölüm değeri şu sınıra kadar sürer:
$$r=\min\left(U(a),\left\lfloor\frac{M}{q}\right\rfloor\right).$$
Dolayısıyla tüm \([\ell,r]\) bloğunun katkısı
$$q\cdot C_a(\ell,r)$$
olur. Esas hızlandırma budur: algoritma her iç değeri tek tek gezmek yerine, taban bölümün değişmediği en büyük aralıklara sıçrar.
Bu durumda
$$A_{\max}=\left\lfloor\frac{\sqrt{61}-1}{2}\right\rfloor=3.$$
\(a=1\) için \(a+1\le s\le 2a-1\) aralığı boştur; katkı yoktur.
\(a=2\) için tek olası değer \(s=3\) olur. Bu sayı \(2\) ile aralarında asaldır ve
$$\left\lfloor\frac{15}{2\cdot 3}\right\rfloor=2$$
katkısını verir.
\(a=3\) için geçerli değerler \(s=4\) ve \(s=5\) tir. Her ikisi de \(3\) ile aralarında asaldır ve katkılar
$$\left\lfloor\frac{15}{3\cdot 4}\right\rfloor=1,\qquad \left\lfloor\frac{15}{3\cdot 5}\right\rfloor=1$$
olur. Sonuç olarak
$$F(15)=2+1+1=4,$$
ki bu da kontrol değeriyle tam uyumludur.
C++, Python ve Java implementasyonları bu matematiksel planı doğrudan uygular. Önce tam sayı karekök ile \(A_{\max}\) hesaplanır; böylece dış aralık tam doğru yerde biter. Sonra bu sınıra kadar Möbius fonksiyonu bir süzgeç ile üretilir ve buradan, aralık formülünde kullanılacak işaretli kareden arınmış bölen listeleri her dış değer için hazırlanır.
Ana toplam sırasında implementasyon geçerli iç aralığı belirler, sabit bölüm blokları arasında ilerler ve her blok için aralarında asal sayı adedini önceden saklanan işaretli bölenlerle hesaplar. C++ ve Java sürümleri bu veriyi daha az bellek tahsisi için düzleştirilmiş kompakt yapıda tutar; Python sürümü ise doğrudan liste yaklaşımı kullanır. Java implementasyonu ayrıca bağımsız dış aralıkları işlemci çekirdekleri arasında paralelleştirir.
\(A=A_{\max}=O(\sqrt L)\) olsun. Möbius fonksiyonunu üretmek süzgeç temelli ve \(A\) cinsinden neredeyse doğrusaldır. Saklanan kareden arınmış bölen verisini kurmanın maliyeti
$$\sum_{\substack{d\le A\\ \mu(d)\ne 0}}\left\lfloor\frac{A}{d}\right\rfloor=O(A\log A)$$
mertebesindedir ve bellek ihtiyacı da aynı düzeydedir.
Toplama aşaması, tüm \(x\) ve \(y\) çiftlerini dolaşan saf yaklaşımdan çok daha küçüktür. İç değerler tek tek işlenmez; bunun yerine bölüm blokları kullanılır ve her blokta yalnızca ilgili dış değerin saklanan bölenleri üzerinde bir Möbius aralık sayımı yapılır. Pratikte davranış \(O(\sqrt L\log L)\) civarındadır. Bu yüzden \(L=10^{12}\) erişilebilir kalır; doğrudan enumerasyon ise tamamen kullanışsız olurdu.
Debemos contar el número \(F(L)\) de ternas de enteros positivos \((x,y,n)\) que cumplen
$$x \lt y \le L,\qquad \frac1x+\frac1y=\frac1n.$$
Como \(x\) e \(y\) aparecen de forma simétrica, la condición estricta \(x \lt y\) elimina la solución reflejada y evita el doble conteo. Los puntos de control publicados son \(F(15)=4\) y \(F(1000)=1069\); el objetivo real es \(F(10^{12})\).
Sea
$$g=\gcd(x,y),\qquad x=g\,u,\qquad y=g\,v,$$
con \(u \lt v\) y \(\gcd(u,v)=1\). Al sustituir en la ecuación recíproca obtenemos
$$n(u+v)=g\,u\,v.$$
Como \(\gcd(u,u+v)=1\) y \(\gcd(v,u+v)=1\), las divisibilidades salen de inmediato. De \(u\mid n(u+v)\) y \(\gcd(u,u+v)=1\) se deduce \(u\mid n\). Del mismo modo, \(v\mid n\). Puesto que \(\gcd(u,v)=1\), concluimos que \(uv\mid n\), así que podemos escribir
$$n=tuv,\qquad t\in\mathbb Z_{>0}.$$
Entonces la identidad se reduce a \(g=t(u+v)\). Por tanto, toda solución válida tiene la forma
$$x=t\,u(u+v),\qquad y=t\,v(u+v),\qquad n=t\,u\,v,$$
con \(u,v,t\in\mathbb Z_{>0}\), \(u \lt v\) y \(\gcd(u,v)=1\). La implicación inversa también vale: toda terna de esta forma satisface \(\frac1x+\frac1y=\frac1n\). Es una parametrización biyectiva.
Introducimos
$$a=v,\qquad s=u+v.$$
Entonces \(u=s-a\). La condición \(1\le u \lt v\) se convierte en
$$a+1\le s\le 2a-1.$$
La coprimalidad se conserva porque
$$\gcd(a,s)=\gcd(v,u+v)=\gcd(u,v)=1.$$
La cota sobre \(y\) pasa a ser
$$y=t\,a\,s\le L.$$
Para un par coprimo fijo \((a,s)\), el número de factores de escala admisibles es exactamente
$$1\le t\le \left\lfloor\frac{L}{as}\right\rfloor,$$
de modo que dicho par aporta \(\left\lfloor L/(as)\right\rfloor\) soluciones. En consecuencia,
$$F(L)=\sum_{a=1}^{A_{\max}}\ \sum_{s=a+1}^{U(a)} \mathbf{1}_{\gcd(a,s)=1}\left\lfloor\frac{L}{as}\right\rfloor,$$
donde
$$U(a)=\min\left(2a-1,\left\lfloor\frac{L}{a}\right\rfloor\right).$$
El bucle exterior puede detenerse cuando ya no exista ningún \(s\) válido. Como el menor valor posible es \(a+1\), necesitamos \(a(a+1)\le L\). Por eso
$$A_{\max}=\left\lfloor\frac{\sqrt{1+4L}-1}{2}\right\rfloor.$$
Para \(a\) fijo, definimos
$$C_a(\alpha,\beta)=\#\left\{s\in\mathbb Z:\alpha\le s\le \beta,\ \gcd(a,s)=1\right\}.$$
La identidad clásica es
$$\mathbf{1}_{\gcd(a,s)=1}=\sum_{d\mid \gcd(a,s)}\mu(d).$$
Al sumarla en un intervalo obtenemos
$$C_a(\alpha,\beta)=\sum_{d\mid a}\mu(d)\left(\left\lfloor\frac{\beta}{d}\right\rfloor-\left\lfloor\frac{\alpha-1}{d}\right\rfloor\right).$$
Solo importan los divisores libres de cuadrados, porque \(\mu(d)=0\) cuando \(d\) contiene un factor primo repetido. Por eso la implementación precomputa, para cada valor exterior, únicamente los divisores con signo de Möbius no nulo y los reutiliza en todas las consultas de intervalo.
Fijemos \(a\) y pongamos
$$M=\left\lfloor\frac{L}{a}\right\rfloor,\qquad q(s)=\left\lfloor\frac{M}{s}\right\rfloor.$$
Cuando \(s\) aumenta, el cociente \(q(s)\) permanece constante en intervalos contiguos. Si el bloque actual empieza en \(\ell\) y
$$q=\left\lfloor\frac{M}{\ell}\right\rfloor,$$
entonces ese mismo valor se mantiene hasta
$$r=\min\left(U(a),\left\lfloor\frac{M}{q}\right\rfloor\right).$$
Todo el bloque \([\ell,r]\) aporta por tanto
$$q\cdot C_a(\ell,r).$$
Esta es la aceleración principal: el algoritmo no recorre cada valor interior de uno en uno, sino que salta entre los rangos máximos donde el cociente entero es constante.
Aquí
$$A_{\max}=\left\lfloor\frac{\sqrt{61}-1}{2}\right\rfloor=3.$$
Para \(a=1\), el intervalo \(a+1\le s\le 2a-1\) está vacío, así que no hay contribución.
Para \(a=2\), el único valor posible es \(s=3\). Es coprimo con \(2\) y aporta
$$\left\lfloor\frac{15}{2\cdot 3}\right\rfloor=2.$$
Para \(a=3\), los valores admisibles son \(s=4\) y \(s=5\). Ambos son coprimos con \(3\), de modo que contribuyen
$$\left\lfloor\frac{15}{3\cdot 4}\right\rfloor=1,\qquad \left\lfloor\frac{15}{3\cdot 5}\right\rfloor=1.$$
Por lo tanto,
$$F(15)=2+1+1=4,$$
exactamente el valor de comprobación.
Las implementaciones en C++, Python y Java siguen exactamente esta derivación. Primero calculan \(A_{\max}\) con una raíz cuadrada entera para que el rango exterior termine en el punto exacto. Después construyen la función de Möbius hasta ese límite mediante un cribado y, a partir de ella, preparan para cada valor exterior la lista con signo de divisores libres de cuadrados que necesita la fórmula de intervalos.
Durante la suma principal, la implementación determina el intervalo interior válido, avanza por bloques máximos de cociente constante y evalúa en cada bloque cuántos valores son coprimos con el exterior fijo usando los divisores ya preparados. Las versiones de C++ y Java almacenan esos datos en forma compacta y aplanada para reducir sobrecostes de memoria, mientras que la versión de Python mantiene una lista directa por valor exterior. La implementación en Java además reparte el rango exterior entre varios núcleos.
Sea \(A=A_{\max}=O(\sqrt L)\). El cálculo de la función de Möbius es de tipo cribado y prácticamente lineal en \(A\). Construir la estructura de divisores libres de cuadrados cuesta
$$\sum_{\substack{d\le A\\ \mu(d)\ne 0}}\left\lfloor\frac{A}{d}\right\rfloor=O(A\log A),$$
y el uso de memoria tiene el mismo orden.
La fase de suma es mucho menor que una búsqueda ingenua sobre todos los pares \((x,y)\). En lugar de procesar cada valor interior por separado, trabaja con bloques de cociente y en cada bloque realiza una cuenta de Möbius sobre los divisores almacenados del valor exterior actual. En la práctica el comportamiento queda cerca de \(O(\sqrt L\log L)\). Por eso \(L=10^{12}\) es tratable, mientras que una enumeración directa sería completamente inviable.
我们要求的是满足
$$x \lt y \le L,\qquad \frac1x+\frac1y=\frac1n$$
的正整数三元组 \((x,y,n)\) 的个数 \(F(L)\)。由于 \(x\) 与 \(y\) 在方程中是对称的，所以加入 \(x \lt y\) 正好去掉镜像重复计数。题目给出的校验值是 \(F(15)=4\) 和 \(F(1000)=1069\)，真正目标是 \(F(10^{12})\)。
设
$$g=\gcd(x,y),\qquad x=g\,u,\qquad y=g\,v,$$
其中 \(u \lt v\) 且 \(\gcd(u,v)=1\)。代入原方程可得
$$n(u+v)=g\,u\,v.$$
因为 \(\gcd(u,u+v)=1\) 且 \(\gcd(v,u+v)=1\)，所以可以直接读出整除关系。由 \(u\mid n(u+v)\) 且 \(u\) 与 \(u+v\) 互素，推出 \(u\mid n\)。同理 \(v\mid n\)。再结合 \(\gcd(u,v)=1\)，便得到 \(uv\mid n\)。因此可以写成
$$n=tuv,\qquad t\in\mathbb Z_{>0}.$$
于是原等式化为 \(g=t(u+v)\)。所以每个解都唯一对应于
$$x=t\,u(u+v),\qquad y=t\,v(u+v),\qquad n=t\,u\,v,$$
其中 \(u,v,t\in\mathbb Z_{>0}\)，\(u \lt v\)，且 \(\gcd(u,v)=1\)。反过来，这种形式的每个三元组也都满足 \(\frac1x+\frac1y=\frac1n\)。因此这是一一对应的参数化。
引入新变量
$$a=v,\qquad s=u+v.$$
则 \(u=s-a\)。条件 \(1\le u \lt v\) 等价于
$$a+1\le s\le 2a-1.$$
互素条件也保持不变，因为
$$\gcd(a,s)=\gcd(v,u+v)=\gcd(u,v)=1.$$
对 \(y\) 的上界现在变成
$$y=t\,a\,s\le L.$$
因此，对固定的互素对 \((a,s)\)，允许的缩放系数个数恰好是
$$1\le t\le \left\lfloor\frac{L}{as}\right\rfloor,$$
于是该对参数贡献 \(\left\lfloor L/(as)\right\rfloor\) 个解。于是有
$$F(L)=\sum_{a=1}^{A_{\max}}\ \sum_{s=a+1}^{U(a)} \mathbf{1}_{\gcd(a,s)=1}\left\lfloor\frac{L}{as}\right\rfloor,$$
其中
$$U(a)=\min\left(2a-1,\left\lfloor\frac{L}{a}\right\rfloor\right).$$
外层循环只需要做到仍然存在合法 \(s\) 为止。因为最小可能值是 \(s=a+1\)，所以必须满足 \(a(a+1)\le L\)。因此
$$A_{\max}=\left\lfloor\frac{\sqrt{1+4L}-1}{2}\right\rfloor.$$
对固定 \(a\)，定义
$$C_a(\alpha,\beta)=\#\left\{s\in\mathbb Z:\alpha\le s\le \beta,\ \gcd(a,s)=1\right\}.$$
经典恒等式为
$$\mathbf{1}_{\gcd(a,s)=1}=\sum_{d\mid \gcd(a,s)}\mu(d).$$
把它在区间上求和，就得到
$$C_a(\alpha,\beta)=\sum_{d\mid a}\mu(d)\left(\left\lfloor\frac{\beta}{d}\right\rfloor-\left\lfloor\frac{\alpha-1}{d}\right\rfloor\right).$$
只有平方因子自由的除数才有贡献，因为一旦 \(d\) 含有重复素因子，就有 \(\mu(d)=0\)。这就是为什么实现会为每个外层参数预先保存所有 Möbius 符号非零的除数，并在后续区间查询中反复复用。
固定 \(a\)，记
$$M=\left\lfloor\frac{L}{a}\right\rfloor,\qquad q(s)=\left\lfloor\frac{M}{s}\right\rfloor.$$
随着 \(s\) 增大，\(q(s)\) 会在连续区间上保持不变。若当前块的左端点为 \(\ell\)，且
$$q=\left\lfloor\frac{M}{\ell}\right\rfloor,$$
则该商值一直保持到
$$r=\min\left(U(a),\left\lfloor\frac{M}{q}\right\rfloor\right).$$
因此整个区间 \([\ell,r]\) 的贡献就是
$$q\cdot C_a(\ell,r).$$
这正是主要优化所在：程序不是逐个扫描内层参数，而是直接跳过整段商值不变的最大区间。
此时
$$A_{\max}=\left\lfloor\frac{\sqrt{61}-1}{2}\right\rfloor=3.$$
当 \(a=1\) 时，区间 \(a+1\le s\le 2a-1\) 为空，因此没有贡献。
当 \(a=2\) 时，只可能有 \(s=3\)。它与 \(2\) 互素，贡献为
$$\left\lfloor\frac{15}{2\cdot 3}\right\rfloor=2.$$
当 \(a=3\) 时，可取 \(s=4,5\)。这两个数都与 \(3\) 互素，于是贡献分别是
$$\left\lfloor\frac{15}{3\cdot 4}\right\rfloor=1,\qquad \left\lfloor\frac{15}{3\cdot 5}\right\rfloor=1.$$
所以
$$F(15)=2+1+1=4,$$
与校验值完全一致。
C++、Python 和 Java 实现都严格遵循上面的推导。首先通过整数平方根算出 \(A_{\max}\)，从而精确截断外层范围。然后用筛法求出直到该上界的 Möbius 函数，并据此为每个外层参数构造后续区间计数所需的带符号平方因子自由除数列表。
在主求和阶段，实现会先求出当前合法的内层区间，再按照整除商不变的最大分块前进，并在每个分块上利用预处理好的带符号除数来计算与当前外层值互素的个数。C++ 和 Java 版本把这些数据压成紧凑的扁平结构，以减少额外分配；Python 版本则为每个外层值保留直接列表。Java 版本还会把彼此独立的外层区间并行分发到多个 CPU 核心上。
记 \(A=A_{\max}=O(\sqrt L)\)。Möbius 函数的计算属于筛法过程，对 \(A\) 来说接近线性。构造并保存平方因子自由除数数据的代价为
$$\sum_{\substack{d\le A\\ \mu(d)\ne 0}}\left\lfloor\frac{A}{d}\right\rfloor=O(A\log A),$$
相应的内存需求也是同一数量级。
真正的求和阶段远小于朴素地遍历所有 \((x,y)\) 对。它不是逐个处理每个内层值，而是只处理整除商分块；每个分块只需要对当前外层参数保存的除数做一次 Möbius 区间计数。实际运行表现接近 \(O(\sqrt L\log L)\)。因此 \(L=10^{12}\) 仍然可计算，而直接枚举会大得无法接受。
Нужно вычислить число \(F(L)\) положительных целочисленных троек \((x,y,n)\), для которых
$$x \lt y \le L,\qquad \frac1x+\frac1y=\frac1n.$$
Так как \(x\) и \(y\) входят в уравнение симметрично, условие \(x \lt y\) просто убирает зеркальный дубликат каждого решения. Контрольные значения: \(F(15)=4\) и \(F(1000)=1069\); настоящая цель задачи — найти \(F(10^{12})\).
Положим
$$g=\gcd(x,y),\qquad x=g\,u,\qquad y=g\,v,$$
где \(u \lt v\) и \(\gcd(u,v)=1\). Подстановка в исходное уравнение даёт
$$n(u+v)=g\,u\,v.$$
Поскольку \(\gcd(u,u+v)=1\) и \(\gcd(v,u+v)=1\), отсюда сразу извлекаются условия делимости. Из \(u\mid n(u+v)\) и взаимной простоты \(u\) и \(u+v\) следует \(u\mid n\). Точно так же \(v\mid n\). А так как \(\gcd(u,v)=1\), получаем \(uv\mid n\). Значит, можно написать
$$n=tuv,\qquad t\in\mathbb Z_{>0}.$$
После этого равенство упрощается до \(g=t(u+v)\). Следовательно, всякое допустимое решение имеет вид
$$x=t\,u(u+v),\qquad y=t\,v(u+v),\qquad n=t\,u\,v,$$
где \(u,v,t\in\mathbb Z_{>0}\), \(u \lt v\), \(\gcd(u,v)=1\). Обратное тоже верно: любая тройка такого вида удовлетворяет \(\frac1x+\frac1y=\frac1n\). Это биективная параметризация всех решений.
Введём новые переменные
$$a=v,\qquad s=u+v.$$
Тогда \(u=s-a\). Условие \(1\le u \lt v\) превращается в
$$a+1\le s\le 2a-1.$$
Взаимная простота сохраняется, потому что
$$\gcd(a,s)=\gcd(v,u+v)=\gcd(u,v)=1.$$
Ограничение на \(y\) теперь имеет вид
$$y=t\,a\,s\le L.$$
Значит, для фиксированной взаимно простой пары \((a,s)\) число допустимых коэффициентов масштаба равно
$$1\le t\le \left\lfloor\frac{L}{as}\right\rfloor,$$
то есть такая пара даёт вклад \(\left\lfloor L/(as)\right\rfloor\). Поэтому
$$F(L)=\sum_{a=1}^{A_{\max}}\ \sum_{s=a+1}^{U(a)} \mathbf{1}_{\gcd(a,s)=1}\left\lfloor\frac{L}{as}\right\rfloor,$$
где
$$U(a)=\min\left(2a-1,\left\lfloor\frac{L}{a}\right\rfloor\right).$$
Внешний цикл можно закончить, как только допустимых \(s\) больше не остаётся. Поскольку минимально возможное значение равно \(a+1\), необходимо \(a(a+1)\le L\). Следовательно,
$$A_{\max}=\left\lfloor\frac{\sqrt{1+4L}-1}{2}\right\rfloor.$$
Для фиксированного \(a\) определим
$$C_a(\alpha,\beta)=\#\left\{s\in\mathbb Z:\alpha\le s\le \beta,\ \gcd(a,s)=1\right\}.$$
Классическое тождество имеет вид
$$\mathbf{1}_{\gcd(a,s)=1}=\sum_{d\mid \gcd(a,s)}\mu(d).$$
Просуммировав его по отрезку, получаем
$$C_a(\alpha,\beta)=\sum_{d\mid a}\mu(d)\left(\left\lfloor\frac{\beta}{d}\right\rfloor-\left\lfloor\frac{\alpha-1}{d}\right\rfloor\right).$$
Вклад дают только квадратсвободные делители, потому что \(\mu(d)=0\), если в \(d\) есть повторяющийся простой множитель. Поэтому реализация заранее хранит для каждого внешнего параметра только те делители, у которых знак Мёбиуса ненулевой, и затем многократно использует их во всех запросах на отрезке.
Зафиксируем \(a\) и обозначим
$$M=\left\lfloor\frac{L}{a}\right\rfloor,\qquad q(s)=\left\lfloor\frac{M}{s}\right\rfloor.$$
При росте \(s\) значение \(q(s)\) остаётся постоянным на целых contiguous-отрезках. Если текущий блок начинается в точке \(\ell\) и
$$q=\left\lfloor\frac{M}{\ell}\right\rfloor,$$
то этот же частный сохраняется до
$$r=\min\left(U(a),\left\lfloor\frac{M}{q}\right\rfloor\right).$$
Следовательно, вклад всего блока \([\ell,r]\) равен
$$q\cdot C_a(\ell,r).$$
Это и есть ключевое ускорение: вместо покоординатного перебора внутреннего параметра алгоритм сразу перескакивает между максимальными диапазонами, где целая часть не меняется.
В этом случае
$$A_{\max}=\left\lfloor\frac{\sqrt{61}-1}{2}\right\rfloor=3.$$
При \(a=1\) интервал \(a+1\le s\le 2a-1\) пуст, поэтому вклада нет.
При \(a=2\) возможен только \(s=3\). Он взаимно прост с \(2\) и даёт вклад
$$\left\lfloor\frac{15}{2\cdot 3}\right\rfloor=2.$$
При \(a=3\) допустимы \(s=4\) и \(s=5\). Оба значения взаимно просты с \(3\), поэтому дают
$$\left\lfloor\frac{15}{3\cdot 4}\right\rfloor=1,\qquad \left\lfloor\frac{15}{3\cdot 5}\right\rfloor=1.$$
Значит,
$$F(15)=2+1+1=4,$$
что точно совпадает с контрольным значением.
Реализации на C++, Python и Java напрямую следуют этому выводу. Сначала они вычисляют \(A_{\max}\) с помощью целочисленного квадратного корня, чтобы точно обрезать внешний диапазон. Затем до этого предела строится функция Мёбиуса с помощью решета, и на её основе для каждого внешнего значения подготавливается знаковый список квадратсвободных делителей, нужный для интервальной формулы.
На основном этапе суммирования реализация определяет допустимый внутренний диапазон, переходит между максимальными блоками с одинаковым частным и для каждого блока вычисляет число взаимно простых значений по заранее подготовленным знаковым делителям. Версии на C++ и Java держат эти данные в компактном уплощённом виде, чтобы уменьшить накладные расходы на память, а версия на Python хранит отдельный список для каждого внешнего значения. Java-реализация дополнительно распараллеливает независимые внешние диапазоны по доступным ядрам процессора.
Пусть \(A=A_{\max}=O(\sqrt L)\). Вычисление функции Мёбиуса имеет решеточный характер и близко к линейному по \(A\). Построение и хранение данных о квадратсвободных делителях стоит
$$\sum_{\substack{d\le A\\ \mu(d)\ne 0}}\left\lfloor\frac{A}{d}\right\rfloor=O(A\log A),$$
и такой же порядок имеет объём памяти.
Сама суммирующая часть намного меньше наивного перебора всех пар \((x,y)\). Вместо обработки каждого внутреннего значения по отдельности она работает по блокам одинакового частного, и в каждом блоке требуется лишь одна интервальная оценка по функции Мёбиуса на сохранённых делителях текущего внешнего значения. На практике поведение близко к \(O(\sqrt L\log L)\). Именно поэтому случай \(L=10^{12}\) остаётся вычислимым, тогда как прямой перебор был бы совершенно непригоден.
نريد حساب عدد الثلاثيات الصحيحة الموجبة \((x,y,n)\) التي تحقق
$$x \lt y \le L,\qquad \frac1x+\frac1y=\frac1n.$$
بما أن \(x\) و\(y\) يظهران بصورة متناظرة، فإن الشرط \(x \lt y\) يزيل النسخة المنعكسة من كل حل ولا يغيّر عدد الأزواج غير المرتبة. قيم التحقق المعطاة هي \(F(15)=4\) و\(F(1000)=1069\)، أما الهدف الفعلي فهو \(F(10^{12})\).
لنكتب
$$g=\gcd(x,y),\qquad x=g\,u,\qquad y=g\,v,$$
بحيث \(u \lt v\) و\(\gcd(u,v)=1\). بالتعويض في المعادلة الأصلية نحصل على
$$n(u+v)=g\,u\,v.$$
وبما أن \(\gcd(u,u+v)=1\) و\(\gcd(v,u+v)=1\)، يمكن استخراج شروط القسمة مباشرة. من \(u\mid n(u+v)\) ومع كون \(u\) أولياً مع \(u+v\)، نستنتج \(u\mid n\). وبنفس الطريقة \(v\mid n\). ومع \(\gcd(u,v)=1\) نحصل على \(uv\mid n\). لذلك يمكن كتابة
$$n=tuv,\qquad t\in\mathbb Z_{>0}.$$
وعندها تختزل العلاقة إلى \(g=t(u+v)\). إذن كل حل مقبول يأخذ الصورة
$$x=t\,u(u+v),\qquad y=t\,v(u+v),\qquad n=t\,u\,v,$$
حيث \(u,v,t\in\mathbb Z_{>0}\)، و\(u \lt v\)، و\(\gcd(u,v)=1\). والعكس صحيح أيضاً: كل ثلاثية من هذا الشكل تحقق \(\frac1x+\frac1y=\frac1n\). أي أن هذا التمثيل يعطي مطابقة تامة مع جميع الحلول.
نعرّف المتغيرين
$$a=v,\qquad s=u+v.$$
إذن \(u=s-a\). والشرط \(1\le u \lt v\) يصبح
$$a+1\le s\le 2a-1.$$
كما أن شرط الأولية النسبية يبقى محفوظاً، لأن
$$\gcd(a,s)=\gcd(v,u+v)=\gcd(u,v)=1.$$
أما قيد الحد الأعلى على \(y\) فيتحول إلى
$$y=t\,a\,s\le L.$$
ومن ثم، إذا ثبّتنا الزوج النسبي الأولي \((a,s)\)، فإن عدد معاملات التوسيع الممكنة يساوي تماماً
$$1\le t\le \left\lfloor\frac{L}{as}\right\rfloor.$$
وبالتالي يساهم هذا الزوج بمقدار \(\left\lfloor L/(as)\right\rfloor\) من الحلول. لذا نحصل على
$$F(L)=\sum_{a=1}^{A_{\max}}\ \sum_{s=a+1}^{U(a)} \mathbf{1}_{\gcd(a,s)=1}\left\lfloor\frac{L}{as}\right\rfloor,$$
حيث
$$U(a)=\min\left(2a-1,\left\lfloor\frac{L}{a}\right\rfloor\right).$$
يمكن إيقاف الحلقة الخارجية عندما لا يعود هناك أي \(s\) صالح. وبما أن أصغر قيمة ممكنة هي \(a+1\)، فلا بد من تحقق \(a(a+1)\le L\). لذلك
$$A_{\max}=\left\lfloor\frac{\sqrt{1+4L}-1}{2}\right\rfloor.$$
لـ \(a\) ثابت، نعرّف
$$C_a(\alpha,\beta)=\#\left\{s\in\mathbb Z:\alpha\le s\le \beta,\ \gcd(a,s)=1\right\}.$$
والهوية القياسية هي
$$\mathbf{1}_{\gcd(a,s)=1}=\sum_{d\mid \gcd(a,s)}\mu(d).$$
وبجمعها على المجال نحصل على
$$C_a(\alpha,\beta)=\sum_{d\mid a}\mu(d)\left(\left\lfloor\frac{\beta}{d}\right\rfloor-\left\lfloor\frac{\alpha-1}{d}\right\rfloor\right).$$
فقط القواسم الخالية من المربعات تساهم في هذا الجمع، لأن \(\mu(d)=0\) متى احتوى \(d\) على عامل أولي مكرر. ولهذا السبب تقوم التنفيذات بتهيئة القواسم ذات إشارة Möbius غير الصفرية لكل قيمة خارجية، ثم تعيد استخدامها في جميع استعلامات المجالات.
ثبّت \(a\)، ثم ضع
$$M=\left\lfloor\frac{L}{a}\right\rfloor,\qquad q(s)=\left\lfloor\frac{M}{s}\right\rfloor.$$
عندما يكبر \(s\)، تبقى القيمة \(q(s)\) ثابتة على فترات متجاورة كاملة. إذا بدأ المقطع الحالي عند \(\ell\) وكان
$$q=\left\lfloor\frac{M}{\ell}\right\rfloor,$$
فإن هذا الخارج نفسه يستمر حتى
$$r=\min\left(U(a),\left\lfloor\frac{M}{q}\right\rfloor\right).$$
وعليه فمساهمة الكتلة كلها \([\ell,r]\) هي
$$q\cdot C_a(\ell,r).$$
وهذه هي نقطة التسريع الأساسية: بدلاً من المرور على كل قيمة داخلية على حدة، تقفز الخوارزمية بين أكبر المجالات التي يبقى فيها خارج القسمة ثابتاً.
لدينا هنا
$$A_{\max}=\left\lfloor\frac{\sqrt{61}-1}{2}\right\rfloor=3.$$
عندما \(a=1\)، يكون المجال \(a+1\le s\le 2a-1\) فارغاً، فلا توجد مساهمة.
وعندما \(a=2\)، تكون القيمة الوحيدة الممكنة هي \(s=3\). وهي أولية نسبياً مع \(2\)، لذا تساهم بمقدار
$$\left\lfloor\frac{15}{2\cdot 3}\right\rfloor=2.$$
وعندما \(a=3\)، تكون القيم المسموح بها هي \(s=4\) و\(s=5\). وكلتاهما أولية نسبياً مع \(3\)، ولذلك نحصل على
$$\left\lfloor\frac{15}{3\cdot 4}\right\rfloor=1,\qquad \left\lfloor\frac{15}{3\cdot 5}\right\rfloor=1.$$
ومن ثم
$$F(15)=2+1+1=4,$$
وهو تماماً مقدار التحقق المعطى.
تتبع التنفيذات في C++ وPython وJava هذا الاشتقاق مباشرة. في البداية تُحسب \(A_{\max}\) باستخدام جذر تربيعي صحيح حتى يتوقف المجال الخارجي في الموضع الصحيح تماماً. بعد ذلك تُبنى دالة Möbius حتى ذلك الحد بواسطة غربال، ومن هذه القيم تُجهَّز لكل قيمة خارجية قائمة موقعة من القواسم الخالية من المربعات اللازمة لصيغة العد على المجالات.
أثناء الجمع الرئيسي، تحدد الشيفرة المجال الداخلي الصالح، ثم تنتقل على شكل كتل ذات خارج قسمة ثابت، وتحسب في كل كتلة عدد القيم الأولية نسبياً باستعمال القواسم الموقعة المحضّرة مسبقاً. نسختا C++ وJava تخزنان هذه البيانات في بنية مسطحة مضغوطة لتقليل كلفة التخصيص، بينما تحتفظ نسخة Python بقائمة مباشرة لكل قيمة خارجية. كما أن تنفيذ Java يوزّع المجالات الخارجية المستقلة على عدة أنوية معالجة.
لنرمز إلى \(A=A_{\max}=O(\sqrt L)\). حساب دالة Möbius ذو طبيعة غربالية وقريب من الخطي بالنسبة إلى \(A\). أما بناء بيانات القواسم الخالية من المربعات فيكلف
$$\sum_{\substack{d\le A\\ \mu(d)\ne 0}}\left\lfloor\frac{A}{d}\right\rfloor=O(A\log A),$$
ولهذا الترتيب نفسه تقريباً متطلب الذاكرة.
مرحلة الجمع أصغر بكثير من بحث ساذج على جميع الأزواج \((x,y)\). فهي لا تمر على كل قيمة داخلية منفردة، بل تعمل على كتل خارج قسمة، وفي كل كتلة يكفي إجراء عدّ Möbius على القواسم المخزنة للقيمة الخارجية الحالية. عملياً يتصرف هذا قريباً من \(O(\sqrt L\log L)\). ولهذا تبقى حالة \(L=10^{12}\) قابلة للحساب، في حين أن التعداد المباشر سيكون ضخماً إلى حد غير عملي.