For positive integers \(a=\prod p^{\alpha_p}\) and \(b=\prod p^{\beta_p}\), define their exponent difference by
$$D(a,b)=\sum_{p}\left|\alpha_p-\beta_p\right|=\sum_{p}\left|\nu_p(a)-\nu_p(b)\right|.$$
The task is to evaluate
$$S(n)=\sum_{1\le a,b\le n} D(a,b)$$
for the Project Euler instance \(n=10^{12}\), with the final result taken modulo \(10^9+7\). A direct scan of all ordered pairs would require \(10^{24}\) pair checks before even factoring the numbers, so the only viable route is to reorganize the sum prime by prime and prime-power by prime-power.
The key simplification is that the contribution of each prime is independent, and the absolute difference of two exponents can be rewritten as a count of divisibility levels.
Because only finitely many primes divide either \(a\) or \(b\), we may interchange the order of summation:
$$S(n)=\sum_{p\le n}\sum_{1\le a,b\le n}\left|\nu_p(a)-\nu_p(b)\right|.$$
For a fixed prime \(p\), define
$$T_p(n)=\sum_{1\le a,b\le n}\left|\nu_p(a)-\nu_p(b)\right|.$$
Then
$$S(n)=\sum_{p\le n} T_p(n).$$
So the entire problem reduces to understanding the contribution of one prime at a time.
For any nonnegative integers \(x\) and \(y\),
$$|x-y|=\sum_{k\ge 1}\left|\mathbf{1}_{x\ge k}-\mathbf{1}_{y\ge k}\right|.$$
Applied to \(x=\nu_p(a)\) and \(y=\nu_p(b)\), this says that level \(k\) contributes \(1\) exactly when one of \(a,b\) is divisible by \(p^k\) and the other is not. Therefore
$$T_p(n)=\sum_{k\ge 1} N_{p,k},$$
where \(N_{p,k}\) counts ordered pairs with different divisibility status by \(p^k\).
Let
$$q_{p,k}=\#\{m\le n: p^k\mid m\}=\left\lfloor\frac{n}{p^k}\right\rfloor.$$
There are \(q_{p,k}\) integers up to \(n\) divisible by \(p^k\), and \(n-q_{p,k}\) that are not. Since the pairs are ordered, the number of mismatches is
$$N_{p,k}=q_{p,k}(n-q_{p,k})+(n-q_{p,k})q_{p,k}=2q_{p,k}(n-q_{p,k}).$$
Hence
$$T_p(n)=\sum_{k\ge 1} 2q_{p,k}(n-q_{p,k}),\qquad q_{p,k}=\left\lfloor\frac{n}{p^k}\right\rfloor.$$
Only finitely many terms are nonzero, because \(q_{p,k}=0\) once \(p^k\gt n\).
Substituting the previous identity yields
$$\boxed{S(n)=2\sum_{p\le n}\sum_{k\ge 1}\left\lfloor\frac{n}{p^k}\right\rfloor\left(n-\left\lfloor\frac{n}{p^k}\right\rfloor\right).}$$
So every prime power \(p^k\le n\) contributes one term that depends only on the quotient \(\left\lfloor n/p^k\right\rfloor\). The algorithm is therefore an efficient way to enumerate those prime powers, or to group them when many primes share the same quotient.
The implementations choose the cutoff
$$B=\min(n,10^8).$$
For each prime \(p\le B\), all powers \(p,p^2,p^3,\dots\le n\) are visited explicitly. For the Project Euler instance \(n=10^{12}\), any prime \(p\gt B\) automatically satisfies \(p^2\gt n\), so such primes contribute only through the first power \(p\) itself.
That observation removes all higher-power work from the large-prime side.
For a large prime \(p\gt B\), the only relevant quantity is
$$q=\left\lfloor\frac{n}{p}\right\rfloor.$$
All primes producing the same quotient \(q\) lie in the interval
$$\frac{n}{q+1}\lt p\le \frac{n}{q}.$$
Instead of iterating through those primes individually, we count how many primes lie in that interval and multiply by their common contribution
$$2q(n-q).$$
The number of primes in any interval \([L,R]\) is
$$\pi(R)-\pi(L-1),$$
so fast evaluation of the prime-counting function \(\pi(x)\) completes the large-prime part.
The primes up to \(10\) are \(2,3,5,7\).
For \(p=2\), the nonzero quotients are \(q_{2,1}=5\), \(q_{2,2}=2\), \(q_{2,3}=1\), giving
$$2\cdot 5\cdot 5+2\cdot 2\cdot 8+2\cdot 1\cdot 9=50+32+18=100.$$
For \(p=3\), we get \(q_{3,1}=3\) and \(q_{3,2}=1\), so
$$2\cdot 3\cdot 7+2\cdot 1\cdot 9=42+18=60.$$
For \(p=5\), only \(q_{5,1}=2\) survives, contributing
$$2\cdot 2\cdot 8=32.$$
For \(p=7\), only \(q_{7,1}=1\) survives, contributing
$$2\cdot 1\cdot 9=18.$$
Therefore
$$S(10)=100+60+32+18=210,$$
which matches the small checkpoint used by the implementations.
The C++, Python, and Java implementations follow the same pipeline. First, they keep all arithmetic modulo \(10^9+7\) and reduce each contribution immediately using
$$2q(n-q)=2(nq-q^2).$$
Second, they generate all primes up to the cutoff \(B\) with a segmented sieve, which avoids storing a full primality table up to \(10^8\). For every small prime, the implementation repeatedly multiplies by the same prime to enumerate \(p^1,p^2,p^3,\dots\) until the next power would exceed \(n\), and each power contributes one quotient term.
Third, the remaining primes \(p\gt B\) are processed in quotient blocks. The implementation loops over possible values of \(\left\lfloor n/p\right\rfloor\), derives the corresponding interval of primes, counts how many primes lie there, and adds that many copies of the same contribution.
Finally, those interval counts come from a Lehmer-style prime-counting routine backed by a small precomputed sieve and memoization. That makes \(\pi(x)\) queries fast enough that the large-prime phase is tiny compared with any direct enumeration of pairs.
Let \(B=\min(n,10^8)\). The segmented sieve over small primes costs \(O(B\log\log B)\) time and uses memory proportional to the segment size, plus the small base sieve needed to mark composites. Enumerating the powers of each small prime adds one step for each prime power \(p^k\le n\) with \(p\le B\).
The large-prime phase iterates quotients up to \(n/(B+1)\). For the actual Project Euler input \(n=10^{12}\), that is only about \(10^4\) outer iterations. Each iteration performs prime-count queries at interval endpoints, and the memoized Lehmer method keeps those queries practical.
Overall, the method is many orders of magnitude faster than the naive \(O(n^2)\) scan over ordered pairs and is easily feasible for the required value \(10^{12}\).
Für positive ganze Zahlen \(a=\prod p^{\alpha_p}\) und \(b=\prod p^{\beta_p}\) sei die Exponentendifferenz definiert als
$$D(a,b)=\sum_{p}\left|\alpha_p-\beta_p\right|=\sum_{p}\left|\nu_p(a)-\nu_p(b)\right|.$$
Gesucht ist
$$S(n)=\sum_{1\le a,b\le n} D(a,b)$$
für den Project-Euler-Fall \(n=10^{12}\), wobei das Endergebnis modulo \(10^9+7\) genommen wird. Ein direkter Durchlauf aller geordneten Paare wäre quadratisch und völlig unpraktikabel, daher ordnet die Lösung die Summe nach Primzahlen und Primzahlpotenzen neu.
Die entscheidende Vereinfachung ist, dass sich jeder Primbeitrag unabhängig behandeln lässt und sich die Betragsdifferenz zweier Exponenten als Zählung von Teilbarkeitsstufen schreiben lässt.
Da nur endlich viele Primzahlen \(a\) oder \(b\) teilen, dürfen wir die Summen vertauschen:
$$S(n)=\sum_{p\le n}\sum_{1\le a,b\le n}\left|\nu_p(a)-\nu_p(b)\right|.$$
Für eine feste Primzahl \(p\) definieren wir
$$T_p(n)=\sum_{1\le a,b\le n}\left|\nu_p(a)-\nu_p(b)\right|.$$
Dann gilt
$$S(n)=\sum_{p\le n} T_p(n).$$
Damit reduziert sich das Gesamtproblem auf die Analyse einer einzelnen Primzahl.
Für nichtnegative ganze Zahlen \(x\) und \(y\) gilt
$$|x-y|=\sum_{k\ge 1}\left|\mathbf{1}_{x\ge k}-\mathbf{1}_{y\ge k}\right|.$$
Setzt man \(x=\nu_p(a)\) und \(y=\nu_p(b)\), dann liefert eine Stufe \(k\) genau dann den Beitrag \(1\), wenn genau eine der beiden Zahlen durch \(p^k\) teilbar ist. Also
$$T_p(n)=\sum_{k\ge 1} N_{p,k},$$
wobei \(N_{p,k}\) die Anzahl geordneter Paare mit unterschiedlichem Teilbarkeitsstatus bezüglich \(p^k\) bezeichnet.
Setze
$$q_{p,k}=\#\{m\le n: p^k\mid m\}=\left\lfloor\frac{n}{p^k}\right\rfloor.$$
Es gibt also \(q_{p,k}\) Zahlen bis \(n\), die durch \(p^k\) teilbar sind, und \(n-q_{p,k}\) Zahlen, die es nicht sind. Da die Paare geordnet sind, ist die Zahl der Abweichungen
$$N_{p,k}=q_{p,k}(n-q_{p,k})+(n-q_{p,k})q_{p,k}=2q_{p,k}(n-q_{p,k}).$$
Folglich
$$T_p(n)=\sum_{k\ge 1} 2q_{p,k}(n-q_{p,k}),\qquad q_{p,k}=\left\lfloor\frac{n}{p^k}\right\rfloor.$$
Sobald \(p^k\gt n\) ist, wird \(q_{p,k}=0\), daher bleiben nur endlich viele Summanden.
Einsetzen liefert
$$\boxed{S(n)=2\sum_{p\le n}\sum_{k\ge 1}\left\lfloor\frac{n}{p^k}\right\rfloor\left(n-\left\lfloor\frac{n}{p^k}\right\rfloor\right).}$$
Jede Primzahlpotenz \(p^k\le n\) liefert also genau einen Term, der nur vom Quotienten \(\left\lfloor n/p^k\right\rfloor\) abhängt. Der Algorithmus muss daher nur diese Primzahlpotenzen effizient erzeugen oder gleichartige Quotienten gemeinsam behandeln.
Die Implementierungen wählen die Schranke
$$B=\min(n,10^8).$$
Für jede Primzahl \(p\le B\) werden alle Potenzen \(p,p^2,p^3,\dots\le n\) direkt durchlaufen. Im Project-Euler-Fall \(n=10^{12}\) gilt für jede Primzahl \(p\gt B\) automatisch \(p^2\gt n\), also kann eine solche Primzahl nur noch über die erste Potenz \(p\) selbst beitragen.
Dadurch verschwindet auf der Seite der großen Primzahlen jede höhere Potenz vollständig aus der Rechnung.
Für eine große Primzahl \(p\gt B\) ist nur noch
$$q=\left\lfloor\frac{n}{p}\right\rfloor$$
relevant. Alle Primzahlen mit demselben Quotienten \(q\) liegen im Intervall
$$\frac{n}{q+1}\lt p\le \frac{n}{q}.$$
Statt diese Primzahlen einzeln aufzulisten, zählt die Implementierung, wie viele Primzahlen in diesem Intervall liegen, und multipliziert mit dem gemeinsamen Beitrag
$$2q(n-q).$$
Die Anzahl der Primzahlen in \([L,R]\) ist
$$\pi(R)-\pi(L-1),$$
weshalb eine schnelle Auswertung der Primzahlzählfunktion \(\pi(x)\) die große-Primzahlen-Phase abschließt.
Die Primzahlen bis \(10\) sind \(2,3,5,7\).
Für \(p=2\) sind die nichtverschwindenden Quotienten \(q_{2,1}=5\), \(q_{2,2}=2\), \(q_{2,3}=1\), also
$$2\cdot 5\cdot 5+2\cdot 2\cdot 8+2\cdot 1\cdot 9=50+32+18=100.$$
Für \(p=3\) erhält man \(q_{3,1}=3\) und \(q_{3,2}=1\), also
$$2\cdot 3\cdot 7+2\cdot 1\cdot 9=42+18=60.$$
Für \(p=5\) bleibt nur \(q_{5,1}=2\), mit Beitrag
$$2\cdot 2\cdot 8=32.$$
Für \(p=7\) bleibt nur \(q_{7,1}=1\), mit Beitrag
$$2\cdot 1\cdot 9=18.$$
Damit
$$S(10)=100+60+32+18=210,$$
genau der kleine Kontrollwert aus den Implementierungen.
Die C++-, Python- und Java-Implementierungen folgen derselben Pipeline. Zuerst wird jede Rechnung modulo \(10^9+7\) geführt, wobei jeder Summand sofort über
$$2q(n-q)=2(nq-q^2)$$
reduziert wird. Danach erzeugt ein segmentiertes Sieb alle Primzahlen bis zur Schranke \(B\), ohne eine vollständige Primtabelle bis \(10^8\) im Speicher zu halten. Für jede kleine Primzahl wird wiederholt mit derselben Primzahl multipliziert, um \(p^1,p^2,p^3,\dots\) zu durchlaufen, bis die nächste Potenz größer als \(n\) wäre.
Anschließend verarbeitet die Implementierung die restlichen Primzahlen \(p\gt B\) in Quotientenblöcken. Sie iteriert über mögliche Werte von \(\left\lfloor n/p\right\rfloor\), bestimmt daraus das zugehörige Primintervall, zählt die Primzahlen in diesem Bereich und addiert so viele identische Beiträge.
Diese Intervallzahlen stammen aus einem Primzahlzählverfahren im Stil von Lehmer, gestützt durch ein kleines vorberechnetes Sieb und Memoisierung. Genau dadurch wird der große-Primzahlen-Teil schnell genug, obwohl \(n\) selbst riesig ist.
Mit \(B=\min(n,10^8)\) kostet das segmentierte Sieb \(O(B\log\log B)\) Zeit und Speicher in der Größenordnung der Segmentgröße, zusätzlich zu einem kleinen Basissieb. Das Durchlaufen der Potenzen kleiner Primzahlen fügt je Primzahlpotenz \(p^k\le n\) mit \(p\le B\) nur einen weiteren Schritt hinzu.
Die Phase für große Primzahlen iteriert Quotienten bis \(n/(B+1)\). Für den tatsächlichen Project-Euler-Wert \(n=10^{12}\) sind das nur ungefähr \(10^4\) äußere Schleifendurchläufe. Jede Iteration braucht Primzahlzählungen an Intervallgrenzen, und die memoisierten Lehmer-Auswertungen halten diese Anfragen praktisch schnell.
Insgesamt ist das Verfahren um viele Größenordnungen schneller als ein naiver \(O(n^2)\)-Durchlauf über alle geordneten Paare und für \(10^{12}\) gut machbar.
Pozitif tamsayılar \(a=\prod p^{\alpha_p}\) ve \(b=\prod p^{\beta_p}\) için üs farkını
$$D(a,b)=\sum_{p}\left|\alpha_p-\beta_p\right|=\sum_{p}\left|\nu_p(a)-\nu_p(b)\right|$$
şeklinde tanımlayalım. İstenen büyüklük
$$S(n)=\sum_{1\le a,b\le n} D(a,b)$$
olup Project Euler girdisinde \(n=10^{12}\) için hesaplanır ve sonuç \(10^9+7\) modunda verilir. Tüm sıralı çiftleri doğrudan gezmek hem çift sayısı hem de çarpanlara ayırma maliyeti yüzünden imkansızdır; bu yüzden çözüm toplamı asal bazında ve asal kuvvet bazında yeniden düzenler.
Temel fikir şudur: her asalın katkısı bağımsızdır ve iki üs arasındaki mutlak fark, bölünebilme seviyelerini sayarak yeniden yazılabilir.
Yalnızca sonlu sayıda asal \(a\) veya \(b\)'yi böldüğü için toplamların sırasını değiştirebiliriz:
$$S(n)=\sum_{p\le n}\sum_{1\le a,b\le n}\left|\nu_p(a)-\nu_p(b)\right|.$$
Sabit bir asal \(p\) için
$$T_p(n)=\sum_{1\le a,b\le n}\left|\nu_p(a)-\nu_p(b)\right|$$
tanımını yapalım. O zaman
$$S(n)=\sum_{p\le n} T_p(n).$$
Böylece tüm problem, tek bir asala ait katkıyı anlamaya indirgenir.
Negatif olmayan tamsayılar \(x\) ve \(y\) için
$$|x-y|=\sum_{k\ge 1}\left|\mathbf{1}_{x\ge k}-\mathbf{1}_{y\ge k}\right|$$
kimliği geçerlidir. Bunu \(x=\nu_p(a)\) ve \(y=\nu_p(b)\) için uygularsak, \(k\). seviye tam olarak şu durumda \(1\) katkı verir: sayılardan biri \(p^k\)'ya bölünürken diğeri bölünmez. Dolayısıyla
$$T_p(n)=\sum_{k\ge 1} N_{p,k},$$
burada \(N_{p,k}\), \(p^k\)'ya göre bölünebilme durumu farklı olan sıralı çiftlerin sayısıdır.
Şunu tanımlayalım:
$$q_{p,k}=\#\{m\le n: p^k\mid m\}=\left\lfloor\frac{n}{p^k}\right\rfloor.$$
Yani \(1\) ile \(n\) arasında \(p^k\)'ya bölünen \(q_{p,k}\) sayı, bölünmeyen ise \(n-q_{p,k}\) sayı vardır. Çiftler sıralı olduğu için uyumsuzluk sayısı
$$N_{p,k}=q_{p,k}(n-q_{p,k})+(n-q_{p,k})q_{p,k}=2q_{p,k}(n-q_{p,k})$$
olur. Demek ki
$$T_p(n)=\sum_{k\ge 1} 2q_{p,k}(n-q_{p,k}),\qquad q_{p,k}=\left\lfloor\frac{n}{p^k}\right\rfloor.$$
\(p^k\gt n\) olduğunda \(q_{p,k}=0\) olur; bu yüzden toplam sonlu sayıda terim içerir.
Önceki eşitliği yerine koyunca
$$\boxed{S(n)=2\sum_{p\le n}\sum_{k\ge 1}\left\lfloor\frac{n}{p^k}\right\rfloor\left(n-\left\lfloor\frac{n}{p^k}\right\rfloor\right)}$$
elde edilir. Böylece her \(p^k\le n\) asal kuvveti, yalnızca \(\left\lfloor n/p^k\right\rfloor\) bölümüne bağlı tek bir katkı üretir. Algoritmanın yapması gereken şey, bu asal kuvvetleri verimli biçimde gezmek veya aynı bölümü verenleri topluca saymaktır.
Uygulamalar
$$B=\min(n,10^8)$$
eşik değerini kullanır. \(p\le B\) olan her asal için \(p,p^2,p^3,\dots\le n\) kuvvetleri doğrudan dolaşılır. Project Euler örneğinde \(n=10^{12}\) olduğu için \(p\gt B\) olan her asal otomatik olarak \(p^2\gt n\) şartını sağlar; yani bu asallar yalnızca birinci kuvvetleriyle katkı verebilir.
Böylece büyük asallar tarafında yüksek kuvvetleri düşünmeye gerek kalmaz.
\(p\gt B\) olan büyük bir asal için tek önemli nicelik
$$q=\left\lfloor\frac{n}{p}\right\rfloor$$
değeridir. Aynı \(q\) değerini veren bütün asallar
$$\frac{n}{q+1}\lt p\le \frac{n}{q}$$
aralığında yer alır. Bu asalları tek tek dolaşmak yerine, bu aralıkta kaç asal olduğunu sayıp ortak katkı
$$2q(n-q)$$
ile çarparız. Herhangi bir \([L,R]\) aralığındaki asal sayısı
$$\pi(R)-\pi(L-1)$$
olduğu için, hızlı bir \(\pi(x)\) hesabı büyük-asal kısmını bitirmeye yeter.
\(10\)'a kadar olan asallar \(2,3,5,7\)'dir.
\(p=2\) için sıfır olmayan bölümler \(q_{2,1}=5\), \(q_{2,2}=2\), \(q_{2,3}=1\) olur ve katkı
$$2\cdot 5\cdot 5+2\cdot 2\cdot 8+2\cdot 1\cdot 9=50+32+18=100$$
çıkar. \(p=3\) için \(q_{3,1}=3\) ve \(q_{3,2}=1\) olduğundan
$$2\cdot 3\cdot 7+2\cdot 1\cdot 9=42+18=60.$$
\(p=5\) için yalnızca \(q_{5,1}=2\) kalır ve katkı
$$2\cdot 2\cdot 8=32$$
olur. \(p=7\) için yalnızca \(q_{7,1}=1\) kalır ve katkı
$$2\cdot 1\cdot 9=18$$
verir. Sonuç olarak
$$S(10)=100+60+32+18=210,$$
ki bu, uygulamaların kullandığı küçük doğrulama değeriyle aynıdır.
C++, Python ve Java uygulamaları aynı akışı izler. Önce tüm aritmetik \(10^9+7\) modunda yürütülür ve her katkı
$$2q(n-q)=2(nq-q^2)$$
biçiminde hemen indirgenir. Sonra segmentli asal eleği ile \(B\)'ye kadar tüm asallar üretilir; böylece \(10^8\)'e kadar tam bir asal dizisi bellekte tutulmaz. Her küçük asal için aynı asal ile tekrar tekrar çarpılarak \(p^1,p^2,p^3,\dots\) kuvvetleri gezilir ve her kuvvet bir bölüm terimi üretir.
Daha sonra \(p\gt B\) olan asallar bölüm blokları halinde işlenir. Uygulama olası \(\left\lfloor n/p\right\rfloor\) değerlerini dolaşır, buna karşılık gelen asal aralığını çıkarır, bu aralıkta kaç asal olduğunu sayar ve aynı katkıdan o kadar ekler.
Bu aralık sayımları Lehmer tarzı bir asal sayma yönteminden gelir; yöntem küçük bir ön elekle ve önbellekleme ile desteklenir. Bu sayede çok büyük \(n\) için bile büyük-asal kısmı pratik kalır.
\(B=\min(n,10^8)\) olsun. Segmentli eleme \(O(B\log\log B)\) zamanda çalışır ve bellek kullanımı segment boyutu ile küçük bir taban eleğinin toplamı mertebesindedir. Küçük asalların kuvvetlerini gezmek de \(p\le B\) ve \(p^k\le n\) olan her asal kuvvet için yalnızca bir ek adım gerektirir.
Büyük-asal aşaması, bölümleri en fazla \(n/(B+1)\) değerine kadar dolaşır. Gerçek Project Euler girdisi \(n=10^{12}\) için bu yaklaşık \(10^4\) dış döngü anlamına gelir. Her döngüde aralık uçlarında asal sayımı yapılır ve önbellekli Lehmer yöntemi bunu pratik hızda tutar.
Sonuç olarak yöntem, sıralı çiftleri kaba kuvvetle \(O(n^2)\) gezmeye göre astronomik derecede daha hızlıdır ve gerekli \(10^{12}\) değeri için rahatlıkla uygulanabilir.
Para enteros positivos \(a=\prod p^{\alpha_p}\) y \(b=\prod p^{\beta_p}\), definimos la diferencia de exponentes como
$$D(a,b)=\sum_{p}\left|\alpha_p-\beta_p\right|=\sum_{p}\left|\nu_p(a)-\nu_p(b)\right|.$$
La cantidad buscada es
$$S(n)=\sum_{1\le a,b\le n} D(a,b),$$
evaluada para el caso de Project Euler \(n=10^{12}\), tomando el resultado final módulo \(10^9+7\). Recorrer directamente todos los pares ordenados sería cuadrático y totalmente inviable, de modo que la solución reorganiza la suma por primo y por potencia prima.
La observación decisiva es que la contribución de cada primo es independiente, y que la diferencia absoluta de exponentes puede reinterpretarse como un conteo de niveles de divisibilidad.
Como solo un número finito de primos divide a \(a\) o a \(b\), podemos intercambiar el orden de las sumas:
$$S(n)=\sum_{p\le n}\sum_{1\le a,b\le n}\left|\nu_p(a)-\nu_p(b)\right|.$$
Para un primo fijo \(p\), definimos
$$T_p(n)=\sum_{1\le a,b\le n}\left|\nu_p(a)-\nu_p(b)\right|.$$
Entonces
$$S(n)=\sum_{p\le n} T_p(n).$$
Así, todo el problema se reduce a entender la contribución de un solo primo.
Para enteros no negativos \(x\) e \(y\), se cumple
$$|x-y|=\sum_{k\ge 1}\left|\mathbf{1}_{x\ge k}-\mathbf{1}_{y\ge k}\right|.$$
Aplicado a \(x=\nu_p(a)\) e \(y=\nu_p(b)\), esto significa que el nivel \(k\) aporta \(1\) exactamente cuando uno de \(a,b\) es divisible por \(p^k\) y el otro no. Por lo tanto,
$$T_p(n)=\sum_{k\ge 1} N_{p,k},$$
donde \(N_{p,k}\) es el número de pares ordenados con distinta divisibilidad por \(p^k\).
Sea
$$q_{p,k}=\#\{m\le n: p^k\mid m\}=\left\lfloor\frac{n}{p^k}\right\rfloor.$$
Hay \(q_{p,k}\) enteros hasta \(n\) divisibles por \(p^k\), y \(n-q_{p,k}\) que no lo son. Como los pares son ordenados, el número de discrepancias es
$$N_{p,k}=q_{p,k}(n-q_{p,k})+(n-q_{p,k})q_{p,k}=2q_{p,k}(n-q_{p,k}).$$
En consecuencia,
$$T_p(n)=\sum_{k\ge 1} 2q_{p,k}(n-q_{p,k}),\qquad q_{p,k}=\left\lfloor\frac{n}{p^k}\right\rfloor.$$
Solo hay un número finito de términos no nulos, porque cuando \(p^k\gt n\), el cociente se vuelve \(0\).
Sustituyendo lo anterior obtenemos
$$\boxed{S(n)=2\sum_{p\le n}\sum_{k\ge 1}\left\lfloor\frac{n}{p^k}\right\rfloor\left(n-\left\lfloor\frac{n}{p^k}\right\rfloor\right).}$$
Así, cada potencia prima \(p^k\le n\) aporta un término que depende solo del cociente \(\left\lfloor n/p^k\right\rfloor\). El algoritmo se limita entonces a enumerar esas potencias de manera eficiente o a agrupar muchas de ellas cuando comparten el mismo cociente.
Las implementaciones eligen el umbral
$$B=\min(n,10^8).$$
Para cada primo \(p\le B\), se recorren explícitamente todas las potencias \(p,p^2,p^3,\dots\le n\). En el caso concreto \(n=10^{12}\), cualquier primo \(p\gt B\) satisface automáticamente \(p^2\gt n\), de modo que esos primos solo pueden contribuir mediante la primera potencia \(p\).
Esa observación elimina todo trabajo de potencias altas en la parte de primos grandes.
Para un primo grande \(p\gt B\), la única cantidad relevante es
$$q=\left\lfloor\frac{n}{p}\right\rfloor.$$
Todos los primos que producen el mismo valor \(q\) están en el intervalo
$$\frac{n}{q+1}\lt p\le \frac{n}{q}.$$
En lugar de iterar uno por uno, contamos cuántos primos hay en ese intervalo y multiplicamos por la contribución común
$$2q(n-q).$$
La cantidad de primos en \([L,R]\) es
$$\pi(R)-\pi(L-1),$$
así que una evaluación rápida de la función contadora de primos \(\pi(x)\) basta para cerrar la parte de primos grandes.
Los primos hasta \(10\) son \(2,3,5,7\).
Para \(p=2\), los cocientes no nulos son \(q_{2,1}=5\), \(q_{2,2}=2\), \(q_{2,3}=1\), y por tanto
$$2\cdot 5\cdot 5+2\cdot 2\cdot 8+2\cdot 1\cdot 9=50+32+18=100.$$
Para \(p=3\), obtenemos \(q_{3,1}=3\) y \(q_{3,2}=1\), así que
$$2\cdot 3\cdot 7+2\cdot 1\cdot 9=42+18=60.$$
Para \(p=5\), solo queda \(q_{5,1}=2\), con contribución
$$2\cdot 2\cdot 8=32.$$
Para \(p=7\), solo queda \(q_{7,1}=1\), con contribución
$$2\cdot 1\cdot 9=18.$$
Por lo tanto,
$$S(10)=100+60+32+18=210,$$
que coincide con la comprobación pequeña usada por las implementaciones.
Las implementaciones en C++, Python y Java siguen el mismo esquema. Primero, toda la aritmética se mantiene módulo \(10^9+7\), y cada contribución se reduce inmediatamente usando
$$2q(n-q)=2(nq-q^2).$$
Después, una criba segmentada genera todos los primos hasta el umbral \(B\), evitando guardar una tabla completa de primalidad hasta \(10^8\). Para cada primo pequeño, la implementación multiplica repetidamente por ese mismo primo para visitar \(p^1,p^2,p^3,\dots\) hasta que la siguiente potencia ya supera \(n\).
Luego, los primos restantes \(p\gt B\) se procesan por bloques de cociente. La implementación recorre los posibles valores de \(\left\lfloor n/p\right\rfloor\), deduce el intervalo correspondiente de primos, cuenta cuántos primos hay allí y suma esa cantidad de contribuciones idénticas.
Esos conteos de intervalos provienen de un método de conteo de primos de tipo Lehmer, apoyado por una criba pequeña precomputada y memoización. Gracias a ello, la fase de primos grandes sigue siendo rápida aunque \(n\) sea enorme.
Sea \(B=\min(n,10^8)\). La criba segmentada sobre los primos pequeños cuesta \(O(B\log\log B)\) tiempo y memoria proporcional al tamaño del segmento, más una criba base pequeña. Enumerar las potencias de cada primo pequeño añade un paso por cada potencia prima \(p^k\le n\) con \(p\le B\).
La fase de primos grandes itera cocientes hasta \(n/(B+1)\). Para el valor real de Project Euler \(n=10^{12}\), eso significa solo unas \(10^4\) iteraciones externas. Cada una necesita consultas de conteo de primos en los extremos del intervalo, y el método de Lehmer con memoización hace que esas consultas sean prácticas.
En conjunto, el método es muchísimos órdenes de magnitud más rápido que la enumeración ingenua \(O(n^2)\) de pares ordenados y resulta perfectamente viable para \(10^{12}\).
对正整数 \(a=\prod p^{\alpha_p}\) 和 \(b=\prod p^{\beta_p}\)，定义它们的指数差为
$$D(a,b)=\sum_{p}\left|\alpha_p-\beta_p\right|=\sum_{p}\left|\nu_p(a)-\nu_p(b)\right|.$$
题目要求计算
$$S(n)=\sum_{1\le a,b\le n} D(a,b),$$
并在 Project Euler 的目标输入 \(n=10^{12}\) 下给出模 \(10^9+7\) 的结果。若直接枚举所有有序对，不仅有 \(10^{24}\) 对需要处理，还要反复分解整数，完全不可行。所以正确做法是把总和拆成“按素数分组”和“按素数幂分组”的形式。
核心观察有两个：第一，不同素数的贡献彼此独立；第二，两个指数的绝对差可以改写成若干“是否被 \(p^k\) 整除”的层级计数。
由于能整除 \(a\) 或 \(b\) 的素数只有有限多个，我们可以交换求和顺序：
$$S(n)=\sum_{p\le n}\sum_{1\le a,b\le n}\left|\nu_p(a)-\nu_p(b)\right|.$$
对固定素数 \(p\)，记
$$T_p(n)=\sum_{1\le a,b\le n}\left|\nu_p(a)-\nu_p(b)\right|.$$
于是有
$$S(n)=\sum_{p\le n} T_p(n).$$
这样一来，整道题就变成了“如何高效求出单个素数 \(p\) 的总贡献”。
对任意非负整数 \(x\) 和 \(y\)，有恒等式
$$|x-y|=\sum_{k\ge 1}\left|\mathbf{1}_{x\ge k}-\mathbf{1}_{y\ge k}\right|.$$
把它代入 \(x=\nu_p(a)\)、\(y=\nu_p(b)\) 后可以看到：当且仅当 \(a,b\) 中恰有一个能被 \(p^k\) 整除时，第 \(k\) 层贡献 \(1\)。因此
$$T_p(n)=\sum_{k\ge 1} N_{p,k},$$
其中 \(N_{p,k}\) 表示在 \(1\le a,b\le n\) 的所有有序对中，对 \(p^k\) 的整除性不同的那些对的数量。
定义
$$q_{p,k}=\#\{m\le n: p^k\mid m\}=\left\lfloor\frac{n}{p^k}\right\rfloor.$$
也就是说，从 \(1\) 到 \(n\) 中，有 \(q_{p,k}\) 个数能被 \(p^k\) 整除，另外 \(n-q_{p,k}\) 个不能。因为这里统计的是有序对，所以整除状态不一致的对数是
$$N_{p,k}=q_{p,k}(n-q_{p,k})+(n-q_{p,k})q_{p,k}=2q_{p,k}(n-q_{p,k}).$$
于是
$$T_p(n)=\sum_{k\ge 1} 2q_{p,k}(n-q_{p,k}),\qquad q_{p,k}=\left\lfloor\frac{n}{p^k}\right\rfloor.$$
当 \(p^k\gt n\) 时，\(q_{p,k}=0\)，因此这实际上是一个有限和。
把上式代回总和，就得到
$$\boxed{S(n)=2\sum_{p\le n}\sum_{k\ge 1}\left\lfloor\frac{n}{p^k}\right\rfloor\left(n-\left\lfloor\frac{n}{p^k}\right\rfloor\right).}$$
现在问题的本质已经很清楚了：每一个满足 \(p^k\le n\) 的素数幂都贡献一个只依赖于 \(\left\lfloor n/p^k\right\rfloor\) 的项。算法要做的，不是枚举所有 \((a,b)\)，而是高效地枚举这些素数幂，或者在很多素数拥有相同商值时把它们合并处理。
实现中选取阈值
$$B=\min(n,10^8).$$
当 \(p\le B\) 时，直接枚举所有 \(p,p^2,p^3,\dots\le n\) 的素数幂即可。对题目的目标值 \(n=10^{12}\) 来说，只要 \(p\gt B\)，就自动满足 \(p^2\gt n\)。这意味着大素数不可能以平方或更高次幂出现，它们的贡献只来自一次幂 \(p\) 本身。
这个事实非常关键，因为它让“大素数部分”从枚举若干幂次，简化成只处理一次整除层级。
对大素数 \(p\gt B\)，唯一需要关心的是
$$q=\left\lfloor\frac{n}{p}\right\rfloor.$$
所有给出同一个 \(q\) 的素数，都落在区间
$$\frac{n}{q+1}\lt p\le \frac{n}{q}$$
里。于是没有必要逐个检查这些素数，只要统计该区间内有多少个素数，再乘上共同贡献
$$2q(n-q)$$
即可。任意区间 \([L,R]\) 中的素数个数都可写成
$$\pi(R)-\pi(L-1),$$
所以整个大素数部分归结为若干次快速的素数计数查询。
\(10\) 以内的素数是 \(2,3,5,7\)。
对 \(p=2\)，非零商值为 \(q_{2,1}=5\)、\(q_{2,2}=2\)、\(q_{2,3}=1\)，因此贡献为
$$2\cdot 5\cdot 5+2\cdot 2\cdot 8+2\cdot 1\cdot 9=50+32+18=100.$$
对 \(p=3\)，有 \(q_{3,1}=3\)、\(q_{3,2}=1\)，所以贡献为
$$2\cdot 3\cdot 7+2\cdot 1\cdot 9=42+18=60.$$
对 \(p=5\)，只有 \(q_{5,1}=2\)，贡献为
$$2\cdot 2\cdot 8=32.$$
对 \(p=7\)，只有 \(q_{7,1}=1\)，贡献为
$$2\cdot 1\cdot 9=18.$$
因此
$$S(10)=100+60+32+18=210,$$
这与实现中的小规模校验完全一致。
C++、Python 和 Java 实现采用的是同一条计算流水线。第一步，所有运算都在模 \(10^9+7\) 下进行，并把每个贡献立刻写成
$$2q(n-q)=2(nq-q^2)$$
的形式来累计。第二步，用分段筛生成阈值 \(B\) 以内的全部素数，这样就不必为 \(10^8\) 范围建立一整张完整布尔表。对每个小素数，再不断乘以自身，依次访问 \(p^1,p^2,p^3,\dots\)，直到下一次乘法会超过 \(n\)。
第三步，对剩余的大素数 \(p\gt B\)，实现并不逐个枚举，而是按 \(\left\lfloor n/p\right\rfloor\) 的相同取值进行分块。每一块先确定对应的素数区间，再查询该区间内的素数个数，然后一次性加入同样的贡献。
最后，这些区间素数个数通过 Lehmer 风格的素数计数算法获得。该算法依赖一个较小的预筛结果，并配合记忆化缓存来重复利用中间结果。正是这一层快速的 \(\pi(x)\) 计算，使得大素数部分在实际输入下非常轻量。
记 \(B=\min(n,10^8)\)。小素数部分的分段筛时间复杂度为 \(O(B\log\log B)\)，内存开销与分段大小成正比，再加上一份较小的基础筛。枚举小素数的幂次时，对每个满足 \(p\le B\)、\(p^k\le n\) 的素数幂只会访问一次。
大素数部分最多枚举到商值 \(n/(B+1)\)。对本题的实际输入 \(n=10^{12}\) 来说，这大约只是 \(10^4\) 次外层循环。每次循环需要做若干个区间端点的素数计数，而带缓存的 Lehmer 方法足以把这部分控制在可接受范围内。
因此，整体算法与朴素的 \(O(n^2)\) 有序对枚举相比快了许多个数量级，完全适用于题目要求的 \(10^{12}\) 规模。
Для положительных целых чисел \(a=\prod p^{\alpha_p}\) и \(b=\prod p^{\beta_p}\) определим разность показателей как
$$D(a,b)=\sum_{p}\left|\alpha_p-\beta_p\right|=\sum_{p}\left|\nu_p(a)-\nu_p(b)\right|.$$
Требуется вычислить
$$S(n)=\sum_{1\le a,b\le n} D(a,b)$$
для случая Project Euler \(n=10^{12}\), а финальный ответ взять по модулю \(10^9+7\). Полный перебор всех упорядоченных пар невозможен: самих пар слишком много, а ещё пришлось бы постоянно разлагать числа на простые множители. Поэтому решение перестраивает сумму по простым числам и по их степеням.
Главная идея состоит в том, что вклад каждого простого можно рассматривать отдельно, а абсолютную разность показателей можно переписать через уровни делимости.
Поскольку только конечное число простых делит \(a\) или \(b\), порядок суммирования можно поменять:
$$S(n)=\sum_{p\le n}\sum_{1\le a,b\le n}\left|\nu_p(a)-\nu_p(b)\right|.$$
Для фиксированного простого \(p\) обозначим
$$T_p(n)=\sum_{1\le a,b\le n}\left|\nu_p(a)-\nu_p(b)\right|.$$
Тогда
$$S(n)=\sum_{p\le n} T_p(n).$$
Итак, нужно научиться быстро считать вклад одной фиксированной простоты.
Для любых неотрицательных целых \(x\) и \(y\) верно
$$|x-y|=\sum_{k\ge 1}\left|\mathbf{1}_{x\ge k}-\mathbf{1}_{y\ge k}\right|.$$
Если подставить сюда \(x=\nu_p(a)\) и \(y=\nu_p(b)\), то уровень \(k\) даёт вклад \(1\) ровно тогда, когда одно из чисел делится на \(p^k\), а другое нет. Следовательно,
$$T_p(n)=\sum_{k\ge 1} N_{p,k},$$
где \(N_{p,k}\) означает число упорядоченных пар с разной делимостью на \(p^k\).
Положим
$$q_{p,k}=\#\{m\le n: p^k\mid m\}=\left\lfloor\frac{n}{p^k}\right\rfloor.$$
Тогда до \(n\) имеется \(q_{p,k}\) чисел, кратных \(p^k\), и \(n-q_{p,k}\) чисел, не кратных \(p^k\). Поскольку пары упорядочены, число несовпадений равно
$$N_{p,k}=q_{p,k}(n-q_{p,k})+(n-q_{p,k})q_{p,k}=2q_{p,k}(n-q_{p,k}).$$
Значит,
$$T_p(n)=\sum_{k\ge 1} 2q_{p,k}(n-q_{p,k}),\qquad q_{p,k}=\left\lfloor\frac{n}{p^k}\right\rfloor.$$
Когда \(p^k\gt n\), получаем \(q_{p,k}=0\), поэтому сумма конечна.
Подставляя это в общую формулу, получаем
$$\boxed{S(n)=2\sum_{p\le n}\sum_{k\ge 1}\left\lfloor\frac{n}{p^k}\right\rfloor\left(n-\left\lfloor\frac{n}{p^k}\right\rfloor\right).}$$
Теперь каждая степень простого \(p^k\le n\) даёт один вклад, зависящий только от частного \(\left\lfloor n/p^k\right\rfloor\). Значит, алгоритм должен лишь эффективно перечислить эти степени или объединить те из них, у которых частное совпадает.
В реализации используется порог
$$B=\min(n,10^8).$$
Для каждого простого \(p\le B\) явно перебираются все степени \(p,p^2,p^3,\dots\le n\). Для целевого значения \(n=10^{12}\) любой простой \(p\gt B\) автоматически удовлетворяет неравенству \(p^2\gt n\), поэтому такой простой может внести вклад только через первую степень \(p\).
Это полностью убирает высокие степени из части с большими простыми.
Для большого простого \(p\gt B\) важна только величина
$$q=\left\lfloor\frac{n}{p}\right\rfloor.$$
Все простые, дающие одно и то же \(q\), лежат в интервале
$$\frac{n}{q+1}\lt p\le \frac{n}{q}.$$
Поэтому вместо поштучного перебора можно просто посчитать число простых в этом интервале и умножить на общий вклад
$$2q(n-q).$$
Число простых на отрезке \([L,R]\) равно
$$\pi(R)-\pi(L-1),$$
и именно поэтому быстрая функция подсчёта простых \(\pi(x)\) закрывает всю часть с большими простыми.
Простые числа до \(10\) равны \(2,3,5,7\).
Для \(p=2\) ненулевые частные таковы: \(q_{2,1}=5\), \(q_{2,2}=2\), \(q_{2,3}=1\), поэтому вклад равен
$$2\cdot 5\cdot 5+2\cdot 2\cdot 8+2\cdot 1\cdot 9=50+32+18=100.$$
Для \(p=3\) получаем \(q_{3,1}=3\) и \(q_{3,2}=1\), значит
$$2\cdot 3\cdot 7+2\cdot 1\cdot 9=42+18=60.$$
Для \(p=5\) остаётся только \(q_{5,1}=2\), что даёт
$$2\cdot 2\cdot 8=32.$$
Для \(p=7\) остаётся только \(q_{7,1}=1\), что даёт
$$2\cdot 1\cdot 9=18.$$
Следовательно,
$$S(10)=100+60+32+18=210,$$
и это совпадает с малой проверкой, используемой в реализациях.
Реализации на C++, Python и Java устроены одинаково. Сначала вся арифметика ведётся по модулю \(10^9+7\), а каждый вклад сразу записывается в виде
$$2q(n-q)=2(nq-q^2).$$
Затем сегментированное решето генерирует все простые до порога \(B\), не требуя хранить полную таблицу простоты до \(10^8\). Для каждого малого простого программа последовательно умножает его само на себя, чтобы посетить \(p^1,p^2,p^3,\dots\) до тех пор, пока следующая степень не превысит \(n\).
После этого оставшиеся простые \(p\gt B\) обрабатываются блоками по одинаковому значению \(\left\lfloor n/p\right\rfloor\). Для каждого такого значения определяется соответствующий интервал простых, подсчитывается число простых внутри интервала и одним действием прибавляется нужное число одинаковых вкладов.
Эти интервальные подсчёты выполняются с помощью алгоритма подсчёта простых в стиле Лемера, основанного на небольшой предварительной решётке и мемоизации промежуточных результатов. Благодаря этому даже часть с большими простыми остаётся очень быстрой.
Пусть \(B=\min(n,10^8)\). Сегментированное решето для малых простых требует \(O(B\log\log B)\) времени и памяти порядка размера сегмента плюс небольшой базовой решётки. Перебор степеней малых простых добавляет по одному шагу на каждую степень \(p^k\le n\) с \(p\le B\).
Часть с большими простыми перебирает частные до \(n/(B+1)\). Для реального значения \(n=10^{12}\) это всего около \(10^4\) внешних итераций. На каждой из них нужны значения функции подсчёта простых на концах интервала, и мемоизированный метод Лемера делает эти запросы практичными.
В результате алгоритм на много порядков быстрее наивного \(O(n^2)\) перебора всех упорядоченных пар и без труда справляется с требуемым значением \(10^{12}\).
لعددين صحيحين موجبين \(a=\prod p^{\alpha_p}\) و\(b=\prod p^{\beta_p}\)، نعرّف فرق الأسس بالصيغة
$$D(a,b)=\sum_{p}\left|\alpha_p-\beta_p\right|=\sum_{p}\left|\nu_p(a)-\nu_p(b)\right|.$$
والمطلوب هو حساب
$$S(n)=\sum_{1\le a,b\le n} D(a,b)$$
في حالة Project Euler حيث \(n=10^{12}\)، ثم أخذ النتيجة النهائية بترديد \(10^9+7\). الفحص المباشر لكل الأزواج المرتبة مستحيل عمليًا، لأن عدد الأزواج وحده هائل، ناهيك عن تكرار التحليل إلى عوامل أولية. لذلك تعيد الحلول ترتيب المجموع بحسب كل عدد أولي وبحسب كل قوة أولية.
الفكرة الحاسمة هي أن مساهمة كل عدد أولي مستقلة عن غيره، وأن القيمة المطلقة لفارق أسين يمكن تحويلها إلى عدّ لمستويات القابلية للقسمة.
بما أن عددًا منتهيًا فقط من الأعداد الأولية يقسم \(a\) أو \(b\)، فيمكن تبديل ترتيب الجمع:
$$S(n)=\sum_{p\le n}\sum_{1\le a,b\le n}\left|\nu_p(a)-\nu_p(b)\right|.$$
ولعدد أولي ثابت \(p\)، نعرّف
$$T_p(n)=\sum_{1\le a,b\le n}\left|\nu_p(a)-\nu_p(b)\right|.$$
وعندئذ
$$S(n)=\sum_{p\le n} T_p(n).$$
وهكذا تُختزل المسألة كلها إلى فهم مساهمة عدد أولي واحد في كل مرة.
لكل عددين صحيحين غير سالبين \(x\) و\(y\) لدينا
$$|x-y|=\sum_{k\ge 1}\left|\mathbf{1}_{x\ge k}-\mathbf{1}_{y\ge k}\right|.$$
وعند تطبيق ذلك على \(x=\nu_p(a)\) و\(y=\nu_p(b)\)، فإن المستوى \(k\) يضيف \(1\) بالضبط عندما يكون أحد العددين قابلًا للقسمة على \(p^k\) والآخر غير قابل لها. لذلك
$$T_p(n)=\sum_{k\ge 1} N_{p,k},$$
حيث \(N_{p,k}\) هو عدد الأزواج المرتبة التي تختلف حالتها بالنسبة للقسمة على \(p^k\).
لنعرّف
$$q_{p,k}=\#\{m\le n: p^k\mid m\}=\left\lfloor\frac{n}{p^k}\right\rfloor.$$
أي أن هناك \(q_{p,k}\) عددًا حتى \(n\) تقبل القسمة على \(p^k\)، و\(n-q_{p,k}\) عددًا لا تقبلها. وبما أن الأزواج مرتبة، فإن عدد حالات الاختلاف يساوي
$$N_{p,k}=q_{p,k}(n-q_{p,k})+(n-q_{p,k})q_{p,k}=2q_{p,k}(n-q_{p,k}).$$
ومن ثم
$$T_p(n)=\sum_{k\ge 1} 2q_{p,k}(n-q_{p,k}),\qquad q_{p,k}=\left\lfloor\frac{n}{p^k}\right\rfloor.$$
وبمجرد أن تصبح \(p^k\gt n\)، يكون \(q_{p,k}=0\)، ولذلك فالمجموع منتهٍ.
بالتعويض نحصل على الصيغة
$$\boxed{S(n)=2\sum_{p\le n}\sum_{k\ge 1}\left\lfloor\frac{n}{p^k}\right\rfloor\left(n-\left\lfloor\frac{n}{p^k}\right\rfloor\right).}$$
الآن كل قوة أولية \(p^k\le n\) تعطي حدًا واحدًا يعتمد فقط على خارج القسمة \(\left\lfloor n/p^k\right\rfloor\). لذلك لا حاجة إلى فحص الأزواج \((a,b)\) نفسها، بل يكفي تعداد القوى الأولية بكفاءة أو تجميع كثير منها عندما تعطي خارج القسمة نفسه.
تختار الحلول الحدّ الفاصل
$$B=\min(n,10^8).$$
لكل عدد أولي \(p\le B\) تُعدّ جميع القوى \(p,p^2,p^3,\dots\le n\) مباشرة. أما في حالة المسألة الفعلية \(n=10^{12}\)، فإن أي عدد أولي \(p\gt B\) يحقق تلقائيًا \(p^2\gt n\)، وهذا يعني أنه لا يساهم إلا عبر القوة الأولى \(p\) فقط.
وهذه الملاحظة تزيل تمامًا الحاجة إلى التعامل مع القوى العليا في جزء الأعداد الأولية الكبيرة.
إذا كان \(p\gt B\)، فإن الكمية الوحيدة المهمة هي
$$q=\left\lfloor\frac{n}{p}\right\rfloor.$$
وكل الأعداد الأولية التي تعطي القيمة نفسها لـ \(q\) تقع في المجال
$$\frac{n}{q+1}\lt p\le \frac{n}{q}.$$
بدل المرور على هذه الأعداد الأولية واحدًا واحدًا، نعدّ كم عددًا أوليًا يوجد في ذلك المجال، ثم نضرب في المساهمة المشتركة
$$2q(n-q).$$
وعدد الأعداد الأولية في أي مجال \([L,R]\) يساوي
$$\pi(R)-\pi(L-1),$$
ولذلك يكفي امتلاك وسيلة سريعة لحساب دالة عدّ الأعداد الأولية \(\pi(x)\).
الأعداد الأولية حتى \(10\) هي \(2,3,5,7\).
بالنسبة إلى \(p=2\)، فإن القيم غير الصفرية هي \(q_{2,1}=5\) و\(q_{2,2}=2\) و\(q_{2,3}=1\)، ومن ثم تكون المساهمة
$$2\cdot 5\cdot 5+2\cdot 2\cdot 8+2\cdot 1\cdot 9=50+32+18=100.$$
وبالنسبة إلى \(p=3\)، نحصل على \(q_{3,1}=3\) و\(q_{3,2}=1\)، ولذلك
$$2\cdot 3\cdot 7+2\cdot 1\cdot 9=42+18=60.$$
أما \(p=5\) فلا يبقى إلا \(q_{5,1}=2\)، فتكون المساهمة
$$2\cdot 2\cdot 8=32.$$
وأما \(p=7\) فلا يبقى إلا \(q_{7,1}=1\)، فتكون المساهمة
$$2\cdot 1\cdot 9=18.$$
إذًا
$$S(10)=100+60+32+18=210,$$
وهو يطابق قيمة الاختبار الصغيرة المستخدمة في التطبيقات.
تتبع تطبيقات C++ وPython وJava المسار نفسه. أولًا، تُجرى كل العمليات بترديد \(10^9+7\)، ويُختزل كل حد مباشرة بصيغة
$$2q(n-q)=2(nq-q^2).$$
بعد ذلك، يولّد غربال مقطعي جميع الأعداد الأولية حتى الحد \(B\)، وبذلك لا حاجة إلى تخزين جدول كامل للأولية حتى \(10^8\). ولكل عدد أولي صغير، تُضرب القيمة في نفسها مرارًا لزيارة \(p^1,p^2,p^3,\dots\) حتى تصبح القوة التالية أكبر من \(n\).
ثم تُعالَج الأعداد الأولية الباقية \(p\gt B\) على شكل كتل تشترك في القيمة نفسها لـ \(\left\lfloor n/p\right\rfloor\). فتُحدَّد حدود المجال الموافق، ثم يُحسب عدد الأعداد الأولية داخل ذلك المجال، ثم تُضاف هذه المساهمة دفعة واحدة.
أما حساب أعداد الأوليات داخل المجالات فيعتمد على خوارزمية عدّ أوليات من نمط Lehmer، مدعومة بغربال صغير مُسبق وبالتخزين المؤقت للنتائج الوسيطة. وهذا ما يجعل جزء الأعداد الأولية الكبيرة سريعًا بما يكفي في الحالة المطلوبة.
إذا وضعنا \(B=\min(n,10^8)\)، فإن الغربال المقطعي للأعداد الأولية الصغيرة يحتاج زمنًا قدره \(O(B\log\log B)\)، بينما تكون الذاكرة في حدود حجم المقطع مع غربال أساس صغير. أما تعداد قوى الأعداد الأولية الصغيرة فيضيف خطوة واحدة لكل قوة أولية \(p^k\le n\) حيث \(p\le B\).
جزء الأعداد الأولية الكبيرة يمر على قيم خارج القسمة حتى \(n/(B+1)\). وبالنسبة إلى القيمة الفعلية \(n=10^{12}\)، فهذا يعني نحو \(10^4\) دورة خارجية فقط. وفي كل دورة تُجرى استعلامات عدّ أوليات عند طرفي المجال، وتتكفل طريقة Lehmer مع التخزين المؤقت بجعل هذه الاستعلامات عملية.
إجمالًا، هذه الطريقة أسرع بكثير جدًا من الفحص الساذج \(O(n^2)\) لكل الأزواج المرتبة، وهي مناسبة تمامًا للحجم المطلوب \(10^{12}\).