For each positive integer \(n\), define
$$f(n)=\left|\left\{(x,y)\in\mathbb{Z}_{>0}^2 : x\le y,\ \operatorname{lcm}(x,y)=n\right\}\right|,$$
and then define the summatory function
$$g(N)=\sum_{n=1}^{N} f(n).$$
The Project Euler task asks for \(g(10^{12})\). A direct search over pairs \((x,y)\) is hopeless, so the solution first derives a closed form for \(f(n)\), then transforms the summatory problem into a Möbius-weighted sum involving the ternary divisor function.
Write
$$n=\prod_{p} p^{a_p}.$$
If \(\operatorname{lcm}(x,y)=n\), then both \(x\) and \(y\) divide \(n\). For each prime \(p^{a_p}\parallel n\), write the \(p\)-adic exponents of \(x\) and \(y\) as \(\alpha_p\) and \(\beta_p\). They must satisfy
$$0\le \alpha_p,\beta_p\le a_p,\qquad \max(\alpha_p,\beta_p)=a_p.$$
For one fixed prime power \(p^{a}\), the number of ordered exponent pairs \((\alpha,\beta)\) is
$$\underbrace{(a+1)^2}_{0\le \alpha,\beta\le a}-\underbrace{a^2}_{0\le \alpha,\beta\le a-1}=2a+1.$$
Different primes act independently, so the total number of ordered pairs \((x,y)\) with \(\operatorname{lcm}(x,y)=n\) is
$$\prod_{p\mid n}(2a_p+1)=\tau(n^2),$$
because \(n^2=\prod p^{2a_p}\) and therefore \(\tau(n^2)=\prod (2a_p+1)\).
The problem does not count ordered pairs; it counts only pairs with \(x\le y\). Every non-diagonal ordered pair appears twice, once as \((x,y)\) and once as \((y,x)\). The only diagonal pair with least common multiple equal to \(n\) is \((n,n)\). Hence
$$\boxed{f(n)=\frac{\tau(n^2)+1}{2}.}$$
Summing the closed form gives
$$g(N)=\sum_{n\le N}\frac{\tau(n^2)+1}{2}=\frac{1}{2}\left(N+\sum_{n\le N}\tau(n^2)\right).$$
So the core task is to evaluate
$$S(N)=\sum_{n\le N}\tau(n^2).$$
Once \(S(N)\) is known, the required value is simply
$$g(N)=\frac{S(N)+N}{2}.$$
Let \(d_3(m)\) denote the ternary divisor function, i.e. the number of ordered factorizations \(m=abc\). Its Dirichlet series is
$$\sum_{m\ge 1}\frac{d_3(m)}{m^s}=\zeta(s)^3.$$
For the square-divisor count we have the classical identity
$$\sum_{n\ge 1}\frac{\tau(n^2)}{n^s}=\frac{\zeta(s)^3}{\zeta(2s)}.$$
Using
$$\frac{1}{\zeta(2s)}=\sum_{d\ge 1}\frac{\mu(d)}{d^{2s}},$$
and comparing coefficients, we obtain
$$\tau(n^2)=\sum_{d^2m=n}\mu(d)\,d_3(m)=\sum_{d^2\mid n}\mu(d)\,d_3\!\left(\frac{n}{d^2}\right).$$
Now sum over \(n\le N\) and interchange the order:
$$S(N)=\sum_{d\le \sqrt N}\mu(d)\,D_3\!\left(\left\lfloor\frac{N}{d^2}\right\rfloor\right),$$
where
$$D_3(X)=\sum_{m\le X} d_3(m)=\left|\left\{(a,b,c)\in\mathbb{Z}_{>0}^3: abc\le X\right\}\right|.$$
This identity is the main bridge from the arithmetic function \(\tau(n^2)\) to something the program can compute quickly.
The helper function summatory_d3 computes \(D_3(X)\). Instead of counting all triples directly, it sorts them by their smallest factor. Let
$$c=\left\lfloor X^{1/3}\right\rfloor,\qquad r_z=\left\lfloor\sqrt{\frac{X}{z}}\right\rfloor.$$
After separating the cases “all three equal”, “exactly two equal”, and “all distinct”, the ordered triple count becomes
$$D_3(X)=3\sum_{z=1}^{c}\left(2\sum_{x=z+1}^{r_z}\left\lfloor\frac{X}{zx}\right\rfloor-r_z^2+\left\lfloor\frac{X}{z^2}\right\rfloor\right)+c^3.$$
The final term \(c^3\) is the compact form of the telescoping contribution coming from triples with repeated coordinates. The inner floor-sum
$$\sum_{x=x_1}^{x_2}\left\lfloor\frac{n}{x}\right\rfloor$$
is exactly what the C++ and Java helper sq_sum_floor evaluates. It uses quotient-difference updates rather than recomputing every division from scratch, which is why the routine stays fast even when called many times.
A direct evaluation of
$$S(N)=\sum_{d\le \sqrt N}\mu(d)\,D_3\!\left(\left\lfloor\frac{N}{d^2}\right\rfloor\right)$$
still runs over too many values of \(d\). The implementation uses the two-scale choice
$$I=\left\lfloor N^{1/3}\right\rfloor,\qquad D=\left\lfloor\sqrt{\frac{N}{I}}\right\rfloor.$$
For \(d\le D\), the code evaluates the terms directly. For \(d>D\), the quantity \(\left\lfloor N/d^2\right\rfloor\) is smaller than \(I\), so many consecutive \(d\)-values share the same small argument. Introduce the Mertens function
$$M(x)=\sum_{n\le x}\mu(n),$$
and the first differences
$$\Delta D_3(i)=D_3(i)-D_3(i-1)=d_3(i).$$
Then the large-\(d\) tail can be regrouped as
$$\sum_{d>D}\mu(d)\,D_3\!\left(\left\lfloor\frac{N}{d^2}\right\rfloor\right)=\sum_{i=1}^{I-1}\Delta D_3(i)\left(M\!\left(\left\lfloor\sqrt{\frac{N}{i}}\right\rfloor\right)-M(D)\right).$$
This is exactly what the arrays small_diff, mertens_small, and mertens_big implement. The subtraction by M(D) * D3(I-1) in the code is the algebraic device that turns prefix values of \(D_3\) into the needed differences \(\Delta D_3(i)\).
The program does not sieve \(\mu(d)\) all the way to \(\sqrt N\). It only sieves up to \(D\) and computes larger Mertens values \(M(v)\) on demand for the arguments \(v=\left\lfloor\sqrt{N/i}\right\rfloor\). The recurrence comes from the identity
$$\sum_{d\le v}\mu(d)\left\lfloor\frac{v}{d}\right\rfloor=1.$$
After separating the range at \(\lfloor\sqrt v\rfloor\), one obtains the form used in the code:
$$M(v)=1-v+\lfloor\sqrt v\rfloor\,M(\lfloor\sqrt v\rfloor)-\sum_{d=2}^{\lfloor\sqrt v\rfloor} M\!\left(\left\lfloor\frac{v}{d}\right\rfloor\right)-\sum_{d=2}^{\lfloor\sqrt v\rfloor}\mu(d)\left\lfloor\frac{v}{d}\right\rfloor.$$
Because the needed arguments are highly structured, storing them in mertens_big[i] is enough; no full sieve up to \(\sqrt N\) is necessary.
Take \(n=12=2^2\cdot 3\). Then
$$\tau(12^2)=\tau(2^4\cdot 3^2)=(4+1)(2+1)=15,$$
so
$$f(12)=\frac{15+1}{2}=8.$$
The checkpoint values embedded in the C++ program are
$$g(10)=29,\qquad g(100)=647,\qquad g(1000)=11751,\qquad g(10^4)=186991,\qquad g(10^6)=37429395,$$
and the optimized method matches all of them before attempting \(g(10^{12})\).
The C++ source is the authoritative implementation. It provides exact integer square-root and cube-root helpers, a fast linear Möbius sieve up to \(D\), the floor-sum helper sq_sum_floor, and the triple-count routine summatory_d3. The main function solve computes the direct small-\(d\) contribution, then adds the grouped tail through Mertens values, finally using
$$g(N)=\frac{S(N)+N}{2}.$$
The functions brute_f_from_lcm and formula_f_from_factorization are used only for checkpoint verification on small inputs. The Java file is a direct translation of the optimized algorithm. The Python file is intentionally different: it is a lightweight bridge that compiles and runs the C++ solver and then parses the printed answer.
The memory footprint is dominated by the Möbius table up to \(D\) and the grouped Mertens arrays up to \(I\), so it is about \(O(N^{1/3})\). The total runtime is far below naive enumeration and, in the implementation-oriented estimate used for this approach, is roughly \(O(N^{2/3})\) arithmetic with fast floor-sum subroutines. That is what makes \(N=10^{12}\) feasible.
Für jede positive ganze Zahl \(n\) sei
$$f(n)=\left|\left\{(x,y)\in\mathbb{Z}_{>0}^2 : x\le y,\ \operatorname{lcm}(x,y)=n\right\}\right|,$$
und dazu die Summenfunktion
$$g(N)=\sum_{n=1}^{N} f(n).$$
Gesucht ist also \(g(10^{12})\). Ein direktes Durchprobieren aller Paare \((x,y)\) ist ausgeschlossen. Die Lösung leitet deshalb zuerst eine geschlossene Formel für \(f(n)\) her und verwandelt die Summation danach in eine Möbius-gewichtete Auswertung der ternären Teilerfunktion.
Schreibe
$$n=\prod_{p} p^{a_p}.$$
Falls \(\operatorname{lcm}(x,y)=n\), dann sind \(x\) und \(y\) Teiler von \(n\). Für jede Primzahl \(p^{a_p}\parallel n\) bezeichnen \(\alpha_p\) und \(\beta_p\) die \(p\)-adischen Exponenten von \(x\) bzw. \(y\). Dann gilt
$$0\le \alpha_p,\beta_p\le a_p,\qquad \max(\alpha_p,\beta_p)=a_p.$$
Für eine feste Primzahlpotenz \(p^a\) ist die Anzahl geordneter Exponentenpaare
$$\underbrace{(a+1)^2}_{0\le \alpha,\beta\le a}-\underbrace{a^2}_{0\le \alpha,\beta\le a-1}=2a+1.$$
Da sich die Primfaktoren unabhängig kombinieren, ist die Gesamtzahl geordneter Paare \((x,y)\) mit \(\operatorname{lcm}(x,y)=n\)
$$\prod_{p\mid n}(2a_p+1)=\tau(n^2),$$
denn aus \(n^2=\prod p^{2a_p}\) folgt \(\tau(n^2)=\prod (2a_p+1)\).
Die Aufgabe zählt jedoch keine geordneten Paare, sondern nur solche mit \(x\le y\). Jedes nichtdiagonale geordnete Paar erscheint genau zweimal, als \((x,y)\) und \((y,x)\). Das einzige diagonale Paar mit kgV gleich \(n\) ist \((n,n)\). Daher
$$\boxed{f(n)=\frac{\tau(n^2)+1}{2}.}$$
Durch Aufsummieren erhält man
$$g(N)=\sum_{n\le N}\frac{\tau(n^2)+1}{2}=\frac{1}{2}\left(N+\sum_{n\le N}\tau(n^2)\right).$$
Damit reduziert sich das Problem auf
$$S(N)=\sum_{n\le N}\tau(n^2).$$
Ist \(S(N)\) bekannt, so folgt unmittelbar
$$g(N)=\frac{S(N)+N}{2}.$$
Sei \(d_3(m)\) die ternäre Teilerfunktion, also die Anzahl geordneter Faktorisierungen \(m=abc\). Ihre Dirichlet-Reihe lautet
$$\sum_{m\ge 1}\frac{d_3(m)}{m^s}=\zeta(s)^3.$$
Für \(\tau(n^2)\) gilt die klassische Identität
$$\sum_{n\ge 1}\frac{\tau(n^2)}{n^s}=\frac{\zeta(s)^3}{\zeta(2s)}.$$
Mit
$$\frac{1}{\zeta(2s)}=\sum_{d\ge 1}\frac{\mu(d)}{d^{2s}}$$
liefert Koeffizientenvergleich
$$\tau(n^2)=\sum_{d^2m=n}\mu(d)\,d_3(m)=\sum_{d^2\mid n}\mu(d)\,d_3\!\left(\frac{n}{d^2}\right).$$
Summiert man über \(n\le N\) und vertauscht die Reihenfolge, erhält man
$$S(N)=\sum_{d\le \sqrt N}\mu(d)\,D_3\!\left(\left\lfloor\frac{N}{d^2}\right\rfloor\right),$$
wobei
$$D_3(X)=\sum_{m\le X} d_3(m)=\left|\left\{(a,b,c)\in\mathbb{Z}_{>0}^3: abc\le X\right\}\right|.$$
Genau diese Formel wird von der Implementierung ausgewertet.
Die Hilfsfunktion summatory_d3 berechnet \(D_3(X)\). Dazu werden die Tripel nicht direkt durchlaufen, sondern nach ihrem kleinsten Faktor geordnet. Setze
$$c=\left\lfloor X^{1/3}\right\rfloor,\qquad r_z=\left\lfloor\sqrt{\frac{X}{z}}\right\rfloor.$$
Nach der Trennung der Fälle „alle drei gleich“, „genau zwei gleich“ und „alle verschieden“ ergibt sich
$$D_3(X)=3\sum_{z=1}^{c}\left(2\sum_{x=z+1}^{r_z}\left\lfloor\frac{X}{zx}\right\rfloor-r_z^2+\left\lfloor\frac{X}{z^2}\right\rfloor\right)+c^3.$$
Der Term \(c^3\) ist die kompakte teleskopische Restform der Beiträge mit mehrfach gleichen Koordinaten. Die innere Floorsumme
$$\sum_{x=x_1}^{x_2}\left\lfloor\frac{n}{x}\right\rfloor$$
wird in C++ und Java von sq_sum_floor ausgewertet. Dort werden Quotientendifferenzen fortgeschrieben, statt jede Division neu zu berechnen.
Auch
$$S(N)=\sum_{d\le \sqrt N}\mu(d)\,D_3\!\left(\left\lfloor\frac{N}{d^2}\right\rfloor\right)$$
ist noch zu groß für eine naive Auswertung. Deshalb benutzt der Code die Zweiskalenwahl
$$I=\left\lfloor N^{1/3}\right\rfloor,\qquad D=\left\lfloor\sqrt{\frac{N}{I}}\right\rfloor.$$
Für \(d\le D\) werden die Terme direkt summiert. Für \(d>D\) ist \(\left\lfloor N/d^2\right\rfloor<I\), also treten nur noch kleine Argumente auf, die sich gruppieren lassen. Mit der Mertens-Funktion
$$M(x)=\sum_{n\le x}\mu(n)$$
und den Differenzen
$$\Delta D_3(i)=D_3(i)-D_3(i-1)=d_3(i)$$
wird der Restblock zu
$$\sum_{d>D}\mu(d)\,D_3\!\left(\left\lfloor\frac{N}{d^2}\right\rfloor\right)=\sum_{i=1}^{I-1}\Delta D_3(i)\left(M\!\left(\left\lfloor\sqrt{\frac{N}{i}}\right\rfloor\right)-M(D)\right).$$
Genau dafür stehen im Programm die Arrays small_diff, mertens_small und mertens_big. Die Subtraktion M(D) * D3(I-1) wandelt Präfixsummen von \(D_3\) algebraisch in die benötigten Differenzen um.
Bis \(\sqrt N\) wird \(\mu\) nicht vollständig gesiebt. Stattdessen siebt das Programm nur bis \(D\) und berechnet größere Werte \(M(v)\) für \(v=\left\lfloor\sqrt{N/i}\right\rfloor\) rekursiv. Ausgangspunkt ist
$$\sum_{d\le v}\mu(d)\left\lfloor\frac{v}{d}\right\rfloor=1.$$
Nach Aufteilung bei \(\lfloor\sqrt v\rfloor\) ergibt sich die im Code verwendete Form
$$M(v)=1-v+\lfloor\sqrt v\rfloor\,M(\lfloor\sqrt v\rfloor)-\sum_{d=2}^{\lfloor\sqrt v\rfloor} M\!\left(\left\lfloor\frac{v}{d}\right\rfloor\right)-\sum_{d=2}^{\lfloor\sqrt v\rfloor}\mu(d)\left\lfloor\frac{v}{d}\right\rfloor.$$
Da nur die speziellen Argumente \(\left\lfloor\sqrt{N/i}\right\rfloor\) gebraucht werden, reicht das zwischengespeicherte Array mertens_big aus.
Nehmen wir \(n=12=2^2\cdot 3\). Dann ist
$$\tau(12^2)=\tau(2^4\cdot 3^2)=(4+1)(2+1)=15,$$
also
$$f(12)=\frac{15+1}{2}=8.$$
Die in der C++-Datei hinterlegten Prüfpunkte lauten
$$g(10)=29,\qquad g(100)=647,\qquad g(1000)=11751,\qquad g(10^4)=186991,\qquad g(10^6)=37429395,$$
und die optimierte Methode reproduziert all diese Werte.
Die C++-Datei ist die maßgebliche Implementierung. Sie enthält exakte ganzzahlige Wurzel-Helfer, ein lineares Möbius-Sieb bis \(D\), die Floorsummen-Routine sq_sum_floor und die Tripelzählung summatory_d3. Die Funktion solve summiert zunächst die kleinen \(d\)-Terme direkt, ergänzt dann den gruppierten Rest via Mertens-Werten und verwendet am Ende
$$g(N)=\frac{S(N)+N}{2}.$$
brute_f_from_lcm und formula_f_from_factorization dienen nur den internen Identitätstests für kleine Werte. Die Java-Datei ist eine direkte Übersetzung des Optimierungswegs. Die Python-Datei implementiert den Algorithmus dagegen nicht selbst, sondern kompiliert und startet den C++-Solver und liest anschließend die ausgegebene Antwort aus.
Der Speicherverbrauch wird von der Möbius-Tabelle bis \(D\) und den gruppierten Mertens-Strukturen bis \(I\) bestimmt und liegt deshalb bei ungefähr \(O(N^{1/3})\). Die Laufzeit ist weit von jeder direkten Enumeration entfernt und wird in dieser implementierungsnahen Analyse grob mit etwa \(O(N^{2/3})\) arithmetischen Operationen angesetzt. Genau dadurch wird \(N=10^{12}\) praktikabel.
Her pozitif tamsayı \(n\) için
$$f(n)=\left|\left\{(x,y)\in\mathbb{Z}_{>0}^2 : x\le y,\ \operatorname{lcm}(x,y)=n\right\}\right|,$$
ve toplam fonksiyonu olarak
$$g(N)=\sum_{n=1}^{N} f(n)$$
tanımlanır. İstenen değer \(g(10^{12})\)'dir. \((x,y)\) çiftlerini doğrudan taramak imkansız olduğundan çözüm önce \(f(n)\) için kapalı form çıkarır, sonra toplamı Möbius ağırlıklı bir \(D_3\) toplamına dönüştürür.
\(n\)'yi
$$n=\prod_{p} p^{a_p}$$
şeklinde yazalım. Eğer \(\operatorname{lcm}(x,y)=n\) ise hem \(x\) hem \(y\), \(n\)'nin bölenidir. Her \(p^{a_p}\parallel n\) için \(x\) ve \(y\)'nin \(p\)-adik üsleri \(\alpha_p\) ve \(\beta_p\) olsun. O zaman
$$0\le \alpha_p,\beta_p\le a_p,\qquad \max(\alpha_p,\beta_p)=a_p$$
olmalıdır.
Tek bir \(p^a\) asal kuvveti için uygun sıralı üs çiftlerinin sayısı
$$\underbrace{(a+1)^2}_{0\le \alpha,\beta\le a}-\underbrace{a^2}_{0\le \alpha,\beta\le a-1}=2a+1$$
olur. Asallar birbirinden bağımsız olduğundan, \(\operatorname{lcm}(x,y)=n\) koşulunu sağlayan sıralı çiftlerin toplam sayısı
$$\prod_{p\mid n}(2a_p+1)=\tau(n^2)$$
çıkar; çünkü \(n^2=\prod p^{2a_p}\) ve dolayısıyla \(\tau(n^2)=\prod(2a_p+1)\).
Problem sıralı çiftleri değil, yalnızca \(x\le y\) olan çiftleri sayar. Köşegen dışında her sıralı çift iki kez görünür: \((x,y)\) ve \((y,x)\). En küçük ortak katı \(n\) olan tek köşegen çifti \((n,n)\)'dir. Bu nedenle
$$\boxed{f(n)=\frac{\tau(n^2)+1}{2}.}$$
Bu kapalı form toplanınca
$$g(N)=\sum_{n\le N}\frac{\tau(n^2)+1}{2}=\frac{1}{2}\left(N+\sum_{n\le N}\tau(n^2)\right)$$
elde edilir. Demek ki asıl hedef
$$S(N)=\sum_{n\le N}\tau(n^2)$$
toplamını hesaplamaktır. Sonra
$$g(N)=\frac{S(N)+N}{2}$$
doğrudan bulunur.
\(d_3(m)\), \(m=abc\) biçimindeki sıralı çarpanlaşmaların sayısı olsun. Bunun Dirichlet serisi
$$\sum_{m\ge 1}\frac{d_3(m)}{m^s}=\zeta(s)^3$$
şeklindedir. Ayrıca klasik özdeşlik
$$\sum_{n\ge 1}\frac{\tau(n^2)}{n^s}=\frac{\zeta(s)^3}{\zeta(2s)}$$
vardır. Burada
$$\frac{1}{\zeta(2s)}=\sum_{d\ge 1}\frac{\mu(d)}{d^{2s}}$$
olduğundan katsayı karşılaştırması ile
$$\tau(n^2)=\sum_{d^2m=n}\mu(d)\,d_3(m)=\sum_{d^2\mid n}\mu(d)\,d_3\!\left(\frac{n}{d^2}\right)$$
elde edilir. Şimdi \(n\le N\) üzerinde toplayalım:
$$S(N)=\sum_{d\le \sqrt N}\mu(d)\,D_3\!\left(\left\lfloor\frac{N}{d^2}\right\rfloor\right),$$
burada
$$D_3(X)=\sum_{m\le X} d_3(m)=\left|\left\{(a,b,c)\in\mathbb{Z}_{>0}^3: abc\le X\right\}\right|.$$
Yani problem, hızlı bir \(D_3\) hesaplama ve etkili bir Möbius toplamı problemine dönüşür.
summatory_d3 yordamı \(D_3(X)\)'i hesaplar. Doğrudan tüm üçlüleri saymak yerine, en küçük çarpana göre düzenlenmiş simetrik bir sayım kullanılır. Şunları tanımlayalım:
$$c=\left\lfloor X^{1/3}\right\rfloor,\qquad r_z=\left\lfloor\sqrt{\frac{X}{z}}\right\rfloor.$$
“Üçü eşit”, “yalnızca ikisi eşit” ve “üçü farklı” durumları ayrıldığında programın kullandığı formül ortaya çıkar:
$$D_3(X)=3\sum_{z=1}^{c}\left(2\sum_{x=z+1}^{r_z}\left\lfloor\frac{X}{zx}\right\rfloor-r_z^2+\left\lfloor\frac{X}{z^2}\right\rfloor\right)+c^3.$$
Buradaki \(c^3\) terimi, tekrarlı koordinatların teleskopik olarak sadeleşmiş kalıntısıdır. İçteki
$$\sum_{x=x_1}^{x_2}\left\lfloor\frac{n}{x}\right\rfloor$$
tipi toplam ise doğrudan sq_sum_floor tarafından hesaplanır. C++ ve Java sürümleri burada bölüm değerlerindeki küçük farkları güncelleyerek çok sayıda bölmeyi yeniden yapmaktan kaçınır.
Doğrudan
$$S(N)=\sum_{d\le \sqrt N}\mu(d)\,D_3\!\left(\left\lfloor\frac{N}{d^2}\right\rfloor\right)$$
hesaplamak yine pahalıdır. Kod bu yüzden
$$I=\left\lfloor N^{1/3}\right\rfloor,\qquad D=\left\lfloor\sqrt{\frac{N}{I}}\right\rfloor$$
seçimini yapar. \(d\le D\) için terimler doğrudan alınır. \(d>D\) olduğunda \(\left\lfloor N/d^2\right\rfloor<I\) olur; yani artık çok sayıda \(d\), aynı küçük \(D_3\) argümanına düşer. Mertens fonksiyonunu
$$M(x)=\sum_{n\le x}\mu(n)$$
ve farkları
$$\Delta D_3(i)=D_3(i)-D_3(i-1)=d_3(i)$$
tanımlarsak, büyük kuyruk
$$\sum_{d>D}\mu(d)\,D_3\!\left(\left\lfloor\frac{N}{d^2}\right\rfloor\right)=\sum_{i=1}^{I-1}\Delta D_3(i)\left(M\!\left(\left\lfloor\sqrt{\frac{N}{i}}\right\rfloor\right)-M(D)\right)$$
olarak gruplanır. Kodda small_diff, mertens_small ve mertens_big dizileri bu eşitliği uygular. Özellikle M(D) * D3(I-1) çıkarımı, prefix değerleri istenen farklara dönüştürmek için yapılır.
Program \(\mu\)'yü \(\sqrt N\)'ye kadar elemez. Sadece \(D\)'ye kadar lineer eleme yapar ve \(v=\left\lfloor\sqrt{N/i}\right\rfloor\) biçimindeki daha büyük \(M(v)\) değerlerini gerektiğinde üretir. Kullanılan temel özdeşlik
$$\sum_{d\le v}\mu(d)\left\lfloor\frac{v}{d}\right\rfloor=1$$
şeklindedir. Aralık \(\lfloor\sqrt v\rfloor\)'de bölündüğünde kodun kullandığı özyinelemeli form gelir:
$$M(v)=1-v+\lfloor\sqrt v\rfloor\,M(\lfloor\sqrt v\rfloor)-\sum_{d=2}^{\lfloor\sqrt v\rfloor} M\!\left(\left\lfloor\frac{v}{d}\right\rfloor\right)-\sum_{d=2}^{\lfloor\sqrt v\rfloor}\mu(d)\left\lfloor\frac{v}{d}\right\rfloor.$$
Gerekli argümanlar çok düzenli olduğundan, bunları mertens_big içinde saklamak yeterlidir.
\(n=12=2^2\cdot 3\) için
$$\tau(12^2)=\tau(2^4\cdot 3^2)=(4+1)(2+1)=15$$
ve dolayısıyla
$$f(12)=\frac{15+1}{2}=8$$
çıkar. C++ kodundaki kontrol değerleri de
$$g(10)=29,\qquad g(100)=647,\qquad g(1000)=11751,\qquad g(10^4)=186991,\qquad g(10^6)=37429395$$
şeklindedir ve optimize çözüm hepsini doğrular.
Asıl uygulama C++ dosyasında bulunur. Burada tam sayı karekök ve küpkök yardımcıları, \(D\)'ye kadar lineer Möbius eleği, sq_sum_floor tabanlı hızlı floor-toplam rutini ve summatory_d3 yer alır. solve önce küçük \(d\) terimlerini toplar, sonra Mertens gruplamasıyla büyük kuyruğu ekler ve sonunda
$$g(N)=\frac{S(N)+N}{2}$$
ilişkisini uygular. brute_f_from_lcm ile formula_f_from_factorization yalnızca küçük değerlerde özdeşlik testi içindir. Java dosyası aynı algoritmanın doğrudan çevirisidir. Python dosyası ise farklı bir rol üstlenir: C++ çözümünü derleyip çalıştıran bir köprü görevi görür ve çıktıyı ayrıştırır.
Bellek tüketimini \(D\)'ye kadar Möbius tablosu ve \(I\)'ye kadar Mertens grupları belirler; bu yüzden yaklaşık \(O(N^{1/3})\) bellek kullanılır. Zaman karmaşıklığı doğrudan taramadan çok daha iyidir ve bu yaklaşım için kullanılan uygulama odaklı tahminde yaklaşık \(O(N^{2/3})\) aritmetik işlem mertebesindedir. Bu sayede \(N=10^{12}\) erişilebilir hale gelir.
Para cada entero positivo \(n\), se define
$$f(n)=\left|\left\{(x,y)\in\mathbb{Z}_{>0}^2 : x\le y,\ \operatorname{lcm}(x,y)=n\right\}\right|,$$
y la función acumulada
$$g(N)=\sum_{n=1}^{N} f(n).$$
La tarea pide calcular \(g(10^{12})\). Un barrido directo de todos los pares \((x,y)\) es imposible, así que la solución primero obtiene una fórmula cerrada para \(f(n)\) y luego transforma la suma en una expresión con la función de Möbius y la función divisor ternaria.
Escribamos
$$n=\prod_{p} p^{a_p}.$$
Si \(\operatorname{lcm}(x,y)=n\), entonces \(x\) e \(y\) son divisores de \(n\). Para cada primo \(p^{a_p}\parallel n\), sean \(\alpha_p\) y \(\beta_p\) los exponentes de \(p\) en \(x\) e \(y\). Deben cumplir
$$0\le \alpha_p,\beta_p\le a_p,\qquad \max(\alpha_p,\beta_p)=a_p.$$
Para una sola potencia prima \(p^a\), el número de pares ordenados \((\alpha,\beta)\) admisibles es
$$\underbrace{(a+1)^2}_{0\le \alpha,\beta\le a}-\underbrace{a^2}_{0\le \alpha,\beta\le a-1}=2a+1.$$
Como los distintos primos son independientes, el número total de pares ordenados \((x,y)\) con \(\operatorname{lcm}(x,y)=n\) es
$$\prod_{p\mid n}(2a_p+1)=\tau(n^2),$$
porque \(n^2=\prod p^{2a_p}\) y, por tanto, \(\tau(n^2)=\prod(2a_p+1)\).
El problema no cuenta pares ordenados, sino solo los que satisfacen \(x\le y\). Todo par no diagonal aparece dos veces, como \((x,y)\) y \((y,x)\). El único par diagonal con mínimo común múltiplo igual a \(n\) es \((n,n)\). Entonces
$$\boxed{f(n)=\frac{\tau(n^2)+1}{2}.}$$
Al sumar la fórmula cerrada obtenemos
$$g(N)=\sum_{n\le N}\frac{\tau(n^2)+1}{2}=\frac{1}{2}\left(N+\sum_{n\le N}\tau(n^2)\right).$$
Por lo tanto, el núcleo del problema es
$$S(N)=\sum_{n\le N}\tau(n^2).$$
Una vez calculado \(S(N)\), se recupera
$$g(N)=\frac{S(N)+N}{2}.$$
Sea \(d_3(m)\) el número de factorizaciones ordenadas \(m=abc\). Su serie de Dirichlet es
$$\sum_{m\ge 1}\frac{d_3(m)}{m^s}=\zeta(s)^3.$$
También se tiene la identidad clásica
$$\sum_{n\ge 1}\frac{\tau(n^2)}{n^s}=\frac{\zeta(s)^3}{\zeta(2s)}.$$
Usando
$$\frac{1}{\zeta(2s)}=\sum_{d\ge 1}\frac{\mu(d)}{d^{2s}},$$
y comparando coeficientes, resulta
$$\tau(n^2)=\sum_{d^2m=n}\mu(d)\,d_3(m)=\sum_{d^2\mid n}\mu(d)\,d_3\!\left(\frac{n}{d^2}\right).$$
Sumando sobre \(n\le N\) e intercambiando el orden de suma, se llega a
$$S(N)=\sum_{d\le \sqrt N}\mu(d)\,D_3\!\left(\left\lfloor\frac{N}{d^2}\right\rfloor\right),$$
donde
$$D_3(X)=\sum_{m\le X} d_3(m)=\left|\left\{(a,b,c)\in\mathbb{Z}_{>0}^3: abc\le X\right\}\right|.$$
La función summatory_d3 calcula \(D_3(X)\). En lugar de enumerar todos los triples, ordena por el menor factor. Definimos
$$c=\left\lfloor X^{1/3}\right\rfloor,\qquad r_z=\left\lfloor\sqrt{\frac{X}{z}}\right\rfloor.$$
Separando los casos “los tres iguales”, “exactamente dos iguales” y “todos distintos”, se obtiene la fórmula usada por el código:
$$D_3(X)=3\sum_{z=1}^{c}\left(2\sum_{x=z+1}^{r_z}\left\lfloor\frac{X}{zx}\right\rfloor-r_z^2+\left\lfloor\frac{X}{z^2}\right\rfloor\right)+c^3.$$
El término \(c^3\) resume de forma telescópica la contribución de los triples con coordenadas repetidas. La suma interior
$$\sum_{x=x_1}^{x_2}\left\lfloor\frac{n}{x}\right\rfloor$$
es precisamente lo que evalúa sq_sum_floor en C++ y Java, usando actualizaciones incrementales del cociente para evitar divisiones redundantes.
La identidad
$$S(N)=\sum_{d\le \sqrt N}\mu(d)\,D_3\!\left(\left\lfloor\frac{N}{d^2}\right\rfloor\right)$$
sigue siendo demasiado costosa si se recorre de manera ingenua. La implementación usa
$$I=\left\lfloor N^{1/3}\right\rfloor,\qquad D=\left\lfloor\sqrt{\frac{N}{I}}\right\rfloor.$$
Para \(d\le D\), los términos se calculan directamente. Para \(d>D\), el valor \(\left\lfloor N/d^2\right\rfloor\) es menor que \(I\), así que muchos valores de \(d\) comparten el mismo argumento pequeño. Introducimos la función de Mertens
$$M(x)=\sum_{n\le x}\mu(n)$$
y las diferencias
$$\Delta D_3(i)=D_3(i)-D_3(i-1)=d_3(i).$$
Entonces la cola se puede reagrupar como
$$\sum_{d>D}\mu(d)\,D_3\!\left(\left\lfloor\frac{N}{d^2}\right\rfloor\right)=\sum_{i=1}^{I-1}\Delta D_3(i)\left(M\!\left(\left\lfloor\sqrt{\frac{N}{i}}\right\rfloor\right)-M(D)\right).$$
Eso es exactamente lo que representan los arreglos small_diff, mertens_small y mertens_big en la implementación.
El programa no criba \(\mu(d)\) hasta \(\sqrt N\). Solo lo hace hasta \(D\), y luego calcula de forma recursiva los valores grandes \(M(v)\) para \(v=\left\lfloor\sqrt{N/i}\right\rfloor\). La identidad base es
$$\sum_{d\le v}\mu(d)\left\lfloor\frac{v}{d}\right\rfloor=1.$$
Separando en \(\lfloor\sqrt v\rfloor\), se obtiene la forma usada en el código:
$$M(v)=1-v+\lfloor\sqrt v\rfloor\,M(\lfloor\sqrt v\rfloor)-\sum_{d=2}^{\lfloor\sqrt v\rfloor} M\!\left(\left\lfloor\frac{v}{d}\right\rfloor\right)-\sum_{d=2}^{\lfloor\sqrt v\rfloor}\mu(d)\left\lfloor\frac{v}{d}\right\rfloor.$$
Como los argumentos requeridos siguen una estructura muy restringida, basta con memorizarlos en mertens_big.
Si \(n=12=2^2\cdot 3\), entonces
$$\tau(12^2)=\tau(2^4\cdot 3^2)=(4+1)(2+1)=15,$$
y por tanto
$$f(12)=\frac{15+1}{2}=8.$$
Los puntos de control incluidos en el archivo C++ son
$$g(10)=29,\qquad g(100)=647,\qquad g(1000)=11751,\qquad g(10^4)=186991,\qquad g(10^6)=37429395,$$
todos verificados por la versión optimizada antes de atacar \(10^{12}\).
La fuente C++ es la implementación principal. Contiene raíces cuadrada y cúbica enteras exactas, una criba lineal para \(\mu\) hasta \(D\), el ayudante de sumas con parte entera sq_sum_floor y la rutina summatory_d3. La función solve suma primero la parte pequeña en \(d\), luego añade la cola agrupada mediante valores de Mertens y finalmente usa
$$g(N)=\frac{S(N)+N}{2}.$$
Las funciones brute_f_from_lcm y formula_f_from_factorization solo sirven para validar la identidad en entradas pequeñas. El archivo Java reproduce el mismo algoritmo optimizado. El archivo Python no reimplementa la matemática: actúa como un puente que compila y ejecuta la solución C++ y extrae la respuesta impresa.
La memoria está dominada por la criba de Möbius hasta \(D\) y por las tablas agrupadas de Mertens hasta \(I\), así que el uso es aproximadamente \(O(N^{1/3})\). El tiempo está muy por debajo de cualquier enumeración directa y, en la estimación práctica empleada para este enfoque, queda alrededor de \(O(N^{2/3})\) operaciones aritméticas con subrutinas rápidas de sumas de pisos.
对每个正整数 \(n\)，定义
$$f(n)=\left|\left\{(x,y)\in\mathbb{Z}_{>0}^2 : x\le y,\ \operatorname{lcm}(x,y)=n\right\}\right|,$$
再定义累加函数
$$g(N)=\sum_{n=1}^{N} f(n).$$
题目要求的是 \(g(10^{12})\)。如果直接枚举所有 \((x,y)\) 对，规模完全不可接受。因此程序先推导 \(f(n)\) 的闭式，再把总和改写成带有 Möbius 函数权重的三元除数和。
把
$$n=\prod_{p} p^{a_p}$$
写成标准质因数分解。若 \(\operatorname{lcm}(x,y)=n\)，那么 \(x\) 与 \(y\) 都必须是 \(n\) 的约数。对每个满足 \(p^{a_p}\parallel n\) 的质数，设 \(x\) 与 \(y\) 中对应的指数分别为 \(\alpha_p,\beta_p\)。它们必须满足
$$0\le \alpha_p,\beta_p\le a_p,\qquad \max(\alpha_p,\beta_p)=a_p.$$
对于单个质数幂 \(p^a\)，合法的有序指数对个数为
$$\underbrace{(a+1)^2}_{0\le \alpha,\beta\le a}-\underbrace{a^2}_{0\le \alpha,\beta\le a-1}=2a+1.$$
不同质数彼此独立，因此满足 \(\operatorname{lcm}(x,y)=n\) 的有序对总数是
$$\prod_{p\mid n}(2a_p+1)=\tau(n^2),$$
因为 \(n^2=\prod p^{2a_p}\)，从而 \(\tau(n^2)=\prod(2a_p+1)\)。
但题目只统计 \(x\le y\) 的对，而不是所有有序对。所有非对角项都会成对出现：\((x,y)\) 和 \((y,x)\)。唯一满足 \(\operatorname{lcm}(x,x)=n\) 的对角项只有 \((n,n)\)。因此
$$\boxed{f(n)=\frac{\tau(n^2)+1}{2}.}$$
对上式求和可得
$$g(N)=\sum_{n\le N}\frac{\tau(n^2)+1}{2}=\frac{1}{2}\left(N+\sum_{n\le N}\tau(n^2)\right).$$
所以核心对象是
$$S(N)=\sum_{n\le N}\tau(n^2).$$
一旦知道 \(S(N)\)，就有
$$g(N)=\frac{S(N)+N}{2}.$$
记 \(d_3(m)\) 为 \(m=abc\) 的有序分解数。它的 Dirichlet 级数为
$$\sum_{m\ge 1}\frac{d_3(m)}{m^s}=\zeta(s)^3.$$
而 \(\tau(n^2)\) 满足经典恒等式
$$\sum_{n\ge 1}\frac{\tau(n^2)}{n^s}=\frac{\zeta(s)^3}{\zeta(2s)}.$$
再利用
$$\frac{1}{\zeta(2s)}=\sum_{d\ge 1}\frac{\mu(d)}{d^{2s}},$$
比较系数得到
$$\tau(n^2)=\sum_{d^2m=n}\mu(d)\,d_3(m)=\sum_{d^2\mid n}\mu(d)\,d_3\!\left(\frac{n}{d^2}\right).$$
于是
$$S(N)=\sum_{d\le \sqrt N}\mu(d)\,D_3\!\left(\left\lfloor\frac{N}{d^2}\right\rfloor\right),$$
其中
$$D_3(X)=\sum_{m\le X} d_3(m)=\left|\left\{(a,b,c)\in\mathbb{Z}_{>0}^3: abc\le X\right\}\right|.$$
这一步把问题从 \(\tau(n^2)\) 的求和转化为一个更适合程序实现的三维计数问题。
辅助函数 summatory_d3 负责计算 \(D_3(X)\)。它并不暴力枚举三元组，而是按最小因子进行对称分类。令
$$c=\left\lfloor X^{1/3}\right\rfloor,\qquad r_z=\left\lfloor\sqrt{\frac{X}{z}}\right\rfloor.$$
把“三个都相等”“恰有两个相等”“三个互不相同”三种情况整理后，可得到代码实现的公式：
$$D_3(X)=3\sum_{z=1}^{c}\left(2\sum_{x=z+1}^{r_z}\left\lfloor\frac{X}{zx}\right\rfloor-r_z^2+\left\lfloor\frac{X}{z^2}\right\rfloor\right)+c^3.$$
最后的 \(c^3\) 项是重复坐标贡献经过望远镜化简后的结果。内部的整除求和
$$\sum_{x=x_1}^{x_2}\left\lfloor\frac{n}{x}\right\rfloor$$
正是 sq_sum_floor 的任务。C++ 和 Java 版本都通过商值差分的增量更新来避免重复做大量整除运算。
即使有了
$$S(N)=\sum_{d\le \sqrt N}\mu(d)\,D_3\!\left(\left\lfloor\frac{N}{d^2}\right\rfloor\right),$$
若仍然逐个遍历 \(d\)，成本也太高。实现中采用
$$I=\left\lfloor N^{1/3}\right\rfloor,\qquad D=\left\lfloor\sqrt{\frac{N}{I}}\right\rfloor.$$
当 \(d\le D\) 时，直接计算对应项。对 \(d>D\)，有 \(\left\lfloor N/d^2\right\rfloor<I\)，因此许多 \(d\) 会映射到同一个较小的 \(D_3\) 参数。引入 Mertens 函数
$$M(x)=\sum_{n\le x}\mu(n)$$
以及差分
$$\Delta D_3(i)=D_3(i)-D_3(i-1)=d_3(i).$$
那么大尾项可以改写为
$$\sum_{d>D}\mu(d)\,D_3\!\left(\left\lfloor\frac{N}{d^2}\right\rfloor\right)=\sum_{i=1}^{I-1}\Delta D_3(i)\left(M\!\left(\left\lfloor\sqrt{\frac{N}{i}}\right\rfloor\right)-M(D)\right).$$
程序中的 small_diff、mertens_small 和 mertens_big 正是对这条公式的直接实现。
程序不会把 \(\mu\) 一直筛到 \(\sqrt N\)。它只筛到 \(D\)，然后对需要的 \(v=\left\lfloor\sqrt{N/i}\right\rfloor\) 递归计算 \(M(v)\)。基础恒等式是
$$\sum_{d\le v}\mu(d)\left\lfloor\frac{v}{d}\right\rfloor=1.$$
在 \(\lfloor\sqrt v\rfloor\) 处分段后，可得到代码里使用的形式：
$$M(v)=1-v+\lfloor\sqrt v\rfloor\,M(\lfloor\sqrt v\rfloor)-\sum_{d=2}^{\lfloor\sqrt v\rfloor} M\!\left(\left\lfloor\frac{v}{d}\right\rfloor\right)-\sum_{d=2}^{\lfloor\sqrt v\rfloor}\mu(d)\left\lfloor\frac{v}{d}\right\rfloor.$$
由于真正用到的参数形状非常规则，只需要把这些值缓存在 mertens_big 中即可。
例如 \(n=12=2^2\cdot 3\)，那么
$$\tau(12^2)=\tau(2^4\cdot 3^2)=(4+1)(2+1)=15,$$
所以
$$f(12)=\frac{15+1}{2}=8.$$
C++ 程序内置的检查点还包括
$$g(10)=29,\qquad g(100)=647,\qquad g(1000)=11751,\qquad g(10^4)=186991,\qquad g(10^6)=37429395,$$
优化算法在处理 \(10^{12}\) 之前会先通过这些验证。
C++ 文件是权威实现。它包含精确的整数平方根与立方根函数、到 \(D\) 为止的线性 Möbius 筛、整除和助手 sq_sum_floor 以及三元计数函数 summatory_d3。主函数 solve 先累加小 \(d\) 的直接部分，再通过分组后的 Mertens 值补上大尾项，最后使用
$$g(N)=\frac{S(N)+N}{2}$$
得到答案。brute_f_from_lcm 和 formula_f_from_factorization 只用于小范围正确性检查。Java 文件是同一算法的直接翻译。Python 文件则不是独立实现，而是一个桥接器：它在需要时编译并运行 C++ 求解器，然后解析输出。
内存主要由到 \(D\) 的 Möbius 表和到 \(I\) 的分组 Mertens 数组占用，因此约为 \(O(N^{1/3})\)。运行时间远优于直接枚举，在这种实现导向的估计下，大致可视为 \(O(N^{2/3})\) 级别，并配合快速整除和子程序使 \(N=10^{12}\) 成为可计算规模。
Для каждого положительного целого \(n\) определим
$$f(n)=\left|\left\{(x,y)\in\mathbb{Z}_{>0}^2 : x\le y,\ \operatorname{lcm}(x,y)=n\right\}\right|,$$
а затем суммирующую функцию
$$g(N)=\sum_{n=1}^{N} f(n).$$
Нужно найти \(g(10^{12})\). Перебор всех пар \((x,y)\) невозможен, поэтому решение сначала выводит явную формулу для \(f(n)\), а затем сводит задачу к сумме с функцией Мёбиуса и троичной функцией делителей.
Пусть
$$n=\prod_{p} p^{a_p}.$$
Если \(\operatorname{lcm}(x,y)=n\), то \(x\) и \(y\) являются делителями \(n\). Для каждого простого \(p^{a_p}\parallel n\) обозначим через \(\alpha_p\) и \(\beta_p\) показатели степени \(p\) в \(x\) и \(y\). Тогда необходимо
$$0\le \alpha_p,\beta_p\le a_p,\qquad \max(\alpha_p,\beta_p)=a_p.$$
Для одной фиксированной степени \(p^a\) число допустимых упорядоченных пар \((\alpha,\beta)\) равно
$$\underbrace{(a+1)^2}_{0\le \alpha,\beta\le a}-\underbrace{a^2}_{0\le \alpha,\beta\le a-1}=2a+1.$$
Поскольку простые множители независимы, общее число упорядоченных пар \((x,y)\) с \(\operatorname{lcm}(x,y)=n\) равно
$$\prod_{p\mid n}(2a_p+1)=\tau(n^2),$$
так как \(n^2=\prod p^{2a_p}\) и, следовательно, \(\tau(n^2)=\prod(2a_p+1)\).
Но задача считает только пары с условием \(x\le y\). Каждая недиагональная упорядоченная пара встречается дважды: как \((x,y)\) и \((y,x)\). Единственная диагональная пара с НОК, равным \(n\), это \((n,n)\). Поэтому
$$\boxed{f(n)=\frac{\tau(n^2)+1}{2}.}$$
Суммируя, получаем
$$g(N)=\sum_{n\le N}\frac{\tau(n^2)+1}{2}=\frac{1}{2}\left(N+\sum_{n\le N}\tau(n^2)\right).$$
Значит, нужно эффективно вычислить
$$S(N)=\sum_{n\le N}\tau(n^2).$$
После этого
$$g(N)=\frac{S(N)+N}{2}$$
получается сразу.
Пусть \(d_3(m)\) обозначает число упорядоченных разложений \(m=abc\). Ее ряд Дирихле имеет вид
$$\sum_{m\ge 1}\frac{d_3(m)}{m^s}=\zeta(s)^3.$$
Для \(\tau(n^2)\) используется тождество
$$\sum_{n\ge 1}\frac{\tau(n^2)}{n^s}=\frac{\zeta(s)^3}{\zeta(2s)}.$$
Так как
$$\frac{1}{\zeta(2s)}=\sum_{d\ge 1}\frac{\mu(d)}{d^{2s}},$$
сравнение коэффициентов дает
$$\tau(n^2)=\sum_{d^2m=n}\mu(d)\,d_3(m)=\sum_{d^2\mid n}\mu(d)\,d_3\!\left(\frac{n}{d^2}\right).$$
Отсюда следует
$$S(N)=\sum_{d\le \sqrt N}\mu(d)\,D_3\!\left(\left\lfloor\frac{N}{d^2}\right\rfloor\right),$$
где
$$D_3(X)=\sum_{m\le X} d_3(m)=\left|\left\{(a,b,c)\in\mathbb{Z}_{>0}^3: abc\le X\right\}\right|.$$
Функция summatory_d3 вычисляет \(D_3(X)\). Вместо полного перебора троек она использует симметричное разбиение по наименьшему множителю. Обозначим
$$c=\left\lfloor X^{1/3}\right\rfloor,\qquad r_z=\left\lfloor\sqrt{\frac{X}{z}}\right\rfloor.$$
После выделения случаев «все три равны», «ровно два равны» и «все различны» получается формула, реализованная в коде:
$$D_3(X)=3\sum_{z=1}^{c}\left(2\sum_{x=z+1}^{r_z}\left\lfloor\frac{X}{zx}\right\rfloor-r_z^2+\left\lfloor\frac{X}{z^2}\right\rfloor\right)+c^3.$$
Член \(c^3\) представляет собой компактную телескопическую запись вклада троек с повторяющимися координатами. Внутренняя сумма вида
$$\sum_{x=x_1}^{x_2}\left\lfloor\frac{n}{x}\right\rfloor$$
вычисляется вспомогательной функцией sq_sum_floor, которая обновляет частные инкрементально, не выполняя каждое деление заново.
Даже формула
$$S(N)=\sum_{d\le \sqrt N}\mu(d)\,D_3\!\left(\left\lfloor\frac{N}{d^2}\right\rfloor\right)$$
остается слишком дорогой при наивном проходе по всем \(d\). Поэтому реализация выбирает
$$I=\left\lfloor N^{1/3}\right\rfloor,\qquad D=\left\lfloor\sqrt{\frac{N}{I}}\right\rfloor.$$
Для \(d\le D\) значения суммируются напрямую. Для \(d>D\) величина \(\left\lfloor N/d^2\right\rfloor\) уже меньше \(I\), поэтому многие соседние \(d\) дают один и тот же малый аргумент. Введем функцию Мертенса
$$M(x)=\sum_{n\le x}\mu(n)$$
и разности
$$\Delta D_3(i)=D_3(i)-D_3(i-1)=d_3(i).$$
Тогда хвост можно переписать так:
$$\sum_{d>D}\mu(d)\,D_3\!\left(\left\lfloor\frac{N}{d^2}\right\rfloor\right)=\sum_{i=1}^{I-1}\Delta D_3(i)\left(M\!\left(\left\lfloor\sqrt{\frac{N}{i}}\right\rfloor\right)-M(D)\right).$$
Это и реализуют массивы small_diff, mertens_small и mertens_big в программе.
Программа не строит решето \(\mu\) до \(\sqrt N\). Она просеивает только до \(D\), а большие значения \(M(v)\) для \(v=\left\lfloor\sqrt{N/i}\right\rfloor\) находит по рекуррентной формуле. Исходное тождество таково:
$$\sum_{d\le v}\mu(d)\left\lfloor\frac{v}{d}\right\rfloor=1.$$
После разделения диапазона на уровне \(\lfloor\sqrt v\rfloor\) получаем форму, использованную в коде:
$$M(v)=1-v+\lfloor\sqrt v\rfloor\,M(\lfloor\sqrt v\rfloor)-\sum_{d=2}^{\lfloor\sqrt v\rfloor} M\!\left(\left\lfloor\frac{v}{d}\right\rfloor\right)-\sum_{d=2}^{\lfloor\sqrt v\rfloor}\mu(d)\left\lfloor\frac{v}{d}\right\rfloor.$$
Поскольку реально нужны только значения специального вида \(\left\lfloor\sqrt{N/i}\right\rfloor\), их достаточно хранить в массиве mertens_big.
Например, для \(n=12=2^2\cdot 3\) имеем
$$\tau(12^2)=\tau(2^4\cdot 3^2)=(4+1)(2+1)=15,$$
поэтому
$$f(12)=\frac{15+1}{2}=8.$$
Контрольные значения, встроенные в C++-файл, таковы:
$$g(10)=29,\qquad g(100)=647,\qquad g(1000)=11751,\qquad g(10^4)=186991,\qquad g(10^6)=37429395,$$
и оптимизированный алгоритм подтверждает их все.
Файл C++ является основной реализацией. В нем есть точные целочисленные функции квадратного и кубического корня, линейное решето Мёбиуса до \(D\), функция sq_sum_floor для быстрых сумм с полом и процедура summatory_d3. Основная функция solve сначала обрабатывает малые \(d\), затем добавляет сгруппированный хвост через значения Мертенса и в конце использует
$$g(N)=\frac{S(N)+N}{2}.$$
Функции brute_f_from_lcm и formula_f_from_factorization нужны только для внутренних проверок на малых входах. Java-файл является прямым переводом того же оптимизированного алгоритма. Python-файл не повторяет математику, а выступает мостом: при необходимости он компилирует и запускает C++-решение и разбирает напечатанный ответ.
Основную память занимают таблица Мёбиуса до \(D\) и массивы сгруппированных значений Мертенса до \(I\), то есть примерно \(O(N^{1/3})\). Временная сложность гораздо лучше любого прямого перебора и в практической оценке для данного подхода составляет около \(O(N^{2/3})\) арифметических операций с быстрыми подпрограммами для сумм целых частей.
لكل عدد صحيح موجب \(n\) نعرّف
$$f(n)=\left|\left\{(x,y)\in\mathbb{Z}_{>0}^2 : x\le y,\ \operatorname{lcm}(x,y)=n\right\}\right|,$$
ثم نعرّف الدالة التجميعية
$$g(N)=\sum_{n=1}^{N} f(n).$$
المطلوب هو حساب \(g(10^{12})\). العد المباشر لكل الأزواج \((x,y)\) غير ممكن عمليًا، لذلك يبدأ الحل باشتقاق صيغة مغلقة لـ \(f(n)\)، ثم يحوّل المجموع إلى صيغة تعتمد على دالة موبيوس وعلى دالة القواسم الثلاثية.
لنكتب
$$n=\prod_{p} p^{a_p}.$$
إذا كان \(\operatorname{lcm}(x,y)=n\)، فلابد أن يكون كل من \(x\) و\(y\) قاسمًا لـ \(n\). ولكل عامل أولي \(p^{a_p}\parallel n\)، نرمز إلى الأسين الموافقين في \(x\) و\(y\) بالرمزين \(\alpha_p\) و\(\beta_p\). عندئذ يجب أن يتحقق
$$0\le \alpha_p,\beta_p\le a_p,\qquad \max(\alpha_p,\beta_p)=a_p.$$
ولعامل أولي واحد \(p^a\)، يكون عدد أزواج الأسس المرتبة الممكنة
$$\underbrace{(a+1)^2}_{0\le \alpha,\beta\le a}-\underbrace{a^2}_{0\le \alpha,\beta\le a-1}=2a+1.$$
وبما أن مساهمات العوامل الأولية مستقلة، فإن عدد الأزواج المرتبة \((x,y)\) التي تحقق \(\operatorname{lcm}(x,y)=n\) هو
$$\prod_{p\mid n}(2a_p+1)=\tau(n^2),$$
لأن \(n^2=\prod p^{2a_p}\)، ومن ثم \(\tau(n^2)=\prod(2a_p+1)\).
لكن المطلوب في المسألة ليس الأزواج المرتبة كلها، بل الأزواج التي تحقق \(x\le y\) فقط. كل زوج غير قطري يظهر مرتين: مرة على هيئة \((x,y)\) ومرة على هيئة \((y,x)\). أما الزوج القطري الوحيد الذي يعطي المضاعف المشترك الأصغر \(n\) فهو \((n,n)\). لذلك نحصل على
$$\boxed{f(n)=\frac{\tau(n^2)+1}{2}.}$$
بجمع الصيغة السابقة نحصل على
$$g(N)=\sum_{n\le N}\frac{\tau(n^2)+1}{2}=\frac{1}{2}\left(N+\sum_{n\le N}\tau(n^2)\right).$$
إذًا المسألة الأساسية تصبح حساب
$$S(N)=\sum_{n\le N}\tau(n^2).$$
وبعدها يكون
$$g(N)=\frac{S(N)+N}{2}.$$
لنرمز بـ \(d_3(m)\) إلى عدد التحليلات المرتبة \(m=abc\). متسلسلتها الديريشليه هي
$$\sum_{m\ge 1}\frac{d_3(m)}{m^s}=\zeta(s)^3.$$
كما أن لدينا الهوية المعروفة
$$\sum_{n\ge 1}\frac{\tau(n^2)}{n^s}=\frac{\zeta(s)^3}{\zeta(2s)}.$$
وباستخدام
$$\frac{1}{\zeta(2s)}=\sum_{d\ge 1}\frac{\mu(d)}{d^{2s}},$$
ثم مقارنة المعاملات، نحصل على
$$\tau(n^2)=\sum_{d^2m=n}\mu(d)\,d_3(m)=\sum_{d^2\mid n}\mu(d)\,d_3\!\left(\frac{n}{d^2}\right).$$
ومن ثم
$$S(N)=\sum_{d\le \sqrt N}\mu(d)\,D_3\!\left(\left\lfloor\frac{N}{d^2}\right\rfloor\right),$$
حيث
$$D_3(X)=\sum_{m\le X} d_3(m)=\left|\left\{(a,b,c)\in\mathbb{Z}_{>0}^3: abc\le X\right\}\right|.$$
وهكذا تتحول المسألة من مجموع على \(\tau(n^2)\) إلى عد ثلاثيات يمكن التعامل معه بسرعة أكبر.
الدالة المساعدة summatory_d3 تحسب \(D_3(X)\). بدلًا من تعداد كل الثلاثيات مباشرة، تستغل التناظر بحسب أصغر عامل. نعرّف
$$c=\left\lfloor X^{1/3}\right\rfloor,\qquad r_z=\left\lfloor\sqrt{\frac{X}{z}}\right\rfloor.$$
بعد فصل الحالات: “الثلاثة متساوية”، و“اثنان فقط متساويان”، و“الثلاثة مختلفة”، نحصل على الصيغة المستعملة في الكود:
$$D_3(X)=3\sum_{z=1}^{c}\left(2\sum_{x=z+1}^{r_z}\left\lfloor\frac{X}{zx}\right\rfloor-r_z^2+\left\lfloor\frac{X}{z^2}\right\rfloor\right)+c^3.$$
الحد \(c^3\) هو الشكل المضغوط للتجميع التلسكوبي الناتج من الثلاثيات ذات الإحداثيات المكررة. أما المجموع الداخلي من النوع
$$\sum_{x=x_1}^{x_2}\left\lfloor\frac{n}{x}\right\rfloor$$
فهو بالضبط ما تحسبه الدالة sq_sum_floor في نسختي C++ وJava، وذلك بتحديث فروق القيم بدل إعادة القسمة في كل مرة.
حتى بعد الوصول إلى
$$S(N)=\sum_{d\le \sqrt N}\mu(d)\,D_3\!\left(\left\lfloor\frac{N}{d^2}\right\rfloor\right),$$
يبقى المرور الساذج على كل \(d\) مكلفًا جدًا. لذلك تستخدم الخوارزمية التقسيم
$$I=\left\lfloor N^{1/3}\right\rfloor,\qquad D=\left\lfloor\sqrt{\frac{N}{I}}\right\rfloor.$$
عندما يكون \(d\le D\)، تُحسب الحدود مباشرة. أما إذا كان \(d>D\)، فإن \(\left\lfloor N/d^2\right\rfloor<I\)، وهذا يعني أن عددًا كبيرًا من قيم \(d\) يعطي القيم الصغيرة نفسها للدالة \(D_3\). نعرّف دالة ميرتنز
$$M(x)=\sum_{n\le x}\mu(n)$$
وكذلك الفروق
$$\Delta D_3(i)=D_3(i)-D_3(i-1)=d_3(i).$$
عندئذ يمكن إعادة كتابة الذيل الكبير بالشكل
$$\sum_{d>D}\mu(d)\,D_3\!\left(\left\lfloor\frac{N}{d^2}\right\rfloor\right)=\sum_{i=1}^{I-1}\Delta D_3(i)\left(M\!\left(\left\lfloor\sqrt{\frac{N}{i}}\right\rfloor\right)-M(D)\right).$$
وهذا بالضبط ما تنفذه المصفوفات small_diff وmertens_small وmertens_big في التطبيق.
لا يقوم البرنامج بعمل غربال كامل لـ \(\mu\) حتى \(\sqrt N\). بل يكتفي بالغربلة حتى \(D\)، ثم يحسب القيم الكبيرة \(M(v)\) المطلوبة عند \(v=\left\lfloor\sqrt{N/i}\right\rfloor\) باستخدام علاقة عودية. الهوية الأساسية هي
$$\sum_{d\le v}\mu(d)\left\lfloor\frac{v}{d}\right\rfloor=1.$$
وبعد فصل المدى عند \(\lfloor\sqrt v\rfloor\)، نحصل على الصيغة المستعملة في الكود:
$$M(v)=1-v+\lfloor\sqrt v\rfloor\,M(\lfloor\sqrt v\rfloor)-\sum_{d=2}^{\lfloor\sqrt v\rfloor} M\!\left(\left\lfloor\frac{v}{d}\right\rfloor\right)-\sum_{d=2}^{\lfloor\sqrt v\rfloor}\mu(d)\left\lfloor\frac{v}{d}\right\rfloor.$$
وبما أن القيم المطلوبة تأتي على صورة منظمة جدًا، فيكفي تخزينها في mertens_big من دون الحاجة إلى غربال ضخم.
إذا أخذنا \(n=12=2^2\cdot 3\)، فإن
$$\tau(12^2)=\tau(2^4\cdot 3^2)=(4+1)(2+1)=15,$$
ومن ثم
$$f(12)=\frac{15+1}{2}=8.$$
كما أن قيم التحقق الموجودة في ملف C++ هي
$$g(10)=29,\qquad g(100)=647,\qquad g(1000)=11751,\qquad g(10^4)=186991,\qquad g(10^6)=37429395,$$
وتتطابق معها الخوارزمية المحسنة قبل الانتقال إلى \(10^{12}\).
ملف C++ هو المرجع الأساسي للتنفيذ. يحتوي على دوال دقيقة للجذر التربيعي والتكعيبي الصحيحين، وغربال خطي لدالة موبيوس حتى \(D\)، والدالة sq_sum_floor لمجاميع الجزء الصحيح، والدالة summatory_d3 لعد الثلاثيات. تقوم الدالة solve أولًا بجمع مساهمات \(d\) الصغيرة مباشرة، ثم تضيف الذيل المجمع باستخدام قيم ميرتنز، وفي النهاية تطبق العلاقة
$$g(N)=\frac{S(N)+N}{2}.$$
أما الدالتان brute_f_from_lcm وformula_f_from_factorization فهما مخصصتان فقط لاختبارات الصحة على القيم الصغيرة. ملف Java هو ترجمة مباشرة للخوارزمية نفسها. أما ملف Python فله دور مختلف: إنه جسر خفيف يقوم بترجمة برنامج C++ وتشغيله ثم استخراج الجواب من المخرجات.
يهيمن على الذاكرة جدول موبيوس حتى \(D\) وبنى ميرتنز المجمعة حتى \(I\)، لذلك يكون الاستهلاك تقريبًا \(O(N^{1/3})\). أما الزمن فهو أفضل بكثير من أي تعداد مباشر، وفي التقدير العملي المستخدم لهذا النهج يمكن اعتباره تقريبًا \(O(N^{2/3})\) من العمليات الحسابية مع روتينات سريعة لمجاميع الأجزاء الصحيحة.