For positive integers \(x\) and \(y\), let \(E(x,y)\) denote the number of divisions performed by Euclid's algorithm until the remainder becomes \(0\). The problem asks for
$$S(N)=\sum_{1\le x,y\le N} E(x,y)$$
at the very large value \(N=5\cdot 10^6\). A direct double loop is hopeless: there are \(N^2\) ordered pairs, and each pair would still require a separate Euclidean run. The solution therefore reorganizes the sum so that most of the work is done by arithmetic formulas, Möbius inversion, and floor-sum recurrences.
When \(x>y\), starting from \((y,x)\) performs one swap before entering exactly the same remainder chain, so
$$E(y,x)=E(x,y)+1 \qquad (x>y).$$
Also, \(E(x,x)=1\) for every diagonal pair. If we define
$$L(N)=\sum_{1\le y<x\le N} E(x,y),$$
then the whole ordered sum becomes
$$\begin{aligned} S(N) &=\sum_{x=1}^{N}E(x,x)+\sum_{1\le y<x\le N}\bigl(E(x,y)+E(y,x)\bigr) \\ &=N+\sum_{1\le y<x\le N}\bigl(2E(x,y)+1\bigr) \\ &=2L(N)+\frac{N(N+1)}{2}. \end{aligned}$$
So it is enough to evaluate the lower-triangular sum \(L(N)\).
Every pair with \(x>y\) needs at least one Euclidean step. There is also an easily recognized second step: if \(x<2y\), then the first quotient is \(1\), so after one subtraction we are left with \((y,x-y)\), which is still a non-terminal pair. Hence
$$E(x,y)=1+\mathbf{1}_{x<2y}+R(x,y),\qquad R(x,y)\ge 0.$$
The easy contribution is therefore
$$B(N)=\sum_{1\le y<x\le N}\left(1+\mathbf{1}_{x<2y}\right).$$
The first part is just the number of pairs below the diagonal:
$$\sum_{1\le y<x\le N}1=\frac{N(N-1)}{2}.$$
For the indicator term, count the pairs with \(y<x<2y\). Writing \(N=2m\) or \(N=2m+1\), one gets
$$\sum_{1\le y<x\le 2m}\mathbf{1}_{x<2y}=m(m-1),$$
$$\sum_{1\le y<x\le 2m+1}\mathbf{1}_{x<2y}=m^2.$$
So the closed base term is
$$B(N)= \begin{cases} (3m-2)m, & N=2m, \\ 3m^2+m, & N=2m+1. \end{cases}$$
This is exactly the closed form used by the implementations before any expensive computation begins.
The remaining term \(R(x,y)\) measures everything beyond those universal first one or two steps. Because Euclid's algorithm is unchanged by scaling both arguments by the same factor,
$$R(dx,dy)=R(x,y)\qquad (d\ge 1).$$
That makes the gcd decomposition natural: the residual depends only on the primitive core of the pair.
Another useful observation is that \(R(x,y)=0\) unless the larger coordinate is at least \(5\). Indeed, the first primitive pairs with more than the peeled-off base behavior are \((5,2)\) and \((5,3)\). This is why the Möbius layer only needs to run up to \(\lfloor N/5\rfloor\).
Let \(\mu\) be the Möbius function and let
$$M(t)=\sum_{d\le t}\mu(d)$$
be its prefix sum. The fast method isolates primitive pairs by Möbius inversion and rewrites the residual part of the lower triangle in the form
$$L(N)=B(N)+2\sum_{d=1}^{\lfloor N/5\rfloor}\mu(d)\,A\!\left(\left\lfloor\frac{N}{d}\right\rfloor\right),$$
where \(A(V)\) is an unrestricted auxiliary counting function. The role of \(\mu(d)\) is to remove contributions coming from non-coprime multiples. This is why all three implementations build Möbius values and their prefix sums before the main accumulation loop.
To compute \(A(V)\), the implementation reverses one step of Euclid's algorithm. If a reduced tail pair is \((u,v)\) with \(u>v\ge 1\), then every predecessor that reaches it after one division has the shape
$$\bigl(ku+v,\;u\bigr),\qquad k\ge 1.$$
Imposing the upper bound \(ku+v\le V\) turns the counting problem into lattice-point counting under linear inequalities. After summing all admissible predecessors of all reduced tails, the required quantities reduce to floor sums of the form
$$\sum_t \left\lfloor\frac{\alpha t+\beta}{\gamma}\right\rfloor,$$
together with a second weighted variant of the same type. The implementation evaluates both of these by a recursive Euclidean-style transform, which is the key reason the method stays fast.
The predecessor lattice is symmetric between two coordinate regions. Because of that symmetry, the auxiliary computation only needs to enumerate reduced tail pairs with larger coordinate at most \(\lfloor\sqrt{V}\rfloor\). One accumulated sum covers the full asymmetric region, another covers the square overlap, and the final combination is a doubled total minus that overlap correction. This explains the integer-square-root bound seen in the code.
The Möbius sum is not evaluated one divisor at a time. The quotient
$$q=\left\lfloor\frac{N}{d}\right\rfloor$$
stays constant on intervals \(l\le d\le r\), where
$$r=\left\lfloor\frac{N}{q}\right\rfloor.$$
So an entire block contributes
$$\sum_{d=l}^{r}\mu(d)\,A(q)=\bigl(M(r)-M(l-1)\bigr)A(q).$$
Only \(O(\sqrt{N})\) distinct quotients occur, so the expensive auxiliary function is called much fewer times than in a linear scan.
For \(N=10\), we have \(m=5\), so the base part is
$$B(10)=(3\cdot 5-2)\cdot 5=65.$$
The residual term computed by the fast Möbius-floor-sum pipeline is \(18\), hence
$$L(10)=65+18=83.$$
Returning to the full ordered sum gives
$$S(10)=2L(10)+\frac{10\cdot 11}{2}=2\cdot 83+55=221,$$
which matches the small checkpoint used to validate the implementation.
The C++, Python, and Java implementations all follow the same plan. They first build the Möbius table and its prefix sums up to \(\lfloor N/5\rfloor\). They then evaluate the closed base term \(B(N)\), generate the harmonic intervals where \(\lfloor N/d\rfloor\) is constant, and for each distinct quotient compute the reverse-Euclid auxiliary kernel via recursive floor sums. Each block is weighted by the corresponding Möbius prefix difference, added to the base term, and finally converted back to the full ordered total through
$$S(N)=2L(N)+\frac{N(N+1)}{2}.$$
A naive method would inspect all \(N^2\) ordered pairs and run Euclid's algorithm separately, which is entirely infeasible for \(N=5\cdot 10^6\). The optimized method uses a linear-style Möbius sieve up to \(\lfloor N/5\rfloor\), only \(O(\sqrt{N})\) distinct quotient blocks in the outer summation, and a recursive floor-sum kernel whose depth is logarithmic in the relevant parameters. In practice this turns a quadratic enumeration into a heavily compressed arithmetic computation. The main memory cost is the Möbius and prefix tables, so the space usage is \(O(N/5)\).
Für positive ganze Zahlen \(x\) und \(y\) bezeichne \(E(x,y)\) die Anzahl der Divisionen des euklidischen Algorithmus, bis der Rest \(0\) wird. Gesucht ist
$$S(N)=\sum_{1\le x,y\le N} E(x,y)$$
für den großen Wert \(N=5\cdot 10^6\). Eine direkte Doppelschleife ist unmöglich: Es gibt \(N^2\) geordnete Paare, und für jedes Paar müsste der Algorithmus nochmals separat laufen. Deshalb formt die Lösung die Summe so um, dass der Großteil der Arbeit durch geschlossene Formeln, Möbius-Inversion und rekursive Floor-Summen erledigt wird.
Für \(x>y\) führt der Start mit \((y,x)\) zunächst nur zu einem Tausch; danach ist die Restkette dieselbe. Daher gilt
$$E(y,x)=E(x,y)+1 \qquad (x>y).$$
Außerdem ist \(E(x,x)=1\). Mit
$$L(N)=\sum_{1\le y<x\le N} E(x,y)$$
erhält man für die volle Summe
$$\begin{aligned} S(N) &=\sum_{x=1}^{N}E(x,x)+\sum_{1\le y<x\le N}\bigl(E(x,y)+E(y,x)\bigr) \\ &=N+\sum_{1\le y<x\le N}\bigl(2E(x,y)+1\bigr) \\ &=2L(N)+\frac{N(N+1)}{2}. \end{aligned}$$
Damit genügt es, die Summe unterhalb der Diagonalen effizient zu berechnen.
Jedes Paar mit \(x>y\) benötigt mindestens einen Schritt. Ein zweiter Schritt ist ebenfalls sofort erkennbar, wenn \(x<2y\): Dann ist der erste Quotient gleich \(1\), also bleibt nach einem Subtraktionsschritt ein nichttriviales Paar übrig. Somit
$$E(x,y)=1+\mathbf{1}_{x<2y}+R(x,y),\qquad R(x,y)\ge 0.$$
Der einfache Anteil ist also
$$B(N)=\sum_{1\le y<x\le N}\left(1+\mathbf{1}_{x<2y}\right).$$
Der erste Summand zählt schlicht alle Paare unterhalb der Diagonalen:
$$\sum_{1\le y<x\le N}1=\frac{N(N-1)}{2}.$$
Für den Indikator muss man die Paare mit \(y<x<2y\) zählen. Schreibt man \(N=2m\) oder \(N=2m+1\), so ergibt sich
$$\sum_{1\le y<x\le 2m}\mathbf{1}_{x<2y}=m(m-1),$$
$$\sum_{1\le y<x\le 2m+1}\mathbf{1}_{x<2y}=m^2.$$
Daraus folgt die im Programm verwendete geschlossene Formel
$$B(N)= \begin{cases} (3m-2)m, & N=2m, \\ 3m^2+m, & N=2m+1. \end{cases}$$
Der Restterm \(R(x,y)\) enthält genau die Schritte, die über diese universellen ersten ein oder zwei Schritte hinausgehen. Da der euklidische Algorithmus sich beim Skalieren beider Argumente mit demselben Faktor nicht ändert, gilt
$$R(dx,dy)=R(x,y)\qquad (d\ge 1).$$
Deshalb lässt sich der Restterm natürlich nach dem größten gemeinsamen Teiler gruppieren; entscheidend ist nur der primitive Kern des Paares.
Wichtig ist auch: \(R(x,y)>0\) tritt erst ab größeren primitiven Paaren auf. Die ersten Beispiele sind \((5,2)\) und \((5,3)\). Daher reicht die Möbius-Schicht nur bis \(\lfloor N/5\rfloor\).
Sei \(\mu\) die Möbius-Funktion und
$$M(t)=\sum_{d\le t}\mu(d)$$
ihre Präfixsumme. Mit Möbius-Inversion isoliert die schnelle Lösung die primitiven Paare und schreibt den Rest der unteren Dreieckssumme als
$$L(N)=B(N)+2\sum_{d=1}^{\lfloor N/5\rfloor}\mu(d)\,A\!\left(\left\lfloor\frac{N}{d}\right\rfloor\right),$$
wobei \(A(V)\) eine unbeschränkte Hilfszählfunktion ist. Die Möbius-Gewichte entfernen genau die Anteile, die nur von nichtteilerfremden Vielfachen stammen. Aus diesem Grund werden in allen drei Sprachen zuerst \(\mu(d)\) und die Präfixsummen aufgebaut.
Um \(A(V)\) auszuwerten, kehrt die Implementierung einen Euklid-Schritt um. Ist \((u,v)\) mit \(u>v\ge 1\) ein reduziertes Endstück, dann hat jeder Vorgänger, der nach einer Division zu diesem Paar führt, die Form
$$\bigl(ku+v,\;u\bigr),\qquad k\ge 1.$$
Die Schranke \(ku+v\le V\) verwandelt das Problem in das Zählen von Gitterpunkten unter linearen Ungleichungen. Nach Aufsummieren aller zulässigen Vorgänger aller reduzierten Endstücke bleiben Floor-Summen der Form
$$\sum_t \left\lfloor\frac{\alpha t+\beta}{\gamma}\right\rfloor$$
sowie eine zweite gewichtete Variante derselben Struktur. Genau diese Summen werden rekursiv mit Euklid-ähnlichen Transformationen berechnet.
Die zugrunde liegende Vorgänger-Geometrie ist zwischen zwei Koordinatenbereichen symmetrisch. Deshalb genügt es, reduzierte Endstücke nur bis zur ganzzahligen Wurzel \(\lfloor\sqrt{V}\rfloor\) zu durchlaufen. Eine akkumulierte Summe deckt den asymmetrischen Bereich ab, eine zweite den quadratischen Überlapp, und am Ende wird verdoppelt und der Überlapp subtrahiert. Daher erscheint im Code die Wurzelgrenze.
Die Möbius-Summe wird nicht für jedes \(d\) einzeln ausgewertet. Der Quotient
$$q=\left\lfloor\frac{N}{d}\right\rfloor$$
ist auf ganzen Intervallen \(l\le d\le r\) konstant, wobei
$$r=\left\lfloor\frac{N}{q}\right\rfloor.$$
Ein kompletter Block liefert also
$$\sum_{d=l}^{r}\mu(d)\,A(q)=\bigl(M(r)-M(l-1)\bigr)A(q).$$
Da nur \(O(\sqrt{N})\) verschiedene Quotienten auftreten, wird die teure Hilfsfunktion nur relativ selten aufgerufen.
Für \(N=10\) gilt \(m=5\), also
$$B(10)=(3\cdot 5-2)\cdot 5=65.$$
Die schnelle Restberechnung liefert \(18\), also
$$L(10)=65+18=83.$$
Damit folgt für die vollständige geordnete Summe
$$S(10)=2L(10)+\frac{10\cdot 11}{2}=2\cdot 83+55=221,$$
genau der kleine Prüfwert, den die Implementierungen ebenfalls reproduzieren.
Die C++-, Python- und Java-Implementierungen arbeiten nach demselben Schema. Zuerst werden die Möbius-Werte und ihre Präfixsummen bis \(\lfloor N/5\rfloor\) erzeugt. Danach wird die geschlossene Basiskomponente \(B(N)\) berechnet. Anschließend werden die harmonischen Intervalle mit konstantem \(\lfloor N/d\rfloor\) gebildet, und für jeden verschiedenen Quotienten wird der Rückwärts-Euklid-Kern über rekursive Floor-Summen ausgewertet. Jeder Block wird mit der passenden Möbius-Präfixdifferenz gewichtet, zum Basisterm addiert und am Ende über
$$S(N)=2L(N)+\frac{N(N+1)}{2}$$
wieder in die volle geordnete Summe zurückgeführt.
Eine naive Lösung würde alle \(N^2\) geordneten Paare untersuchen und für jedes Paar den euklidischen Algorithmus separat ausführen. Das ist für \(N=5\cdot 10^6\) ausgeschlossen. Die optimierte Methode verwendet ein lineares Möbius-Sieb bis \(\lfloor N/5\rfloor\), nur \(O(\sqrt{N})\) verschiedene Quotientenblöcke in der äußeren Summe und einen rekursiven Floor-Sum-Kern mit logarithmischer Tiefe. In der Praxis ersetzt diese Struktur eine quadratische Enumeration durch eine stark komprimierte arithmetische Berechnung. Der Speicherbedarf wird von Möbius- und Präfixtabellen dominiert und liegt deshalb bei \(O(N/5)\).
Pozitif tamsayılar \(x\) ve \(y\) için, \(E(x,y)\) Öklid algoritmasının kalan \(0\) olana kadar yaptığı bölme sayısı olsun. İstenen büyüklük
$$S(N)=\sum_{1\le x,y\le N} E(x,y)$$
ifadesidir ve burada \(N=5\cdot 10^6\) gibi çok büyük bir değerdir. Doğrudan çift döngü yaklaşımı uygulanamaz; çünkü \(N^2\) tane sıralı çift vardır ve her biri için algoritmayı ayrıca yürütmek gerekir. Çözüm bu yüzden toplamı, kapalı formüller, Möbius tersleme ve floor-sum özyinelemesiyle hesaplanabilecek bir yapıya dönüştürür.
\(x>y\) iken \((y,x)\) ile başlamak önce sadece bir yer değiştirme yapar, sonra aynı kalan zincirine girer. Bu nedenle
$$E(y,x)=E(x,y)+1 \qquad (x>y).$$
Ayrıca köşegen üzerinde \(E(x,x)=1\) olur. Şimdi
$$L(N)=\sum_{1\le y<x\le N} E(x,y)$$
tanımını yaparsak tüm sıralı toplam
$$\begin{aligned} S(N) &=\sum_{x=1}^{N}E(x,x)+\sum_{1\le y<x\le N}\bigl(E(x,y)+E(y,x)\bigr) \\ &=N+\sum_{1\le y<x\le N}\bigl(2E(x,y)+1\bigr) \\ &=2L(N)+\frac{N(N+1)}{2} \end{aligned}$$
şeklini alır. Böylece yalnızca alt üçgendeki toplamı hesaplamak yeterlidir.
\(x>y\) olan her çift en az bir Öklid adımı gerektirir. Ayrıca \(x<2y\) ise ilk bölüm katsayısı \(1\) olur; yani ilk çıkarımdan sonra elimizde hâlâ sonlu olmayan bir çift kalır ve ikinci adım da zorunludur. O halde
$$E(x,y)=1+\mathbf{1}_{x<2y}+R(x,y),\qquad R(x,y)\ge 0.$$
Kolay kısım
$$B(N)=\sum_{1\le y<x\le N}\left(1+\mathbf{1}_{x<2y}\right)$$
olur. İlk toplam sadece köşegen altındaki çift sayısıdır:
$$\sum_{1\le y<x\le N}1=\frac{N(N-1)}{2}.$$
İndikatör terimi için \(y<x<2y\) çiftlerini saymak gerekir. \(N=2m\) veya \(N=2m+1\) yazılırsa
$$\sum_{1\le y<x\le 2m}\mathbf{1}_{x<2y}=m(m-1),$$
$$\sum_{1\le y<x\le 2m+1}\mathbf{1}_{x<2y}=m^2$$
elde edilir. Sonuç olarak kodda kullanılan kapalı formül
$$B(N)= \begin{cases} (3m-2)m, & N=2m, \\ 3m^2+m, & N=2m+1 \end{cases}$$
şeklindedir.
\(R(x,y)\) terimi, bu zorunlu ilk bir veya iki adımdan sonra kalan katkıyı temsil eder. Öklid algoritması her iki argüman aynı katsayıyla çarpıldığında değişmediği için
$$R(dx,dy)=R(x,y)\qquad (d\ge 1)$$
olur. Bu yüzden artık katkı doğal olarak ortak bölen sınıflarına ayrılır; önemli olan sadece çiftin asal olmayan ortak kısmı çıkarıldıktan sonra kalan ilkel gövdesidir.
Bir başka önemli nokta şudur: \(R(x,y)>0\) olması için büyük koordinatın en az \(5\) olması gerekir. İlk örnekler \((5,2)\) ve \((5,3)\) çiftleridir. Bu da Möbius katmanının neden \(\lfloor N/5\rfloor\) ile sınırlı olduğunu açıklar.
\(\mu\) Möbius fonksiyonu olsun ve
$$M(t)=\sum_{d\le t}\mu(d)$$
ön-ek toplamını tanımlayalım. Hızlı yöntem, ilkel çiftleri Möbius tersleme ile ayırarak alt üçgen toplamının artık kısmını
$$L(N)=B(N)+2\sum_{d=1}^{\lfloor N/5\rfloor}\mu(d)\,A\!\left(\left\lfloor\frac{N}{d}\right\rfloor\right)$$
biçimine getirir. Burada \(A(V)\), tüm ortak bölenler ayrıştırılmadan önce kullanılan yardımcı sayım fonksiyonudur. \(\mu(d)\) katsayılarının görevi, ilkel olmayan çoklu katkıları tam olarak silmektir. Bu yüzden üç dildeki uygulamaların hepsi önce Möbius dizisini ve onun ön-ek toplamını kurar.
\(A(V)\) fonksiyonunu hesaplamak için çözüm bir Öklid adımını tersine çevirir. Eğer indirgenmiş kuyruk çifti \((u,v)\) ve \(u>v\ge 1\) ise, bir bölme sonra buna ulaşan her önceki çift
$$\bigl(ku+v,\;u\bigr),\qquad k\ge 1$$
şeklindedir. \(ku+v\le V\) üst sınırı bu problemi doğrusal eşitsizliklerin altındaki kafes noktalarını sayma problemine dönüştürür. Tüm kuyruk çiftleri üzerinden toplandığında ihtiyaç duyulan miktarlar
$$\sum_t \left\lfloor\frac{\alpha t+\beta}{\gamma}\right\rfloor$$
tipi floor-sum ifadelerine ve onların ağırlıklı bir varyantına indirgenir. Uygulama tam olarak bu çekirdeği özyinelemeli, Öklid-benzeri dönüşümlerle hesaplar.
Yardımcı sayımın arkasındaki önceki-adım kafesi iki koordinat bölgesi arasında simetriktir. Bu nedenle indirgenmiş kuyruk çiftlerini yalnızca \(\lfloor\sqrt{V}\rfloor\) sınırına kadar dolaşmak yeterlidir. Bir toplam asimetrik bölgeyi, ikinci toplam kare biçimli çakışmayı temsil eder; sonda birini ikiyle çarpıp diğerini çıkararak tam sonuç elde edilir. Koddaki karekök sınırı bundan gelir.
Möbius toplamı her \(d\) için tek tek hesaplanmaz. Çünkü
$$q=\left\lfloor\frac{N}{d}\right\rfloor$$
değeri bazı aralıklarda sabittir. Eğer \(l\le d\le r\) için sabit kalıyorsa
$$r=\left\lfloor\frac{N}{q}\right\rfloor$$
olur ve tam blok katkısı
$$\sum_{d=l}^{r}\mu(d)\,A(q)=\bigl(M(r)-M(l-1)\bigr)A(q)$$
şeklinde yazılır. Farklı bölüm değeri sayısı yalnızca \(O(\sqrt{N})\) olduğundan, pahalı yardımcı fonksiyon çok daha az çağrılır.
\(N=10\) için \(m=5\) olur ve taban katkısı
$$B(10)=(3\cdot 5-2)\cdot 5=65$$
çıkar. Hızlı artık hesaplama \(18\) verdiği için
$$L(10)=65+18=83$$
elde edilir. Sonra tam sıralı toplam
$$S(10)=2L(10)+\frac{10\cdot 11}{2}=2\cdot 83+55=221$$
olur; bu da küçük doğrulama değeriyle aynıdır.
C++, Python ve Java uygulamaları aynı akışı izler. Önce \(\lfloor N/5\rfloor\) değerine kadar Möbius tablosu ve ön-ek toplamı oluşturulur. Ardından kapalı taban terimi \(B(N)\) hesaplanır. Sonraki adımda \(\lfloor N/d\rfloor\) sabit kalan harmonik aralıklar kurulup her farklı bölüm değeri için ters-Öklid yardımcı çekirdeği özyinelemeli floor-sum formülleriyle değerlendirilir. Her blok uygun Möbius ön-ek farkıyla ağırlanır, taban katkısına eklenir ve en sonunda
$$S(N)=2L(N)+\frac{N(N+1)}{2}$$
bağıntısıyla tam sıralı toplama dönülür.
Naif bir yöntem \(N^2\) sıralı çifti dolaşıp her biri için Öklid algoritmasını ayrıca çalıştırır; bu, \(N=5\cdot 10^6\) için mümkün değildir. Optimizasyonlu yöntem \(\lfloor N/5\rfloor\) boyutuna kadar doğrusal tarzda bir Möbius eleği, dış toplamda yalnızca \(O(\sqrt{N})\) farklı bölüm bloğu ve parametrelere göre logaritmik derinlikte bir floor-sum çekirdeği kullanır. Pratikte bu yapı karesel taramayı sıkıştırılmış bir aritmetik toplama problemine dönüştürür. Bellek maliyeti esas olarak Möbius ve ön-ek dizileridir; bu yüzden alan karmaşıklığı \(O(N/5)\) düzeyindedir.
Para enteros positivos \(x\) e \(y\), sea \(E(x,y)\) el número de divisiones que ejecuta el algoritmo de Euclides hasta que el resto se hace \(0\). Hay que calcular
$$S(N)=\sum_{1\le x,y\le N} E(x,y)$$
para el valor muy grande \(N=5\cdot 10^6\). Un doble bucle directo no sirve: hay \(N^2\) pares ordenados y cada uno requeriría su propia ejecución del algoritmo. La solución reordena la suma para que casi todo el trabajo quede en fórmulas cerradas, inversión de Möbius y recurrencias de floor-sum.
Si \(x>y\), comenzar con \((y,x)\) hace primero un intercambio y luego sigue exactamente la misma cadena de restos. Por tanto,
$$E(y,x)=E(x,y)+1 \qquad (x>y).$$
Además, \(E(x,x)=1\) para todo \(x\). Definimos
$$L(N)=\sum_{1\le y<x\le N} E(x,y).$$
Entonces la suma total se convierte en
$$\begin{aligned} S(N) &=\sum_{x=1}^{N}E(x,x)+\sum_{1\le y<x\le N}\bigl(E(x,y)+E(y,x)\bigr) \\ &=N+\sum_{1\le y<x\le N}\bigl(2E(x,y)+1\bigr) \\ &=2L(N)+\frac{N(N+1)}{2}. \end{aligned}$$
Así, basta con calcular eficientemente la parte con \(x>y\).
Todo par con \(x>y\) aporta al menos un paso. También hay un segundo paso fácil de identificar: si \(x<2y\), el primer cociente es \(1\), así que después de una resta todavía queda un par no terminal. De este modo,
$$E(x,y)=1+\mathbf{1}_{x<2y}+R(x,y),\qquad R(x,y)\ge 0.$$
La contribución sencilla es
$$B(N)=\sum_{1\le y<x\le N}\left(1+\mathbf{1}_{x<2y}\right).$$
La primera suma cuenta simplemente los pares bajo la diagonal:
$$\sum_{1\le y<x\le N}1=\frac{N(N-1)}{2}.$$
Para el indicador hay que contar los pares con \(y<x<2y\). Si escribimos \(N=2m\) o \(N=2m+1\), obtenemos
$$\sum_{1\le y<x\le 2m}\mathbf{1}_{x<2y}=m(m-1),$$
$$\sum_{1\le y<x\le 2m+1}\mathbf{1}_{x<2y}=m^2.$$
Por lo tanto, el término base cerrado es
$$B(N)= \begin{cases} (3m-2)m, & N=2m, \\ 3m^2+m, & N=2m+1. \end{cases}$$
El término \(R(x,y)\) recoge todo lo que queda después de esos primeros uno o dos pasos universales. Como el algoritmo de Euclides no cambia al multiplicar ambos argumentos por el mismo factor,
$$R(dx,dy)=R(x,y)\qquad (d\ge 1).$$
Eso hace natural agrupar por el máximo común divisor: la parte difícil depende sólo del núcleo primitivo del par.
Además, \(R(x,y)>0\) no puede aparecer antes de que la coordenada mayor llegue a \(5\). Los primeros ejemplos primitivos son \((5,2)\) y \((5,3)\). Por eso la capa de Möbius sólo necesita llegar hasta \(\lfloor N/5\rfloor\).
Sea \(\mu\) la función de Möbius y
$$M(t)=\sum_{d\le t}\mu(d)$$
su suma prefija. El método rápido aísla los pares coprimos mediante inversión de Möbius y reescribe la parte residual del triángulo inferior como
$$L(N)=B(N)+2\sum_{d=1}^{\lfloor N/5\rfloor}\mu(d)\,A\!\left(\left\lfloor\frac{N}{d}\right\rfloor\right),$$
donde \(A(V)\) es una función auxiliar sin restricción de coprimalidad. Los pesos \(\mu(d)\) eliminan exactamente las contribuciones procedentes de múltiplos no primitivos. Por eso las implementaciones construyen primero \(\mu\) y sus prefijos.
Para evaluar \(A(V)\), la implementación invierte un paso del algoritmo de Euclides. Si \((u,v)\) con \(u>v\ge 1\) es una cola reducida, entonces todo predecesor que llega a ella tras una división tiene la forma
$$\bigl(ku+v,\;u\bigr),\qquad k\ge 1.$$
La restricción \(ku+v\le V\) transforma el problema en contar puntos de una red bajo desigualdades lineales. Al sumar todos los predecesores admisibles de todas las colas reducidas, aparecen sumas del tipo
$$\sum_t \left\lfloor\frac{\alpha t+\beta}{\gamma}\right\rfloor,$$
junto con una segunda variante ponderada. Ese es precisamente el núcleo matemático que las implementaciones evalúan recursivamente.
La geometría de predecesores es simétrica entre dos regiones de coordenadas. Gracias a ello, sólo hace falta enumerar colas reducidas cuya coordenada mayor sea como mucho \(\lfloor\sqrt{V}\rfloor\). Una suma acumulada cubre la región asimétrica completa, otra cubre la superposición cuadrada, y el resultado final se obtiene duplicando la primera y restando la corrección de solapamiento.
La suma con Möbius no recorre cada divisor por separado. El cociente
$$q=\left\lfloor\frac{N}{d}\right\rfloor$$
permanece constante en intervalos \(l\le d\le r\), donde
$$r=\left\lfloor\frac{N}{q}\right\rfloor.$$
Así, todo un bloque aporta
$$\sum_{d=l}^{r}\mu(d)\,A(q)=\bigl(M(r)-M(l-1)\bigr)A(q).$$
Como sólo aparecen \(O(\sqrt{N})\) cocientes distintos, la función auxiliar costosa se invoca muchas menos veces que en un recorrido lineal.
Para \(N=10\), tenemos \(m=5\), luego
$$B(10)=(3\cdot 5-2)\cdot 5=65.$$
La parte residual calculada por el método rápido vale \(18\), de modo que
$$L(10)=65+18=83.$$
Finalmente,
$$S(10)=2L(10)+\frac{10\cdot 11}{2}=2\cdot 83+55=221,$$
que coincide con el valor pequeño de validación.
Las implementaciones en C++, Python y Java siguen exactamente la misma estrategia. Primero construyen la tabla de Möbius y sus prefijos hasta \(\lfloor N/5\rfloor\). Después calculan el término base cerrado \(B(N)\), forman los bloques armónicos en los que \(\lfloor N/d\rfloor\) es constante y, para cada cociente distinto, evalúan el núcleo auxiliar basado en Euclides inverso mediante fórmulas recursivas de floor-sum. Cada bloque se pondera con la diferencia correspondiente de prefijos de Möbius, se suma al término base y, al final, se vuelve a la suma total ordenada usando
$$S(N)=2L(N)+\frac{N(N+1)}{2}.$$
Un método ingenuo examinaría los \(N^2\) pares ordenados y ejecutaría Euclides por separado en cada uno, algo imposible para \(N=5\cdot 10^6\). El método optimizado utiliza un cribado tipo lineal para \(\mu\) hasta \(\lfloor N/5\rfloor\), sólo \(O(\sqrt{N})\) bloques distintos en la suma exterior y un núcleo recursivo de floor-sum con profundidad logarítmica. En la práctica, esto reemplaza una enumeración cuadrática por un cálculo aritmético muy comprimido. La memoria está dominada por las tablas de Möbius y de prefijos, así que el espacio es \(O(N/5)\).
对正整数 \(x,y\)，记 \(E(x,y)\) 为欧几里得算法从 \((x,y)\) 开始直到余数变成 \(0\) 所执行的除法次数。题目要求计算
$$S(N)=\sum_{1\le x,y\le N} E(x,y)$$
其中 \(N=5\cdot 10^6\)。如果直接枚举，必须处理 \(N^2\) 个有序对，而且每个有序对还要单独跑一次欧几里得算法，这在规模上完全不可行。实现采用的办法是把总和拆成几个可以通过闭式、Möbius 反演和 floor-sum 递归来处理的部分。
当 \(x>y\) 时，从 \((y,x)\) 开始会先做一次交换，然后进入与 \((x,y)\) 完全相同的余数链，因此
$$E(y,x)=E(x,y)+1 \qquad (x>y).$$
对角线上有 \(E(x,x)=1\)。定义
$$L(N)=\sum_{1\le y<x\le N} E(x,y).$$
那么完整的有序和就是
$$\begin{aligned} S(N) &=\sum_{x=1}^{N}E(x,x)+\sum_{1\le y<x\le N}\bigl(E(x,y)+E(y,x)\bigr) \\ &=N+\sum_{1\le y<x\le N}\bigl(2E(x,y)+1\bigr) \\ &=2L(N)+\frac{N(N+1)}{2}. \end{aligned}$$
所以真正需要高效求的是下三角部分 \(L(N)\)。
对任意 \(x>y\)，欧几里得算法至少要走一步。如果 \(x<2y\)，那么第一次商一定是 \(1\)，做完这一步后还没有结束，所以第二步也是必然的。因此可以写成
$$E(x,y)=1+\mathbf{1}_{x<2y}+R(x,y),\qquad R(x,y)\ge 0.$$
容易计算的基底项是
$$B(N)=\sum_{1\le y<x\le N}\left(1+\mathbf{1}_{x<2y}\right).$$
第一部分只是对角线下方点的个数：
$$\sum_{1\le y<x\le N}1=\frac{N(N-1)}{2}.$$
指示函数部分需要统计满足 \(y<x<2y\) 的点。写 \(N=2m\) 或 \(N=2m+1\) 后可得
$$\sum_{1\le y<x\le 2m}\mathbf{1}_{x<2y}=m(m-1),$$
$$\sum_{1\le y<x\le 2m+1}\mathbf{1}_{x<2y}=m^2.$$
因此闭式基底项为
$$B(N)= \begin{cases} (3m-2)m, & N=2m, \\ 3m^2+m, & N=2m+1. \end{cases}$$
\(R(x,y)\) 表示在这前一到两步之外还剩下多少贡献。由于把两个参数同时乘以同一个正整数并不会改变欧几里得算法的步数，故有
$$R(dx,dy)=R(x,y)\qquad (d\ge 1).$$
这说明困难部分天然应该按 \(\gcd\) 分层；真正决定复杂度的是约去公因数之后得到的原始互素核。
还可以看出 \(R(x,y)>0\) 的最早原始例子是 \((5,2)\) 和 \((5,3)\)，也就是说只有当较大坐标至少为 \(5\) 时，剩余项才可能非零。这正是实现只需要处理到 \(\lfloor N/5\rfloor\) 的原因。
设 \(\mu\) 为 Möbius 函数，并定义前缀和
$$M(t)=\sum_{d\le t}\mu(d).$$
快速算法用 Möbius 反演把互素对单独提出来，并把下三角中的剩余部分写成
$$L(N)=B(N)+2\sum_{d=1}^{\lfloor N/5\rfloor}\mu(d)\,A\!\left(\left\lfloor\frac{N}{d}\right\rfloor\right),$$
其中 \(A(V)\) 是一个不要求互素的辅助计数函数。这里 \(\mu(d)\) 的作用是把所有非原始倍数层的贡献精确消掉。因此三种语言的实现都会先构造 Möbius 数组和其前缀和。
为了计算 \(A(V)\)，程序把欧几里得算法反过来看。若某个约化后的尾部对为 \((u,v)\)，其中 \(u>v\ge 1\)，那么所有在一次除法后落到它的前驱都形如
$$\bigl(ku+v,\;u\bigr),\qquad k\ge 1.$$
加上上界 \(ku+v\le V\) 后，问题变成线性不等式下的格点计数。把所有约化尾部对的合法前驱都加起来之后，所需的数量归结为
$$\sum_t \left\lfloor\frac{\alpha t+\beta}{\gamma}\right\rfloor$$
这类 floor-sum，以及一个与之配套的加权版本。实现正是用递归的欧几里得式变换同时求出这两个量。
辅助格点区域在两个坐标方向上具有对称性。因此只需要枚举较大坐标不超过 \(\lfloor\sqrt{V}\rfloor\) 的约化尾部对即可。一部分累加覆盖整个非对称区域，另一部分累加负责方形重叠区，最后通过“二倍主区域减去重叠修正”的方式恢复完整计数。这就是代码里出现整数平方根上界的原因。
Möbius 求和并不是对每个 \(d\) 单独做的。因为
$$q=\left\lfloor\frac{N}{d}\right\rfloor$$
会在区间 \(l\le d\le r\) 上保持不变，而
$$r=\left\lfloor\frac{N}{q}\right\rfloor.$$
于是整个区间的贡献就是
$$\sum_{d=l}^{r}\mu(d)\,A(q)=\bigl(M(r)-M(l-1)\bigr)A(q).$$
不同的商值只有 \(O(\sqrt{N})\) 个，所以昂贵的辅助函数调用次数远少于线性扫描。
当 \(N=10\) 时，\(m=5\)，因此基底项为
$$B(10)=(3\cdot 5-2)\cdot 5=65.$$
快速残差部分给出 \(18\)，所以
$$L(10)=65+18=83.$$
最后回到完整有序和：
$$S(10)=2L(10)+\frac{10\cdot 11}{2}=2\cdot 83+55=221,$$
这与实现中使用的小规模验证值一致。
C++、Python 和 Java 三个实现的结构完全一致。首先建立到 \(\lfloor N/5\rfloor\) 为止的 Möbius 表及其前缀和。然后直接计算闭式基底项 \(B(N)\)，再把所有使 \(\lfloor N/d\rfloor\) 相同的 \(d\) 合并成调和区间。对每个不同的商值，程序用递归 floor-sum 公式求出反向欧几里得辅助核，再乘上相应的 Möbius 前缀差并累加到基底项上。最后使用
$$S(N)=2L(N)+\frac{N(N+1)}{2}$$
把下三角总和还原为完整的有序对总和。
朴素方法需要检查 \(N^2\) 个有序对，并对每一对单独执行欧几里得算法；在 \(N=5\cdot 10^6\) 时这是不可能完成的。优化后的方法只需要做到 \(\lfloor N/5\rfloor\) 的 Möbius 线性筛，外层只有 \(O(\sqrt{N})\) 个不同的商值区间，而 floor-sum 递归核的深度对参数来说是对数级。实际运行中，这把一个二次枚举问题压缩成了一个高度结构化的算术求和问题。主要内存消耗来自 Möbius 与前缀数组，因此空间复杂度为 \(O(N/5)\)。
Пусть для положительных целых \(x\) и \(y\) величина \(E(x,y)\) обозначает число делений в алгоритме Евклида до первого нулевого остатка. Требуется вычислить
$$S(N)=\sum_{1\le x,y\le N} E(x,y)$$
для очень большого значения \(N=5\cdot 10^6\). Прямой перебор нереален: существует \(N^2\) упорядоченных пар, и для каждой пришлось бы отдельно запускать алгоритм. Поэтому решение перестраивает сумму так, чтобы основная работа свелась к замкнутым формулам, инверсии Мёбиуса и рекурсивным floor-sum вычислениям.
Если \(x>y\), то старт из \((y,x)\) сначала делает лишь перестановку аргументов, а затем проходит ту же цепочку остатков. Поэтому
$$E(y,x)=E(x,y)+1 \qquad (x>y).$$
На диагонали \(E(x,x)=1\). Введем
$$L(N)=\sum_{1\le y<x\le N} E(x,y).$$
Тогда полная сумма принимает вид
$$\begin{aligned} S(N) &=\sum_{x=1}^{N}E(x,x)+\sum_{1\le y<x\le N}\bigl(E(x,y)+E(y,x)\bigr) \\ &=N+\sum_{1\le y<x\le N}\bigl(2E(x,y)+1\bigr) \\ &=2L(N)+\frac{N(N+1)}{2}. \end{aligned}$$
Значит, достаточно быстро вычислить сумму ниже диагонали.
Каждая пара с \(x>y\) требует как минимум одного шага. Второй шаг тоже можно выделить сразу: если \(x<2y\), то первый частный равен \(1\), а после первого вычитания процесс еще не завершен. Значит,
$$E(x,y)=1+\mathbf{1}_{x<2y}+R(x,y),\qquad R(x,y)\ge 0.$$
Простая часть равна
$$B(N)=\sum_{1\le y<x\le N}\left(1+\mathbf{1}_{x<2y}\right).$$
Первое слагаемое просто считает пары ниже диагонали:
$$\sum_{1\le y<x\le N}1=\frac{N(N-1)}{2}.$$
Для индикатора нужно посчитать пары с \(y<x<2y\). Если \(N=2m\) или \(N=2m+1\), то
$$\sum_{1\le y<x\le 2m}\mathbf{1}_{x<2y}=m(m-1),$$
$$\sum_{1\le y<x\le 2m+1}\mathbf{1}_{x<2y}=m^2.$$
Следовательно, замкнутая базовая формула имеет вид
$$B(N)= \begin{cases} (3m-2)m, & N=2m, \\ 3m^2+m, & N=2m+1. \end{cases}$$
Термин \(R(x,y)\) содержит все шаги сверх этих универсальных первых одного или двух шагов. Поскольку алгоритм Евклида не меняется при умножении обоих аргументов на один и тот же множитель, имеем
$$R(dx,dy)=R(x,y)\qquad (d\ge 1).$$
Поэтому естественно группировать вклад по значениям общего делителя: трудная часть определяется только примитивным ядром пары.
Кроме того, \(R(x,y)>0\) впервые возникает лишь при старших координатах не меньше \(5\); первые примитивные примеры это \((5,2)\) и \((5,3)\). Именно поэтому слой с функцией Мёбиуса ограничивается значением \(\lfloor N/5\rfloor\).
Пусть \(\mu\) — функция Мёбиуса, а
$$M(t)=\sum_{d\le t}\mu(d)$$
— ее префиксная сумма. Быстрый метод выделяет взаимно простые пары с помощью инверсии Мёбиуса и записывает остаточную часть нижнего треугольника так:
$$L(N)=B(N)+2\sum_{d=1}^{\lfloor N/5\rfloor}\mu(d)\,A\!\left(\left\lfloor\frac{N}{d}\right\rfloor\right),$$
где \(A(V)\) — вспомогательная функция без ограничения взаимной простоты. Коэффициенты \(\mu(d)\) убирают вклад непервичных кратных, поэтому во всех реализациях сначала строятся значения функции Мёбиуса и их префиксы.
Чтобы вычислить \(A(V)\), реализация обращает один шаг алгоритма Евклида. Если редуцированный хвост равен \((u,v)\), где \(u>v\ge 1\), то каждый предок, который после одного деления переходит в него, имеет вид
$$\bigl(ku+v,\;u\bigr),\qquad k\ge 1.$$
Ограничение \(ku+v\le V\) превращает задачу в подсчет решеточных точек под линейными границами. После суммирования по всем допустимым предкам всех редуцированных хвостов возникают floor-sum выражения типа
$$\sum_t \left\lfloor\frac{\alpha t+\beta}{\gamma}\right\rfloor,$$
а также вторая взвешенная сумма такого же типа. Именно эти величины вычисляются рекурсивным евклидовым преобразованием.
Геометрия области предков симметрична относительно двух координатных областей. Поэтому достаточно перечислять редуцированные хвосты лишь до \(\lfloor\sqrt{V}\rfloor\). Одна накопленная сумма покрывает всю несимметричную часть, вторая отвечает за квадратное пересечение, а в конце берется удвоение первой минус поправка на пересечение. Отсюда и граница по целой квадратной корню в коде.
Сумма по Мёбиусу вычисляется не по одному делителю. Частное
$$q=\left\lfloor\frac{N}{d}\right\rfloor$$
постоянно на интервале \(l\le d\le r\), где
$$r=\left\lfloor\frac{N}{q}\right\rfloor.$$
Значит, вклад целого блока равен
$$\sum_{d=l}^{r}\mu(d)\,A(q)=\bigl(M(r)-M(l-1)\bigr)A(q).$$
Различных значений частного всего \(O(\sqrt{N})\), поэтому дорогая вспомогательная функция вызывается сравнительно редко.
При \(N=10\) имеем \(m=5\), и базовая часть равна
$$B(10)=(3\cdot 5-2)\cdot 5=65.$$
Быстрая остаточная часть дает \(18\), следовательно,
$$L(10)=65+18=83.$$
После возврата к полной упорядоченной сумме получаем
$$S(10)=2L(10)+\frac{10\cdot 11}{2}=2\cdot 83+55=221,$$
что совпадает с малым контрольным значением, используемым для проверки реализации.
Реализации на C++, Python и Java устроены одинаково. Сначала строится таблица значений функции Мёбиуса и ее префиксных сумм до \(\lfloor N/5\rfloor\). Затем сразу вычисляется замкнутая базовая часть \(B(N)\). После этого формируются гармонические интервалы, на которых \(\lfloor N/d\rfloor\) постоянно, и для каждого различного частного вычисляется вспомогательное ядро обратного алгоритма Евклида через рекурсивные floor-sum формулы. Каждый блок умножается на соответствующую разность префиксов Мёбиуса, добавляется к базе, а в конце нижнетреугольная сумма переводится обратно в полную через
$$S(N)=2L(N)+\frac{N(N+1)}{2}.$$
Наивный подход проверял бы все \(N^2\) упорядоченных пар и отдельно запускал алгоритм Евклида для каждой из них, что невозможно при \(N=5\cdot 10^6\). Оптимизированный метод использует линейное по стилю решето для функции Мёбиуса до \(\lfloor N/5\rfloor\), только \(O(\sqrt{N})\) различных блоков частных во внешней сумме и рекурсивное floor-sum ядро логарифмической глубины. На практике это заменяет квадратичный перебор на сильно сжатое арифметическое вычисление. Основная память уходит на таблицы Мёбиуса и префиксов, поэтому пространственная сложность равна \(O(N/5)\).
لأي عددين صحيحين موجبين \(x\) و\(y\)، نرمز بـ \(E(x,y)\) إلى عدد القسَمات التي تنفذها خوارزمية إقليدس حتى يصبح الباقي \(0\). المطلوب هو حساب
$$S(N)=\sum_{1\le x,y\le N} E(x,y)$$
عند القيمة الكبيرة جدًا \(N=5\cdot 10^6\). الحل المباشر غير ممكن عمليًا، لأن لدينا \(N^2\) زوجًا مرتبًا، ولكل زوج يجب تشغيل الخوارزمية على حدة. لذلك يعيد الحل تنظيم المجموع بحيث تتحول أغلب العملية إلى صيغ مغلقة، وقلب Möbius، ونواة floor-sum سريعة ذات بنية تكرارية.
إذا كان \(x>y\)، فإن البدء من \((y,x)\) يضيف فقط خطوة تبادل أولية، ثم يدخل في سلسلة البواقي نفسها تمامًا. لذلك
$$E(y,x)=E(x,y)+1 \qquad (x>y).$$
وعلى القطر نحصل دائمًا على \(E(x,x)=1\). إذا عرّفنا
$$L(N)=\sum_{1\le y<x\le N} E(x,y),$$
فإن المجموع الكامل يصبح
$$\begin{aligned} S(N) &=\sum_{x=1}^{N}E(x,x)+\sum_{1\le y<x\le N}\bigl(E(x,y)+E(y,x)\bigr) \\ &=N+\sum_{1\le y<x\le N}\bigl(2E(x,y)+1\bigr) \\ &=2L(N)+\frac{N(N+1)}{2}. \end{aligned}$$
إذن يكفي أن نحسب مجموع الجزء الذي فيه \(x>y\).
كل زوج يحقق \(x>y\) يساهم بخطوة واحدة على الأقل. وهناك خطوة ثانية يمكن التعرف عليها مباشرة: إذا كان \(x<2y\)، فإن خارج القسمة الأول يساوي \(1\)، وبالتالي تبقى بعد الخطوة الأولى حالة غير منتهية، فتكون الخطوة الثانية مضمونة أيضًا. لذلك
$$E(x,y)=1+\mathbf{1}_{x<2y}+R(x,y),\qquad R(x,y)\ge 0.$$
والجزء السهل يساوي
$$B(N)=\sum_{1\le y<x\le N}\left(1+\mathbf{1}_{x<2y}\right).$$
الحد الأول هو ببساطة عدد النقاط تحت القطر:
$$\sum_{1\le y<x\le N}1=\frac{N(N-1)}{2}.$$
أما الحد المؤشري فيتطلب عدّ الأزواج التي تحقق \(y<x<2y\). إذا كتبنا \(N=2m\) أو \(N=2m+1\)، نحصل على
$$\sum_{1\le y<x\le 2m}\mathbf{1}_{x<2y}=m(m-1),$$
$$\sum_{1\le y<x\le 2m+1}\mathbf{1}_{x<2y}=m^2.$$
ومن ثم تكون الصيغة المغلقة للحد الأساسي
$$B(N)= \begin{cases} (3m-2)m, & N=2m, \\ 3m^2+m, & N=2m+1. \end{cases}$$
الحد \(R(x,y)\) يجمع كل ما يأتي بعد تلك الخطوات الأولى العامة. وبما أن خوارزمية إقليدس لا تتغير إذا ضربنا المدخلين في العامل نفسه، فإن
$$R(dx,dy)=R(x,y)\qquad (d\ge 1).$$
لهذا السبب يصبح تنظيم الأزواج بحسب القاسم المشترك الأكبر أمرًا طبيعيًا؛ فالجزء الصعب يعتمد فقط على اللبّ الأولي للزوج بعد إزالة العامل المشترك.
كما أن \(R(x,y)>0\) لا يظهر قبل أن تصل الإحداثية الكبرى إلى \(5\). أول مثالين أوليين هما \((5,2)\) و\((5,3)\). وهذا يفسر لماذا يكفي المرور في طبقة Möbius حتى \(\lfloor N/5\rfloor\).
لتكن \(\mu\) هي دالة Möbius، ولتكن
$$M(t)=\sum_{d\le t}\mu(d)$$
هي مجاميعها السابقة. الطريقة السريعة تعزل الأزواج المتباينة أوليًا باستخدام قلب Möbius، وتكتب الجزء المتبقي من مجموع المثلث السفلي على الصورة
$$L(N)=B(N)+2\sum_{d=1}^{\lfloor N/5\rfloor}\mu(d)\,A\!\left(\left\lfloor\frac{N}{d}\right\rfloor\right).$$
هنا \(A(V)\) دالة مساعدة لا تفرض شرط التباين الأولي. وظيفة الأوزان \(\mu(d)\) هي إزالة مساهمات المضاعفات غير الأولية بدقة، ولهذا تبدأ جميع التطبيقات ببناء قيم Möbius ومجاميعها السابقة.
لحساب \(A(V)\)، تقلب الخوارزمية خطوة واحدة من خوارزمية إقليدس. إذا كانت النهاية المختزلة هي \((u,v)\) مع \(u>v\ge 1\)، فإن كل سلف يصل إليها بعد قسمة واحدة يأخذ الشكل
$$\bigl(ku+v,\;u\bigr),\qquad k\ge 1.$$
الشرط \(ku+v\le V\) يحول المسألة إلى عدّ نقاط شبكية تحت متراجحات خطية. وبعد جمع جميع الأسلاف المسموح بها لكل النهايات المختزلة، نصل إلى مجاميع من النوع
$$\sum_t \left\lfloor\frac{\alpha t+\beta}{\gamma}\right\rfloor,$$
مع نسخة موزونة ثانية من البنية نفسها. هذه هي النواة الرياضية التي تحسبها التطبيقات بطريقة سريعة تعتمد على التكرار.
منطقة العد في الفضاء الشبكي متماثلة بين جهتين من الإحداثيات. لذلك يكفي تعداد النهايات المختزلة التي لا تتجاوز إحداثيتها الكبرى \(\lfloor\sqrt{V}\rfloor\). أحد المجاميع يغطي المنطقة غير المتماثلة، ومجموع ثانٍ يعالج منطقة التداخل المربعة، ثم يستعاد الجواب الكامل عن طريق مضاعفة الأول وطرح تصحيح التداخل. لهذا يظهر حد الجذر التربيعي في التنفيذ.
مجموع Möbius لا يُحسب حدًا حدًا. فالقيمة
$$q=\left\lfloor\frac{N}{d}\right\rfloor$$
تبقى ثابتة على فترات من الشكل \(l\le d\le r\)، حيث
$$r=\left\lfloor\frac{N}{q}\right\rfloor.$$
وعليه تكون مساهمة الكتلة كاملة
$$\sum_{d=l}^{r}\mu(d)\,A(q)=\bigl(M(r)-M(l-1)\bigr)A(q).$$
وبما أن عدد القيم المختلفة لهذا الخارج هو فقط \(O(\sqrt{N})\)، فإن عدد استدعاءات الدالة المساعدة المكلفة ينخفض كثيرًا.
عند \(N=10\) يكون \(m=5\)، ومن ثم
$$B(10)=(3\cdot 5-2)\cdot 5=65.$$
ويعطي الجزء المتبقي المحسوب سريعًا القيمة \(18\)، ولذلك
$$L(10)=65+18=83.$$
ثم نحصل على المجموع الكامل
$$S(10)=2L(10)+\frac{10\cdot 11}{2}=2\cdot 83+55=221,$$
وهو نفس مقدار التحقق الصغير المستخدم في التنفيذ.
تتبع تطبيقات C++ وPython وJava الخطة نفسها. أولًا تُبنى دالة Möbius ومجاميعها السابقة حتى \(\lfloor N/5\rfloor\). بعد ذلك يُحسب الحد الأساسي المغلق \(B(N)\)، ثم تُجمع القيم في كتل يكون فيها \(\lfloor N/d\rfloor\) ثابتًا. ولكل خارج مختلف تُحسب نواة المساعدة المبنية على عكس خطوات إقليدس بواسطة صيغ floor-sum سريعة. تُوزن كل كتلة بفارق المجاميع السابقة المناسب، وتُضاف إلى الحد الأساسي، ثم يُعاد بناء المجموع الكامل باستخدام
$$S(N)=2L(N)+\frac{N(N+1)}{2}.$$
الطريقة الساذجة تفحص جميع الأزواج المرتبة وعددها \(N^2\)، وتشغّل خوارزمية إقليدس لكل زوج على حدة، وهذا غير ممكن عند \(N=5\cdot 10^6\). أما الطريقة المحسّنة فتستخدم غربال Möbius حتى \(\lfloor N/5\rfloor\)، وعددًا لا يتجاوز \(O(\sqrt{N})\) من كتل القيم المختلفة في المجموع الخارجي، ونواة floor-sum ذات عمق لوغاريتمي. عمليًا يتحول التعداد التربيعي إلى حساب حسابي مضغوط جدًا. الذاكرة تهيمن عليها جداول Möbius والمجاميع السابقة، ولذلك يكون استهلاك الذاكرة من رتبة \(O(N/5)\).