We must evaluate
$$S(N)=\sum_{n=2}^{N}\gcd(n,n'),$$
where \(n'\) is the arithmetic derivative, defined by \(p'=1\) for every prime \(p\) and by the product rule \((ab)'=a'b+ab'\). For \(n=\prod p^{e_p}\), this implies
$$n'=n\sum_{p^{e_p}\parallel n}\frac{e_p}{p}.$$
The target value uses \(N=5\times 10^{15}\), so factoring every integer up to \(N\) is completely infeasible. The solution therefore turns \(\gcd(n,n')\) into a multiplicative quantity and sums it by grouping integers according to their repeated prime factors.
Write
$$g(n)=\gcd(n,n').$$
The whole problem is to compute \(\sum_{n\le N} g(n)\) and then exclude \(n=1\).
Fix a prime \(p\), and write \(n=p^e m\) with \(p\nmid m\). Using the product rule,
$$n'=(p^e m)'=e\,p^{e-1}m+p^e m'=p^{e-1}(e\,m+p\,m').$$
Because \(m\) is not divisible by \(p\), the \(p\)-adic valuation of \(n'\) behaves in a simple way:
$$v_p(n')=\begin{cases} e-1,& p\nmid e,\\ \ge e,& p\mid e. \end{cases}$$
Therefore the \(p\)-part of \(\gcd(n,n')\) is
$$v_p(g(n))=\min\!\bigl(e,v_p(n')\bigr)=\begin{cases} e-1,& p\nmid e,\\ e,& p\mid e. \end{cases}$$
Equivalently, for one prime power,
$$g(p^e)=\begin{cases} p^{e-1},& p\nmid e,\\ p^e,& p\mid e. \end{cases}$$
Since each prime valuation is determined independently, \(g\) is multiplicative:
$$g(n)=\prod_{p^e\parallel n} g(p^e).$$
If a prime appears only once, then \(g(p)=1\). So primes of exponent \(1\) do not change \(g(n)\). This suggests writing every integer uniquely as
$$n=u\,s,$$
where \(u\) is the product of all prime powers \(p^e\) with \(e\ge 2\), and \(s\) is the product of the remaining primes of exponent \(1\).
Then:
$$u \text{ is powerful},\qquad s \text{ is squarefree},\qquad \gcd(s,u)=1,$$
and, crucially,
$$g(n)=g(u).$$
If we define the radical of \(u\) by
$$\operatorname{rad}(u)=\prod_{p\mid u} p,$$
then for each fixed powerful part \(u\), the admissible tails are exactly the squarefree integers \(s\le N/u\) that are coprime to \(\operatorname{rad}(u)\). Hence
$$S(N)+1=\sum_{\substack{u\le N\\ u\ \mathrm{powerful}}} g(u)\,Q\!\left(\left\lfloor\frac{N}{u}\right\rfloor,\operatorname{rad}(u)\right),$$
where
$$Q(x,r)=\#\{\,s\le x:\ s\text{ is squarefree and }\gcd(s,r)=1\,\}.$$
The extra \(1\) comes from \(n=1\), which is handled separately at the end.
The indicator of squarefreeness is
$$1_{\mathrm{sqf}}(s)=\sum_{d^2\mid s}\mu(d),$$
so Möbius inversion gives
$$Q(x,r)=\sum_{\substack{d\le \sqrt{x}\\ \gcd(d,r)=1}}\mu(d)\,\Phi_r\!\left(\left\lfloor\frac{x}{d^2}\right\rfloor\right),$$
where \(\Phi_r(y)\) counts integers up to \(y\) that are coprime to \(r\):
$$\Phi_r(y)=\#\{\,m\le y:\gcd(m,r)=1\,\}.$$
Because \(r=\operatorname{rad}(u)\) is squarefree, inclusion-exclusion over the primes dividing \(r\) gives
$$\Phi_r(y)=\sum_{t\mid r}\mu(t)\left\lfloor\frac{y}{t}\right\rfloor.$$
This is the exact counting formula used by the full implementation: one outer Möbius sum for squarefree numbers, and one inner divisor sum for the coprimality condition.
Every powerful number can be built by choosing distinct primes and assigning each chosen prime an exponent at least \(2\). A depth-first search over increasing primes therefore visits each powerful part \(u\le N\) exactly once.
At each node, the algorithm knows three pieces of information:
\(u\) itself, the already accumulated value \(g(u)\), and the prime set of \(\operatorname{rad}(u)\). With those data, it immediately adds
$$g(u)\,Q\!\left(\left\lfloor\frac{N}{u}\right\rfloor,\operatorname{rad}(u)\right)$$
to the answer, then extends \(u\) by appending a new prime with exponent \(2,3,4,\dots\) as long as the product stays within \(N\).
The shorter implementations encode the same idea in an equivalent telescoping recursion. If \(G_e(p)=g(p^e)\), then each increase
$$\Delta_e(p)=G_e(p)-G_{e-1}(p)$$
is multiplied by the number of integers whose powerful part first acquires the exponent \(e\) at that prime, and the recursion continues with larger primes only. This compact recursion and the explicit \(Q(x,r)\) formulation compute the same sum.
The powerful parts not exceeding \(20\) are
$$1,\ 4,\ 8,\ 9,\ 16.$$
Now count the squarefree tails for each one:
$$\begin{aligned} u=1&:&&g(u)=1,\quad Q(20,1)=13,&&\text{contribution }13,\\ u=4&:&&g(u)=4,\quad Q(5,2)=3,&&\text{contribution }12,\\ u=8&:&&g(u)=4,\quad Q(2,2)=1,&&\text{contribution }4,\\ u=9&:&&g(u)=3,\quad Q(2,3)=2,&&\text{contribution }6,\\ u=16&:&&g(u)=16,\quad Q(1,2)=1,&&\text{contribution }16. \end{aligned}$$
Therefore
$$\sum_{n=1}^{20} g(n)=13+12+4+6+16=51,$$
so after removing the \(n=1\) term we get
$$\sum_{n=2}^{20}\gcd(n,n')=50.$$
This small checkpoint matches the multiplicative formula perfectly.
The C++, Python, and Java implementations all begin from the same prime-power formula for \(g(p^e)\). The C++ implementation first validates that formula on small ranges, then sieves primes and Möbius values up to \(\sqrt{N}\). It enumerates powerful parts by depth-first search and, for each one, counts squarefree coprime tails with the Möbius-plus-inclusion-exclusion formula above. The same file also contains a more compact equivalent recursion, and that is the version mirrored by the Python and Java implementations: they recurse over increasing primes, let the exponent at one prime grow from \(2\) upward, and whenever the local factor of \(g\) increases they add that increment times the number of admissible remaining multiples. All three languages therefore compute the same decomposition, just with different levels of explicit counting infrastructure.
Let \(M=\lfloor\sqrt{N}\rfloor\). The sieve stage costs \(O(M\log\log M)\) time and \(O(M)\) memory. The search space of powerful numbers is sparse: there are far fewer powerful integers up to \(N\) than ordinary integers, so the recursion visits only a tiny fraction of the interval \([1,N]\). For each visited powerful part, the squarefree-coprime count sums over squarefree \(d\le \sqrt{N/u}\), while the inner inclusion-exclusion runs over divisors of \(\operatorname{rad}(u)\), whose prime support is small. In practice this is many orders of magnitude faster than iterating over every \(n\le N\), which is why \(N=5\times 10^{15}\) becomes feasible.
Gesucht ist
$$S(N)=\sum_{n=2}^{N}\gcd(n,n'),$$
wobei \(n'\) die arithmetische Ableitung ist. Sie wird durch \(p'=1\) für Primzahlen und die Produktregel \((ab)'=a'b+ab'\) definiert. Für
$$n=\prod p^{e_p}$$
folgt daraus
$$n'=n\sum_{p^{e_p}\parallel n}\frac{e_p}{p}.$$
Da \(N=5\times 10^{15}\) ist, scheidet direkte Enumeration vollständig aus. Die Lösung macht aus \(\gcd(n,n')\) eine multiplikative Größe und summiert dann nach der Struktur der wiederholten Primfaktoren.
Setze
$$g(n)=\gcd(n,n').$$
Dann ist die Aufgabe, \(\sum_{n\le N} g(n)\) zu berechnen und anschließend den Beitrag von \(n=1\) abzuziehen.
Fixiere eine Primzahl \(p\) und schreibe \(n=p^e m\) mit \(p\nmid m\). Aus der Produktregel erhält man
$$n'=(p^e m)'=e\,p^{e-1}m+p^e m'=p^{e-1}(e\,m+p\,m').$$
Weil \(m\) nicht durch \(p\) teilbar ist, ergibt sich für die \(p\)-adische Bewertung von \(n'\):
$$v_p(n')=\begin{cases} e-1,& p\nmid e,\\ \ge e,& p\mid e. \end{cases}$$
Damit ist der \(p\)-Anteil von \(\gcd(n,n')\)
$$v_p(g(n))=\min\!\bigl(e,v_p(n')\bigr)=\begin{cases} e-1,& p\nmid e,\\ e,& p\mid e. \end{cases}$$
Also gilt für eine einzelne Primzahlpotenz
$$g(p^e)=\begin{cases} p^{e-1},& p\nmid e,\\ p^e,& p\mid e. \end{cases}$$
Da jede Primzahl separat behandelt werden kann, ist \(g\) multiplikativ:
$$g(n)=\prod_{p^e\parallel n} g(p^e).$$
Tritt eine Primzahl nur mit Exponent \(1\) auf, dann ist \(g(p)=1\). Solche Primfaktoren beeinflussen \(g(n)\) also nicht. Deshalb schreibt man jedes \(n\) eindeutig als
$$n=u\,s,$$
wobei \(u\) alle Primzahlpotenzen \(p^e\) mit \(e\ge 2\) enthält und \(s\) aus den verbleibenden Primfaktoren mit Exponent \(1\) besteht.
Dann gilt:
$$u \text{ ist powerful},\qquad s \text{ ist quadratfrei},\qquad \gcd(s,u)=1,$$
und vor allem
$$g(n)=g(u).$$
Mit
$$\operatorname{rad}(u)=\prod_{p\mid u} p$$
werden für ein festes \(u\) genau die quadratfreien Zahlen \(s\le N/u\) benötigt, die zu \(\operatorname{rad}(u)\) teilerfremd sind. Daher
$$S(N)+1=\sum_{\substack{u\le N\\ u\ \mathrm{powerful}}} g(u)\,Q\!\left(\left\lfloor\frac{N}{u}\right\rfloor,\operatorname{rad}(u)\right),$$
wobei
$$Q(x,r)=\#\{\,s\le x:\ s\text{ quadratfrei und }\gcd(s,r)=1\,\}.$$
Das zusätzliche \(1\) entspricht dem Fall \(n=1\), der erst ganz am Ende entfernt wird.
Die Kennfunktion für Quadratfreiheit lautet
$$1_{\mathrm{sqf}}(s)=\sum_{d^2\mid s}\mu(d).$$
Damit folgt per Möbius-Inversion
$$Q(x,r)=\sum_{\substack{d\le \sqrt{x}\\ \gcd(d,r)=1}}\mu(d)\,\Phi_r\!\left(\left\lfloor\frac{x}{d^2}\right\rfloor\right),$$
wobei \(\Phi_r(y)\) die Anzahl der Zahlen \(\le y\) bezeichnet, die zu \(r\) teilerfremd sind:
$$\Phi_r(y)=\#\{\,m\le y:\gcd(m,r)=1\,\}.$$
Da \(r=\operatorname{rad}(u)\) quadratfrei ist, liefert Inklusion-Exklusion über seine Primteiler
$$\Phi_r(y)=\sum_{t\mid r}\mu(t)\left\lfloor\frac{y}{t}\right\rfloor.$$
Genau diese doppelte Struktur verwendet die ausführliche Implementierung: außen Möbius für Quadratfreiheit, innen Teiler-Summen für die Teilerfremdheit.
Jede powerful Zahl entsteht, indem man verschiedene Primzahlen wählt und jeder gewählten Primzahl einen Exponenten mindestens \(2\) gibt. Eine Tiefensuche über wachsende Primzahlen besucht deshalb jedes \(u\le N\) genau einmal.
Zu jedem Suchzustand kennt der Algorithmus \(u\), den schon aufgebauten Wert \(g(u)\) und die Primzahlmenge von \(\operatorname{rad}(u)\). Damit kann sofort
$$g(u)\,Q\!\left(\left\lfloor\frac{N}{u}\right\rfloor,\operatorname{rad}(u)\right)$$
zum Ergebnis addiert werden, bevor neue Primzahlpotenzen angehängt werden.
Die kürzeren Implementierungen schreiben dieselbe Rechnung als teleskopische Rekursion. Setzt man \(G_e(p)=g(p^e)\), dann wird jede Erhöhung
$$\Delta_e(p)=G_e(p)-G_{e-1}(p)$$
mit der Anzahl der noch möglichen Vielfachen multipliziert, und die Rekursion läuft nur mit größeren Primzahlen weiter. Inhaltlich ist das äquivalent zur expliziten Summation über powerful-Teile.
Die powerful-Teile bis \(20\) sind
$$1,\ 4,\ 8,\ 9,\ 16.$$
Für jeden davon zählt man die zulässigen quadratfreien Reste:
$$\begin{aligned} u=1&:&&g(u)=1,\quad Q(20,1)=13,&&\text{Beitrag }13,\\ u=4&:&&g(u)=4,\quad Q(5,2)=3,&&\text{Beitrag }12,\\ u=8&:&&g(u)=4,\quad Q(2,2)=1,&&\text{Beitrag }4,\\ u=9&:&&g(u)=3,\quad Q(2,3)=2,&&\text{Beitrag }6,\\ u=16&:&&g(u)=16,\quad Q(1,2)=1,&&\text{Beitrag }16. \end{aligned}$$
Also
$$\sum_{n=1}^{20} g(n)=13+12+4+6+16=51,$$
und nach Abzug von \(n=1\) bleibt
$$\sum_{n=2}^{20}\gcd(n,n')=50.$$
Die C++-, Python- und Java-Implementierungen gehen alle von derselben Primzahlpotenzformel für \(g(p^e)\) aus. Die C++-Implementierung prüft diese Formel zunächst auf kleinen Bereichen, siebt dann Primzahlen und Möbius-Werte bis \(\sqrt{N}\) und enumeriert alle powerful-Teile per Tiefensuche. Für jeden Zustand wird die Zahl der quadratfreien, zu \(\operatorname{rad}(u)\) teilerfremden Reste mit der oben hergeleiteten Möbius- plus Inklusion-Exklusion-Formel berechnet. Dieselbe Datei enthält außerdem eine kompakte äquivalente Rekursion; genau diese kürzere Sicht wird in Python und Java umgesetzt. Dort wächst der Exponent einer Primzahl stufenweise, und immer wenn sich der lokale ggT-Faktor erhöht, wird dieser Zuwachs mit der Zahl der verbleibenden zulässigen Vielfachen und der Rekursion über größere Primzahlen multipliziert. Alle drei Sprachen berechnen also dieselbe Summenzerlegung.
Mit \(M=\lfloor\sqrt{N}\rfloor\) kostet das Sieb \(O(M\log\log M)\) Zeit und \(O(M)\) Speicher. Die Zahl der powerful-Teile bis \(N\) ist sehr klein im Vergleich zu \(N\) selbst, daher betrachtet die Rekursion nur einen dünnen Suchraum. Für jedes besuchte \(u\) summiert die Zählfunktion über quadratfreie \(d\le \sqrt{N/u}\); die innere Inklusion-Exklusion läuft nur über Teiler von \(\operatorname{rad}(u)\), also über eine kleine Primzahlmenge. Insgesamt ist das Verfahren um Größenordnungen schneller als ein vollständiger Durchlauf über alle \(n\le N\) und damit für \(N=5\times 10^{15}\) praktikabel.
Hesaplanması gereken toplam şudur:
$$S(N)=\sum_{n=2}^{N}\gcd(n,n').$$
Burada \(n'\) aritmetik türevdir; her asal için \(p'=1\) alınır ve çarpım kuralı \((ab)'=a'b+ab'\) kullanılır. Eğer
$$n=\prod p^{e_p}$$
ise, bundan
$$n'=n\sum_{p^{e_p}\parallel n}\frac{e_p}{p}$$
sonucu çıkar. Gerçek hedefte \(N=5\times 10^{15}\) olduğundan, tüm sayıları tek tek incelemek mümkün değildir. Çözüm, \(\gcd(n,n')\) ifadesini çarpımsal bir fonksiyona dönüştürüp toplamı tekrarlı asal çarpan yapısına göre gruplar.
Şunu tanımlayalım:
$$g(n)=\gcd(n,n').$$
Amaç önce \(\sum_{n\le N} g(n)\) toplamını bulmak, sonra da \(n=1\) terimini çıkarmaktır.
Bir asal \(p\) seçelim ve \(n=p^e m\), \(p\nmid m\) yazalım. Çarpım kuralı ile
$$n'=(p^e m)'=e\,p^{e-1}m+p^e m'=p^{e-1}(e\,m+p\,m')$$
elde edilir. \(m\), \(p\) ile bölünmediği için \(n'\) içindeki \(p\)-üssü şu şekilde davranır:
$$v_p(n')=\begin{cases} e-1,& p\nmid e,\\ \ge e,& p\mid e. \end{cases}$$
Dolayısıyla \(\gcd(n,n')\) içindeki \(p\)-katkısı
$$v_p(g(n))=\min\!\bigl(e,v_p(n')\bigr)=\begin{cases} e-1,& p\nmid e,\\ e,& p\mid e. \end{cases}$$
olur. Yani tek bir asal kuvvet için
$$g(p^e)=\begin{cases} p^{e-1},& p\nmid e,\\ p^e,& p\mid e. \end{cases}$$
Her asal ayrı ayrı davrandığından \(g\) çarpımsaldır:
$$g(n)=\prod_{p^e\parallel n} g(p^e).$$
Bir asalın üssü \(1\) ise \(g(p)=1\) olur. Demek ki üsleri yalnızca \(1\) olan asallar, \(g(n)\)'yi değiştirmez. Bu yüzden her sayı tekil biçimde
$$n=u\,s$$
şeklinde yazılır; burada \(u\), üssü en az \(2\) olan tüm asal kuvvetleri içerir ve \(s\), üssü \(1\) olan geri kalan asalların çarpımıdır.
Böylece
$$u \text{ powerful'dır},\qquad s \text{ karekök-özgürdür},\qquad \gcd(s,u)=1,$$
ve en önemli nokta
$$g(n)=g(u)$$
olur. Ayrıca
$$\operatorname{rad}(u)=\prod_{p\mid u} p$$
tanımıyla, sabit bir \(u\) için izin verilen kuyruklar tam olarak \(\operatorname{rad}(u)\) ile aralarında asal olan ve \(N/u\)'yu aşmayan karekök-özgür sayılardır. Bu nedenle
$$S(N)+1=\sum_{\substack{u\le N\\ u\ \mathrm{powerful}}} g(u)\,Q\!\left(\left\lfloor\frac{N}{u}\right\rfloor,\operatorname{rad}(u)\right),$$
burada
$$Q(x,r)=\#\{\,s\le x:\ s\text{ karekök-özgür ve }\gcd(s,r)=1\,\}.$$
Eklenen \(1\), sonunda ayrı çıkarılan \(n=1\) teriminden gelir.
Karekök-özgür olma göstergesi
$$1_{\mathrm{sqf}}(s)=\sum_{d^2\mid s}\mu(d)$$
olduğundan Möbius terslemesi ile
$$Q(x,r)=\sum_{\substack{d\le \sqrt{x}\\ \gcd(d,r)=1}}\mu(d)\,\Phi_r\!\left(\left\lfloor\frac{x}{d^2}\right\rfloor\right)$$
elde edilir. Burada \(\Phi_r(y)\), \(y\)'ye kadar olan ve \(r\) ile aralarında asal olan sayıların sayısıdır:
$$\Phi_r(y)=\#\{\,m\le y:\gcd(m,r)=1\,\}.$$
\(r=\operatorname{rad}(u)\) karekök-özgür olduğu için, \(r\)'yi bölen asallar üzerinde dahil et-çıkar uygulanır:
$$\Phi_r(y)=\sum_{t\mid r}\mu(t)\left\lfloor\frac{y}{t}\right\rfloor.$$
Geniş C++ uygulaması tam olarak bu iki katmanlı sayımı kullanır: dışta karekök-özgürlük için Möbius, içte aralarında asallık için bölen toplamı.
Her powerful sayı, farklı asallar seçilip her birine en az \(2\) üssü verilerek kurulur. Artan asallar üzerinde yapılan bir DFS bu yüzden her \(u\le N\) değerini tam bir kez ziyaret eder.
Her düğümde algoritma \(u\)'yu, o ana kadar biriken \(g(u)\) değerini ve \(\operatorname{rad}(u)\)'nun asal kümesini bilir. Böylece hemen
$$g(u)\,Q\!\left(\left\lfloor\frac{N}{u}\right\rfloor,\operatorname{rad}(u)\right)$$
katkısı eklenir; sonra yeni bir asal, üsleri \(2,3,4,\dots\) olacak şekilde eklenerek arama derinleşir.
Daha kısa uygulamalar aynı hesabı teleskopik bir özyineleme olarak yazar. \(G_e(p)=g(p^e)\) denirse, her artış
$$\Delta_e(p)=G_e(p)-G_{e-1}(p)$$
kalan izinli katların sayısıyla çarpılır ve özyineleme yalnız daha büyük asallarla devam eder. Yani kısa sürüm ile açık \(Q(x,r)\) formülü aynı toplamı hesaplar.
\(20\)'yi aşmayan powerful kısımlar şunlardır:
$$1,\ 4,\ 8,\ 9,\ 16.$$
Her biri için uygun karekök-özgür kuyruklar sayılır:
$$\begin{aligned} u=1&:&&g(u)=1,\quad Q(20,1)=13,&&\text{katkı }13,\\ u=4&:&&g(u)=4,\quad Q(5,2)=3,&&\text{katkı }12,\\ u=8&:&&g(u)=4,\quad Q(2,2)=1,&&\text{katkı }4,\\ u=9&:&&g(u)=3,\quad Q(2,3)=2,&&\text{katkı }6,\\ u=16&:&&g(u)=16,\quad Q(1,2)=1,&&\text{katkı }16. \end{aligned}$$
Böylece
$$\sum_{n=1}^{20} g(n)=13+12+4+6+16=51$$
ve \(n=1\) çıkarılınca
$$\sum_{n=2}^{20}\gcd(n,n')=50$$
elde edilir.
C++, Python ve Java uygulamalarının çıkış noktası aynıdır: \(g(p^e)\) için bulunan asal-kuvvet formülü. C++ uygulaması önce bu formülü küçük aralıklarda doğrular, sonra \(\sqrt{N}\)'ye kadar asal ve Möbius değerlerini üretir. Powerful kısımlar DFS ile gezilir; her durumda karekök-özgür ve \(\operatorname{rad}(u)\) ile aralarında asal kuyrukların sayısı yukarıdaki Möbius artı dahil et-çıkar formülüyle hesaplanır. Aynı dosyada bunun daha sıkı yazılmış eşdeğer bir özyinelemeli biçimi de vardır; Python ve Java sürümleri o kompakt yaklaşımı izler. Orada bir asalın üssü adım adım artırılır, yerel \(\gcd\) katsayısı ne zaman büyürse bu artış kalan izinli çarpan sayısıyla ve daha büyük asallar üzerindeki özyineleme ile çarpılarak toplanır. Üç dil de aynı matematiksel ayrışımı hesaplar.
\(M=\lfloor\sqrt{N}\rfloor\) olsun. Elek aşaması \(O(M\log\log M)\) zaman ve \(O(M)\) bellek kullanır. Powerful sayıların sayısı \(N\)'ye göre çok seyrek olduğundan, arama yalnız ince bir alt kümede dolaşır. Her ziyaret edilen \(u\) için kuyruk sayımı, \(\sqrt{N/u}\)'ya kadar karekök-özgür \(d\) değerleri üzerinde döner; iç dahil et-çıkar ise \(\operatorname{rad}(u)\)'nun bölenleriyle sınırlı olduğundan küçüktür. Sonuç olarak yöntem, tüm \(n\le N\) değerlerini taramaya kıyasla dramatik biçimde daha hızlıdır ve \(N=5\times 10^{15}\) için uygulanabilir hale gelir.
Debemos calcular
$$S(N)=\sum_{n=2}^{N}\gcd(n,n').$$
Aquí \(n'\) es la derivada aritmética, definida por \(p'=1\) para todo primo \(p\) y por la regla del producto \((ab)'=a'b+ab'\). Si
$$n=\prod p^{e_p},$$
entonces
$$n'=n\sum_{p^{e_p}\parallel n}\frac{e_p}{p}.$$
Como el valor objetivo es \(N=5\times 10^{15}\), no se puede factorizar cada entero hasta \(N\). La idea es convertir \(\gcd(n,n')\) en una función multiplicativa y reagrupar la suma según la parte con exponentes repetidos.
Definamos
$$g(n)=\gcd(n,n').$$
La tarea es calcular \(\sum_{n\le N} g(n)\) y después retirar el caso \(n=1\).
Fijemos un primo \(p\) y escribamos \(n=p^e m\) con \(p\nmid m\). Por la regla del producto,
$$n'=(p^e m)'=e\,p^{e-1}m+p^e m'=p^{e-1}(e\,m+p\,m').$$
Como \(m\) no es divisible por \(p\), la valuación \(p\)-ádica de \(n'\) satisface
$$v_p(n')=\begin{cases} e-1,& p\nmid e,\\ \ge e,& p\mid e. \end{cases}$$
Por tanto, la contribución de \(p\) a \(\gcd(n,n')\) es
$$v_p(g(n))=\min\!\bigl(e,v_p(n')\bigr)=\begin{cases} e-1,& p\nmid e,\\ e,& p\mid e. \end{cases}$$
Es decir, para una sola potencia prima,
$$g(p^e)=\begin{cases} p^{e-1},& p\nmid e,\\ p^e,& p\mid e. \end{cases}$$
Como esto se decide primo por primo, \(g\) es multiplicativa:
$$g(n)=\prod_{p^e\parallel n} g(p^e).$$
Si un primo aparece con exponente \(1\), entonces \(g(p)=1\). Por eso los factores con exponente \(1\) no cambian el valor de \(g(n)\). Toda \(n\) puede escribirse de forma única como
$$n=u\,s,$$
donde \(u\) reúne todas las potencias primas con exponente al menos \(2\), y \(s\) es el producto de los primos restantes con exponente \(1\).
Así se obtiene
$$u \text{ es powerful},\qquad s \text{ es libre de cuadrados},\qquad \gcd(s,u)=1,$$
y además
$$g(n)=g(u).$$
Si definimos
$$\operatorname{rad}(u)=\prod_{p\mid u} p,$$
entonces, para cada \(u\) fijo, las colas admisibles son exactamente los \(s\le N/u\) que son libres de cuadrados y coprimos con \(\operatorname{rad}(u)\). Por lo tanto,
$$S(N)+1=\sum_{\substack{u\le N\\ u\ \mathrm{powerful}}} g(u)\,Q\!\left(\left\lfloor\frac{N}{u}\right\rfloor,\operatorname{rad}(u)\right),$$
con
$$Q(x,r)=\#\{\,s\le x:\ s\text{ es libre de cuadrados y }\gcd(s,r)=1\,\}.$$
El término adicional \(1\) corresponde a \(n=1\), que se elimina al final.
La función indicadora de los enteros libres de cuadrados es
$$1_{\mathrm{sqf}}(s)=\sum_{d^2\mid s}\mu(d).$$
Aplicando inversión de Möbius,
$$Q(x,r)=\sum_{\substack{d\le \sqrt{x}\\ \gcd(d,r)=1}}\mu(d)\,\Phi_r\!\left(\left\lfloor\frac{x}{d^2}\right\rfloor\right),$$
donde \(\Phi_r(y)\) cuenta los enteros \(\le y\) coprimos con \(r\):
$$\Phi_r(y)=\#\{\,m\le y:\gcd(m,r)=1\,\}.$$
Como \(r=\operatorname{rad}(u)\) es libre de cuadrados, la inclusión-exclusión sobre sus primos produce
$$\Phi_r(y)=\sum_{t\mid r}\mu(t)\left\lfloor\frac{y}{t}\right\rfloor.$$
Esa es exactamente la estructura del conteo completo: una suma exterior de Möbius para imponer libertad de cuadrados y una suma interior sobre divisores para imponer coprimalidad.
Todo número powerful se construye eligiendo primos distintos y asignando a cada uno un exponente al menos \(2\). Una búsqueda en profundidad sobre primos crecientes enumera así cada \(u\le N\) exactamente una vez.
En cada estado, el algoritmo conoce \(u\), el valor acumulado \(g(u)\) y el conjunto de primos de \(\operatorname{rad}(u)\). Entonces añade de inmediato
$$g(u)\,Q\!\left(\left\lfloor\frac{N}{u}\right\rfloor,\operatorname{rad}(u)\right)$$
y luego prueba a incorporar un nuevo primo con exponente \(2,3,4,\dots\), mientras el producto siga siendo \(\le N\).
Las implementaciones más cortas usan una recursión telescópica equivalente. Si \(G_e(p)=g(p^e)\), cada incremento
$$\Delta_e(p)=G_e(p)-G_{e-1}(p)$$
se multiplica por la cantidad de múltiplos todavía posibles, y la recursión continúa solo con primos mayores. Esa versión compacta y la suma explícita sobre \(Q(x,r)\) son dos formas de la misma descomposición.
Las partes powerful no mayores que \(20\) son
$$1,\ 4,\ 8,\ 9,\ 16.$$
Para cada una contamos las colas cuadrado-libres admisibles:
$$\begin{aligned} u=1&:&&g(u)=1,\quad Q(20,1)=13,&&\text{aporte }13,\\ u=4&:&&g(u)=4,\quad Q(5,2)=3,&&\text{aporte }12,\\ u=8&:&&g(u)=4,\quad Q(2,2)=1,&&\text{aporte }4,\\ u=9&:&&g(u)=3,\quad Q(2,3)=2,&&\text{aporte }6,\\ u=16&:&&g(u)=16,\quad Q(1,2)=1,&&\text{aporte }16. \end{aligned}$$
Así,
$$\sum_{n=1}^{20} g(n)=13+12+4+6+16=51,$$
y quitando \(n=1\) obtenemos
$$\sum_{n=2}^{20}\gcd(n,n')=50.$$
Las implementaciones en C++, Python y Java parten de la misma fórmula local para \(g(p^e)\). La implementación en C++ primero verifica esa fórmula en rangos pequeños, luego genera los primos y los valores de Möbius hasta \(\sqrt{N}\). Después recorre en profundidad todas las partes powerful y, para cada una, evalúa el número de colas libres de cuadrados y coprimas mediante la fórmula de Möbius más inclusión-exclusión descrita arriba. El mismo archivo contiene además una recursión equivalente mucho más compacta, y esa es la versión que reflejan Python y Java: avanzan primo por primo, dejan crecer el exponente desde \(2\) en adelante, y cada vez que aumenta el factor local de \(g\) suman ese incremento multiplicado por el número de múltiplos todavía admisibles y por la continuación recursiva con primos mayores. Las tres implementaciones calculan exactamente la misma suma.
Si \(M=\lfloor\sqrt{N}\rfloor\), la criba cuesta \(O(M\log\log M)\) tiempo y \(O(M)\) memoria. El conjunto de números powerful es muy ralo, así que la búsqueda visita una fracción diminuta de \([1,N]\). Para cada \(u\) visitado, el conteo de la cola suma sobre valores cuadrado-libres \(d\le \sqrt{N/u}\); la inclusión-exclusión interior solo recorre divisores de \(\operatorname{rad}(u)\), cuyo soporte primo es pequeño. En la práctica esto es muchísimo más rápido que procesar uno por uno todos los enteros hasta \(N\), y por eso \(N=5\times 10^{15}\) resulta alcanzable.
要求计算
$$S(N)=\sum_{n=2}^{N}\gcd(n,n').$$
其中 \(n'\) 是算术导数。它由两条规则定义：对任意素数 \(p\)，有 \(p'=1\)；同时满足乘法求导法则 \((ab)'=a'b+ab'\)。因此若
$$n=\prod p^{e_p},$$
就有
$$n'=n\sum_{p^{e_p}\parallel n}\frac{e_p}{p}.$$
本题的目标规模是 \(N=5\times 10^{15}\)。如果逐个整数分解质因数并直接求 \(\gcd(n,n')\)，时间上根本不可行。真正可行的思路是先把 \(\gcd(n,n')\) 化成一个按质因数分离的乘法函数，再按照“重复质因子部分”和“平方自由尾部”来重组总和。
记
$$g(n)=\gcd(n,n').$$
那么问题就是先求 \(\sum_{n\le N} g(n)\)，最后再把 \(n=1\) 的那一项去掉。
固定一个素数 \(p\)，把 \(n\) 写成 \(n=p^e m\)，其中 \(p\nmid m\)。由乘法法则可得
$$n'=(p^e m)'=e\,p^{e-1}m+p^e m'=p^{e-1}(e\,m+p\,m').$$
因为 \(m\) 不被 \(p\) 整除，所以 \(n'\) 的 \(p\)-进赋值非常简单：
$$v_p(n')=\begin{cases} e-1,& p\nmid e,\\ \ge e,& p\mid e. \end{cases}$$
因此 \(\gcd(n,n')\) 中关于 \(p\) 的指数就是
$$v_p(g(n))=\min\!\bigl(e,v_p(n')\bigr)=\begin{cases} e-1,& p\nmid e,\\ e,& p\mid e. \end{cases}$$
也就是说，对单个素数幂有
$$g(p^e)=\begin{cases} p^{e-1},& p\nmid e,\\ p^e,& p\mid e. \end{cases}$$
每个素数的贡献彼此独立，所以 \(g\) 是乘法函数：
$$g(n)=\prod_{p^e\parallel n} g(p^e).$$
如果某个素数只出现一次，那么 \(g(p)=1\)。这说明指数恰好为 \(1\) 的素因子不会改变 \(g(n)\)。于是每个正整数都可以唯一地写成
$$n=u\,s,$$
其中 \(u\) 收集了所有指数至少为 \(2\) 的素数幂，而 \(s\) 则是其余指数为 \(1\) 的素数之积。
于是就有
$$u \text{ 是 powerful 数},\qquad s \text{ 是平方自由数},\qquad \gcd(s,u)=1,$$
并且最关键的是
$$g(n)=g(u).$$
再定义
$$\operatorname{rad}(u)=\prod_{p\mid u} p,$$
那么对于固定的 \(u\)，所有允许的尾部 \(s\) 恰好是满足 \(s\le N/u\)、平方自由、并且与 \(\operatorname{rad}(u)\) 互素的那些整数。因此
$$S(N)+1=\sum_{\substack{u\le N\\ u\ \mathrm{powerful}}} g(u)\,Q\!\left(\left\lfloor\frac{N}{u}\right\rfloor,\operatorname{rad}(u)\right),$$
其中
$$Q(x,r)=\#\{\,s\le x:\ s\text{ 为平方自由且 }\gcd(s,r)=1\,\}.$$
这里多出来的 \(1\) 对应 \(n=1\)，实际答案最后再把它减掉。
平方自由指示函数满足
$$1_{\mathrm{sqf}}(s)=\sum_{d^2\mid s}\mu(d).$$
因此用 Möbius 反演可以得到
$$Q(x,r)=\sum_{\substack{d\le \sqrt{x}\\ \gcd(d,r)=1}}\mu(d)\,\Phi_r\!\left(\left\lfloor\frac{x}{d^2}\right\rfloor\right),$$
其中 \(\Phi_r(y)\) 表示不超过 \(y\) 且与 \(r\) 互素的整数个数：
$$\Phi_r(y)=\#\{\,m\le y:\gcd(m,r)=1\,\}.$$
由于这里的 \(r=\operatorname{rad}(u)\) 本身就是平方自由的，所以对 \(r\) 的所有素因子做容斥即可：
$$\Phi_r(y)=\sum_{t\mid r}\mu(t)\left\lfloor\frac{y}{t}\right\rfloor.$$
这正是完整实现中出现的两层计数结构：外层 Möbius 和用于强制平方自由的 \(d^2\) 约束，内层是对 \(\operatorname{rad}(u)\) 的因子做容斥，从而保证互素条件。
任意一个 powerful 数都可以通过“选若干个不同素数，并给每个素数至少赋指数 \(2\)”来构造。因此只要按素数递增做深度优先搜索，就能把每个 \(u\le N\) 恰好访问一次。
在搜索的每个节点，算法都掌握三类信息：当前的 \(u\)、已经累积好的 \(g(u)\)，以及 \(\operatorname{rad}(u)\) 对应的素数集合。于是该节点的贡献立刻就是
$$g(u)\,Q\!\left(\left\lfloor\frac{N}{u}\right\rfloor,\operatorname{rad}(u)\right).$$
然后继续尝试加入一个新的素数，并把它的指数依次取为 \(2,3,4,\dots\)，只要乘积仍然不超过 \(N\) 就继续递归。
较短的实现把同一件事写成了一个等价的“望远镜式”递归。如果记 \(G_e(p)=g(p^e)\)，那么每当指数从 \(e-1\) 增长到 \(e\) 时，就会产生一个增量
$$\Delta_e(p)=G_e(p)-G_{e-1}(p).$$
这个增量乘上“还能放进剩余范围中的倍数个数”，再乘上更大素数的递归贡献，就得到和显式 powerful 分解完全相同的总和。Python 和 Java 版本采用的正是这种紧凑写法。
不超过 \(20\) 的 powerful 部分只有
$$1,\ 4,\ 8,\ 9,\ 16.$$
分别计算对应的平方自由尾部个数：
$$\begin{aligned} u=1&:&&g(u)=1,\quad Q(20,1)=13,&&\text{贡献 }13,\\ u=4&:&&g(u)=4,\quad Q(5,2)=3,&&\text{贡献 }12,\\ u=8&:&&g(u)=4,\quad Q(2,2)=1,&&\text{贡献 }4,\\ u=9&:&&g(u)=3,\quad Q(2,3)=2,&&\text{贡献 }6,\\ u=16&:&&g(u)=16,\quad Q(1,2)=1,&&\text{贡献 }16. \end{aligned}$$
因此
$$\sum_{n=1}^{20} g(n)=13+12+4+6+16=51,$$
去掉 \(n=1\) 后得到
$$\sum_{n=2}^{20}\gcd(n,n')=50.$$
这个小例子清楚地展示了“先枚举 powerful 部分，再统计平方自由尾部”的总和结构。
C++、Python 和 Java 实现都建立在同一个局部公式 \(g(p^e)\) 上。C++ 实现先在小范围内验证这个公式，然后筛出 \(\sqrt{N}\) 以内的素数和 Möbius 值，再通过深度优先搜索枚举所有 possible powerful 部分。对每个枚举到的部分，都用上面的 Möbius 反演加容斥公式计算平方自由且互素的尾部数量。C++ 文件中还包含一个更紧凑但完全等价的递归写法；Python 和 Java 采用的就是这个版本。它们按素数递增递归，逐步提升某个素数的指数，只要局部 \(\gcd\) 因子发生增加，就把这个增量乘以剩余可用倍数的数量，再乘上更大素数的后续贡献。三种语言只是写法不同，计算的数学对象完全相同。
设 \(M=\lfloor\sqrt{N}\rfloor\)。筛法阶段需要 \(O(M\log\log M)\) 时间和 \(O(M)\) 空间。真正被枚举的 powerful 数远少于全部整数，因此递归搜索只访问很稀疏的一部分状态。对每个访问到的 \(u\)，尾部计数需要遍历不超过 \(\sqrt{N/u}\) 的平方自由 \(d\)，而内层容斥只在 \(\operatorname{rad}(u)\) 的因子上运行，这个集合通常很小。所以总体运行量与直接遍历 \(1\) 到 \(N\) 相比小了几个数量级，才使得 \(N=5\times 10^{15}\) 变得可计算。
Нужно вычислить сумму
$$S(N)=\sum_{n=2}^{N}\gcd(n,n').$$
Здесь \(n'\) обозначает арифметическую производную. Она задается правилами \(p'=1\) для любого простого \(p\) и \((ab)'=a'b+ab'\). Поэтому при
$$n=\prod p^{e_p}$$
получаем
$$n'=n\sum_{p^{e_p}\parallel n}\frac{e_p}{p}.$$
Для \(N=5\times 10^{15}\) полный перебор невозможен. Решение сводит \(\gcd(n,n')\) к мультипликативной функции и затем суммирует вклад не по всем числам подряд, а по их powerful-части и квадратсвободному хвосту.
Обозначим
$$g(n)=\gcd(n,n').$$
Тогда нужно сначала найти \(\sum_{n\le N} g(n)\), а затем исключить член \(n=1\).
Зафиксируем простое \(p\) и представим \(n\) в виде \(n=p^e m\), где \(p\nmid m\). По правилу произведения
$$n'=(p^e m)'=e\,p^{e-1}m+p^e m'=p^{e-1}(e\,m+p\,m').$$
Так как \(m\) не делится на \(p\), \(p\)-адическая оценка \(n'\) имеет вид
$$v_p(n')=\begin{cases} e-1,& p\nmid e,\\ \ge e,& p\mid e. \end{cases}$$
Отсюда показатель \(p\) в \(\gcd(n,n')\) равен
$$v_p(g(n))=\min\!\bigl(e,v_p(n')\bigr)=\begin{cases} e-1,& p\nmid e,\\ e,& p\mid e. \end{cases}$$
Следовательно, для одной простой степени
$$g(p^e)=\begin{cases} p^{e-1},& p\nmid e,\\ p^e,& p\mid e. \end{cases}$$
Поскольку эта формула определяется независимо по каждому простому, функция \(g\) мультипликативна:
$$g(n)=\prod_{p^e\parallel n} g(p^e).$$
Если простой множитель входит с показателем \(1\), то \(g(p)=1\). Значит, такие множители вообще не меняют значение \(g(n)\). Поэтому каждое число единственным образом раскладывается как
$$n=u\,s,$$
где \(u\) содержит все простые степени с показателем не меньше \(2\), а \(s\) является произведением оставшихся простых с показателем \(1\).
Тогда
$$u \text{ powerful},\qquad s \text{ квадратсвободно},\qquad \gcd(s,u)=1,$$
и, что важнее всего,
$$g(n)=g(u).$$
Если обозначить
$$\operatorname{rad}(u)=\prod_{p\mid u} p,$$
то при фиксированном \(u\) допустимые хвосты \(s\) — это ровно квадратсвободные числа \(s\le N/u\), взаимно простые с \(\operatorname{rad}(u)\). Значит,
$$S(N)+1=\sum_{\substack{u\le N\\ u\ \mathrm{powerful}}} g(u)\,Q\!\left(\left\lfloor\frac{N}{u}\right\rfloor,\operatorname{rad}(u)\right),$$
где
$$Q(x,r)=\#\{\,s\le x:\ s\text{ квадратсвободно и }\gcd(s,r)=1\,\}.$$
Лишняя единица соответствует \(n=1\), который вычитается в самом конце.
Индикатор квадратсвободности равен
$$1_{\mathrm{sqf}}(s)=\sum_{d^2\mid s}\mu(d).$$
Отсюда по формуле Мёбиуса следует
$$Q(x,r)=\sum_{\substack{d\le \sqrt{x}\\ \gcd(d,r)=1}}\mu(d)\,\Phi_r\!\left(\left\lfloor\frac{x}{d^2}\right\rfloor\right),$$
где \(\Phi_r(y)\) — число целых \(m\le y\), взаимно простых с \(r\):
$$\Phi_r(y)=\#\{\,m\le y:\gcd(m,r)=1\,\}.$$
Так как \(r=\operatorname{rad}(u)\) квадратсвободно, внутренний счет выполняется обычным включением-исключением по простым делителям \(r\):
$$\Phi_r(y)=\sum_{t\mid r}\mu(t)\left\lfloor\frac{y}{t}\right\rfloor.$$
Именно эта двухуровневая структура и используется в полном решении: внешняя сумма Мёбиуса навязывает квадратсвободность, а внутренняя сумма по делителям отвечает за условие взаимной простоты.
Любое powerful-число получается выбором различных простых и назначением каждому выбранному простому показателя не меньше \(2\). Поэтому поиск в глубину по возрастающим простым перечисляет каждое \(u\le N\) ровно один раз.
В каждом состоянии алгоритм знает текущее \(u\), уже накопленное значение \(g(u)\) и множество простых в \(\operatorname{rad}(u)\). Значит, можно сразу добавить
$$g(u)\,Q\!\left(\left\lfloor\frac{N}{u}\right\rfloor,\operatorname{rad}(u)\right).$$
После этого к \(u\) присоединяется новый простой множитель со степенью \(2,3,4,\dots\), пока произведение не превышает \(N\).
Короткие реализации записывают ту же идею в виде телескопической рекурсии. Если обозначить \(G_e(p)=g(p^e)\), то каждый прирост
$$\Delta_e(p)=G_e(p)-G_{e-1}(p)$$
умножается на количество допустимых оставшихся кратных, а затем рекурсия продолжается только по большим простым. Это полностью эквивалентно явной сумме по powerful-частям.
Powerful-части, не превосходящие \(20\), таковы:
$$1,\ 4,\ 8,\ 9,\ 16.$$
Для каждой из них считаем допустимые квадратсвободные хвосты:
$$\begin{aligned} u=1&:&&g(u)=1,\quad Q(20,1)=13,&&\text{вклад }13,\\ u=4&:&&g(u)=4,\quad Q(5,2)=3,&&\text{вклад }12,\\ u=8&:&&g(u)=4,\quad Q(2,2)=1,&&\text{вклад }4,\\ u=9&:&&g(u)=3,\quad Q(2,3)=2,&&\text{вклад }6,\\ u=16&:&&g(u)=16,\quad Q(1,2)=1,&&\text{вклад }16. \end{aligned}$$
Получаем
$$\sum_{n=1}^{20} g(n)=13+12+4+6+16=51,$$
а после исключения \(n=1\)
$$\sum_{n=2}^{20}\gcd(n,n')=50.$$
Реализации на C++, Python и Java исходят из одной и той же локальной формулы для \(g(p^e)\). В версии на C++ эта формула сначала проверяется на малых диапазонах, затем строится решето простых чисел и значений функции Мёбиуса до \(\sqrt{N}\). После этого алгоритм в глубину перечисляет powerful-части и для каждой из них вычисляет число квадратсвободных хвостов, взаимно простых с \(\operatorname{rad}(u)\), с помощью приведенной выше формулы Мёбиуса и включения-исключения. В том же файле есть и более компактная, но эквивалентная рекурсия; именно она отражена в версиях на Python и Java. Там показатель у текущего простого последовательно растет, и каждый раз, когда локальный множитель \(\gcd\) увеличивается, этот прирост умножается на число еще допустимых кратных и на рекурсивный вклад больших простых. Все три реализации вычисляют одну и ту же математическую сумму.
Пусть \(M=\lfloor\sqrt{N}\rfloor\). Построение решета требует \(O(M\log\log M)\) времени и \(O(M)\) памяти. Количество powerful-чисел до \(N\) крайне мало по сравнению с \(N\), поэтому рекурсивный перебор посещает только редкое множество состояний. Для каждого посещенного \(u\) подсчет хвоста суммирует по квадратсвободным \(d\le \sqrt{N/u}\), а внутренняя часть включения-исключения идет лишь по делителям \(\operatorname{rad}(u)\), то есть по небольшому набору простых. Благодаря этому алгоритм работает на многие порядки быстрее полного перебора всех \(n\le N\), что и делает возможным случай \(N=5\times 10^{15}\).
المطلوب هو حساب
$$S(N)=\sum_{n=2}^{N}\gcd(n,n').$$
حيث إن \(n'\) هو الاشتقاق الحسابي. يُعرَّف هذا الاشتقاق بشرطين: \(p'=1\) لكل عدد أولي \(p\)، وقاعدة الضرب \((ab)'=a'b+ab'\). ولذلك إذا كان
$$n=\prod p^{e_p},$$
فإن
$$n'=n\sum_{p^{e_p}\parallel n}\frac{e_p}{p}.$$
بما أن القيمة المطلوبة تستخدم \(N=5\times 10^{15}\)، فلا يمكن فحص كل عدد حتى \(N\) مباشرة. الفكرة الصحيحة هي تحويل \(\gcd(n,n')\) إلى دالة مضاعِفية ثم إعادة ترتيب الجمع وفق الجزء الذي يحتوي على الأسس المتكررة والذيل الخالي من المربعات.
لنكتب
$$g(n)=\gcd(n,n').$$
وعندئذ تصبح المهمة هي حساب \(\sum_{n\le N} g(n)\) ثم حذف حد \(n=1\) في النهاية.
ثبّت عددًا أوليًا \(p\)، واكتب \(n=p^e m\) مع \(p\nmid m\). من قاعدة الضرب نحصل على
$$n'=(p^e m)'=e\,p^{e-1}m+p^e m'=p^{e-1}(e\,m+p\,m').$$
وبما أن \(m\) غير قابل للقسمة على \(p\)، فإن تقييم \(p\)-الأدي لـ \(n'\) يأخذ الصورة
$$v_p(n')=\begin{cases} e-1,& p\nmid e,\\ \ge e,& p\mid e. \end{cases}$$
إذًا أس \(p\) داخل \(\gcd(n,n')\) يساوي
$$v_p(g(n))=\min\!\bigl(e,v_p(n')\bigr)=\begin{cases} e-1,& p\nmid e,\\ e,& p\mid e. \end{cases}$$
ومن ثم، لقوة أولية منفردة،
$$g(p^e)=\begin{cases} p^{e-1},& p\nmid e,\\ p^e,& p\mid e. \end{cases}$$
وبما أن مساهمة كل أولي تُحدَّد مستقلّة عن غيرها، فإن \(g\) دالة مضاعِفية:
$$g(n)=\prod_{p^e\parallel n} g(p^e).$$
إذا ظهر أولي بأس \(1\) فقط، فإن \(g(p)=1\). لهذا لا تؤثر العوامل ذات الأس \(1\) في قيمة \(g(n)\). لذلك يُكتب كل عدد بشكل وحيد على الصورة
$$n=u\,s,$$
حيث يحتوي \(u\) على كل القوى الأولية ذات الأس \(e\ge 2\)، بينما \(s\) هو حاصل ضرب الأوليات الباقية ذات الأس \(1\).
وعندئذ
$$u\ \mathrm{powerful},\qquad s\ \mathrm{squarefree},\qquad \gcd(s,u)=1.$$
والأهم من ذلك أن
$$g(n)=g(u).$$
إذا عرّفنا
$$\operatorname{rad}(u)=\prod_{p\mid u} p,$$
فإن الذيول المسموح بها عند تثبيت \(u\) هي الأعداد \(s\le N/u\) الخالية من المربعات والمتباينة أوليًا مع \(\operatorname{rad}(u)\). لذلك نحصل على
$$S(N)+1=\sum_{\substack{u\le N\\ u\ \mathrm{powerful}}} g(u)\,Q\!\left(\left\lfloor\frac{N}{u}\right\rfloor,\operatorname{rad}(u)\right),$$
حيث
$$Q(x,r)=\#\{\,s\le x:\ s\ \mathrm{squarefree},\ \gcd(s,r)=1\,\}.$$
الواحد الزائد هنا يأتي من الحد \(n=1\)، ويُطرح في نهاية الحساب.
دالة المؤشر لكون العدد خاليًا من المربعات هي
$$1_{\mathrm{sqf}}(s)=\sum_{d^2\mid s}\mu(d).$$
ومن ثم تعطي معكوسية Möbius الصيغة
$$Q(x,r)=\sum_{\substack{d\le \sqrt{x}\\ \gcd(d,r)=1}}\mu(d)\,\Phi_r\!\left(\left\lfloor\frac{x}{d^2}\right\rfloor\right),$$
حيث تمثل \(\Phi_r(y)\) عدد الأعداد حتى \(y\) التي تكون متباينة أوليًا مع \(r\):
$$\Phi_r(y)=\#\{\,m\le y:\gcd(m,r)=1\,\}.$$
وبما أن \(r=\operatorname{rad}(u)\) عدد خالٍ من المربعات، فيمكن حساب \(\Phi_r(y)\) بالاحتواء والاستبعاد على قواسم \(r\):
$$\Phi_r(y)=\sum_{t\mid r}\mu(t)\left\lfloor\frac{y}{t}\right\rfloor.$$
وهذه هي البنية نفسها الموجودة في التنفيذ الكامل: مجموع Möbius خارجي لفرض شرط الخلو من المربعات، ومجموع داخلي على القواسم لفرض شرط التباين الأولي.
كل عدد powerful يمكن بناؤه باختيار أوليات مختلفة وإعطاء كل واحد منها أسًا لا يقل عن \(2\). لذا فإن بحثًا عمقيًا عبر الأوليات بترتيب متزايد يزور كل \(u\le N\) مرة واحدة فقط.
في كل عقدة من هذا البحث، يعرف التنفيذ القيمة الحالية لـ \(u\)، وقيمة \(g(u)\) المتراكمة، ومجموعة الأوليات التي تكوّن \(\operatorname{rad}(u)\). ومن ثم يضيف مباشرة
$$g(u)\,Q\!\left(\left\lfloor\frac{N}{u}\right\rfloor,\operatorname{rad}(u)\right)$$
إلى الجواب، ثم يحاول ضم أولي جديد مع أس \(2,3,4,\dots\) ما دام الناتج لا يتجاوز \(N\).
النسختان الأقصر تكتبان الفكرة نفسها على شكل عودية تلسكوبية مكافئة. فإذا وضعنا \(G_e(p)=g(p^e)\)، فإن كل زيادة
$$\Delta_e(p)=G_e(p)-G_{e-1}(p)$$
تُضرب في عدد المضاعفات التي ما زالت ممكنة، ثم تستمر العودية على الأوليات الأكبر فقط. هذه الصياغة المختصرة والصياغة الصريحة عبر \(Q(x,r)\) تعطيان المجموع نفسه تمامًا.
الأجزاء powerful التي لا تتجاوز \(20\) هي
$$1,\ 4,\ 8,\ 9,\ 16.$$
ثم نحسب عدد الذيول المسموح بها لكل جزء:
$$\begin{aligned} u=1&:&&g(u)=1,\quad Q(20,1)=13,\quad c=13,\\ u=4&:&&g(u)=4,\quad Q(5,2)=3,\quad c=12,\\ u=8&:&&g(u)=4,\quad Q(2,2)=1,\quad c=4,\\ u=9&:&&g(u)=3,\quad Q(2,3)=2,\quad c=6,\\ u=16&:&&g(u)=16,\quad Q(1,2)=1,\quad c=16. \end{aligned}$$
إذن
$$\sum_{n=1}^{20} g(n)=13+12+4+6+16=51,$$
وبعد حذف حد \(n=1\) نحصل على
$$\sum_{n=2}^{20}\gcd(n,n')=50.$$
تعتمد تطبيقات C++ وPython وJava كلها على الصيغة المحلية نفسها لـ \(g(p^e)\). يبدأ تنفيذ C++ بالتحقق من صحة هذه الصيغة على مجالات صغيرة، ثم يبني غربالًا للأوليات وقيم Möbius حتى \(\sqrt{N}\). بعد ذلك يُعدِّد الأجزاء powerful ببحث عمقي، ولكل جزء يحسب عدد الذيول الخالية من المربعات والمتباينة أوليًا مع \(\operatorname{rad}(u)\) باستخدام صيغة Möbius مع الاحتواء والاستبعاد المذكورة أعلاه. ويحتوي التنفيذ نفسه أيضًا على نسخة عودية أقصر لكنها مكافئة رياضيًا؛ وهذه هي الفكرة التي تعتمدها نسختا Python وJava. هناك يُزاد أس أولي واحد تدريجيًا، وكلما زاد العامل المحلي داخل \(\gcd\) أُضيف مقدار هذه الزيادة مضروبًا في عدد المضاعفات الباقية الممكنة وفي المساهمة العودية للأوليات الأكبر. اختلاف اللغات هنا هو اختلاف في شكل التنفيذ فقط، لا في المضمون الرياضي.
إذا وضعنا \(M=\lfloor\sqrt{N}\rfloor\)، فإن مرحلة الغربلة تحتاج إلى \(O(M\log\log M)\) زمنًا و\(O(M)\) ذاكرة. أما فضاء البحث الحقيقي فهو فضاء الأعداد powerful، وهو نادر جدًا مقارنةً بجميع الأعداد حتى \(N\)، لذلك لا تزور العودية إلا جزءًا صغيرًا جدًا من المجال الكامل. ولكل \(u\) مُزار، فإن عدّ الذيول يجري على قيم مربعة-حرة \(d\le \sqrt{N/u}\)، بينما يعمل الاحتواء والاستبعاد الداخلي فقط على قواسم \(\operatorname{rad}(u)\)، أي على مجموعة صغيرة من الأوليات. لهذا يكون الأداء أسرع بمراحل من المرور على كل \(n\le N\)، وهو ما يجعل \(N=5\times 10^{15}\) ممكنًا عمليًا.