For each positive integer \(k\), define
$$S_k(N)=\sum_{n\le N}\operatorname{rad}(n)^k,$$
where \(\operatorname{rad}(n)\) is the product of the distinct prime divisors of \(n\). The problem asks for
$$\sum_{k=1}^{50} S_k(10^{12}) \pmod{10^9+7}.$$
The direct approach would require evaluating \(\operatorname{rad}(n)\) for every \(n\le 10^{12}\), which is hopeless. The implementation instead starts from the easier power sum \(\sum n^k\), then corrects it prime by prime until the local behavior matches \(\operatorname{rad}(n)^k\).
Fix one value of \(k\). The task is to compute \(S_k(N)\), and then repeat that computation for \(k=1,2,\dots,50\).
The multiplicative function \(n^k\) has local values
$$n^k=\prod_{p^a\parallel n} p^{ak}.$$
By contrast,
$$\operatorname{rad}(n)^k=\prod_{p^a\parallel n} p^k.$$
So the only difference is what happens when a prime appears with exponent \(a\ge 2\). If a prime occurs only once, both functions contribute the same factor \(p^k\).
The implementation therefore begins with the summatory function
$$P_k(x)=\sum_{n\le x} n^k,$$
which can be evaluated quickly for many different \(x\) by a Faulhaber polynomial modulo \(10^9+7\).
Suppose we have a temporary multiplicative weight in which the primes already processed behave like \(\operatorname{rad}(n)^k\), while the remaining primes still behave like \(n^k\). For the next prime \(p\), only prime powers \(p^a\) with \(a\ge 2\) need correction.
For one local factor we want to replace
$$p^{ak}\quad\text{by}\quad p^k \qquad (a\ge 2).$$
The difference is
$$p^{ak}-p^k=p^k\left(p^{(a-1)k}-1\right).$$
Using a geometric series, this becomes
$$p^{ak}-p^k=p^k(p^k-1)\left(1+p^k+p^{2k}+\cdots+p^{(a-2)k}\right).$$
Therefore the update coefficient for prime \(p\) is
$$t_p=p^k(p^k-1).$$
Let \(F(x)\) be the current summatory function before correcting prime \(p\). Then subtracting
$$t_p\sum_{a\ge 2} F\!\left(\left\lfloor\frac{x}{p^a}\right\rfloor\right)$$
exactly changes every local factor \(p^{ak}\) with \(a\ge 2\) into \(p^k\).
To see why, consider an integer \(n=p^e m\) with \(p\nmid m\). Its old local contribution at \(p\) is \(p^{ek}\). The subtraction contributes
$$t_p\sum_{a=2}^{e} p^{(e-a)k}=p^k(p^k-1)\sum_{j=0}^{e-2} p^{jk}=p^{ek}-p^k.$$
After subtraction, the new local contribution is therefore
$$p^{ek}-(p^{ek}-p^k)=p^k,$$
which is exactly what \(\operatorname{rad}(p^e)^k\) requires. If \(e=0\) or \(e=1\), no subtraction occurs, and the factor is already correct.
So the prime-by-prime transition is
$$F_{\text{new}}(x)=F_{\text{old}}(x)-t_p\sum_{a\ge 2} F_{\text{old}}\!\left(\left\lfloor\frac{x}{p^a}\right\rfloor\right).$$
If \(p>\sqrt{N}\), then \(p^2>N\). Such a prime can appear in any \(n\le N\) only with exponent \(0\) or \(1\). But for exponents \(0\) and \(1\), the weights \(n^k\) and \(\operatorname{rad}(n)^k\) already agree. Hence corrections are needed only for primes
$$p\le \sqrt{N}.$$
This is why the implementation sieves primes only up to \(\lfloor\sqrt{N}\rfloor\).
The recurrence never asks for arbitrary values of \(x\). Starting from \(N\), every transition divides by some \(p^a\) with \(a\ge 2\). After several corrections, every queried argument has the form
$$\left\lfloor\frac{N}{m}\right\rfloor,$$
where \(m\) is a powerful number, meaning that every prime exponent in \(m\) is either \(0\) or at least \(2\).
This is the crucial compression step. Instead of storing values for all \(1\le x\le N\), the implementation stores values only for the distinct quotients produced by powerful denominators. These quotients are sorted once, indexed once, and then reused for all fifty values of \(k\).
For each fixed \(k\), the algorithm starts from \(P_k(x)\), applies the prime corrections in increasing prime order, and obtains
$$S_k(N)=\sum_{n\le N}\operatorname{rad}(n)^k \pmod{10^9+7}.$$
The final answer is simply
$$\sum_{k=1}^{50} S_k(N)\pmod{10^9+7}.$$
For \(k=1\), the base sum is
$$P_1(10)=1+2+\cdots+10=55.$$
Only primes up to \(\sqrt{10}\) need correction, so only \(p=2\) and \(p=3\) matter.
For \(p=2\),
$$t_2=2(2-1)=2.$$
The correction uses \(2^2\) and \(2^3\):
$$2\left(P_1\!\left(\left\lfloor\frac{10}{4}\right\rfloor\right)+P_1\!\left(\left\lfloor\frac{10}{8}\right\rfloor\right)\right)=2\left(P_1(2)+P_1(1)\right)=2(3+1)=8.$$
So the total becomes \(55-8=47\).
For \(p=3\),
$$t_3=3(3-1)=6,$$
and only \(3^2\) contributes:
$$6\,P_1\!\left(\left\lfloor\frac{10}{9}\right\rfloor\right)=6\,P_1(1)=6.$$
Therefore
$$S_1(10)=55-8-6=41.$$
This matches the direct check
$$\operatorname{rad}(1),\dots,\operatorname{rad}(10)=1,2,3,2,5,6,7,2,3,10,$$
whose sum is \(41\).
The C++, Python, and Java implementations follow the same structure. First they generate every prime up to \(\lfloor\sqrt{N}\rfloor\). Then they precompute Faulhaber-style polynomial coefficients modulo \(10^9+7\), so that \(\sum_{n\le x} n^k\) can be evaluated quickly for any needed quotient \(x\).
Next they enumerate the powerful denominators that can arise in the recurrence. From those denominators they build the distinct quotient states \(\left\lfloor N/m\right\rfloor\), sort them in descending order, and create index tables so every future division lands on an already known state.
For each \(k\in\{1,\dots,50\}\), every stored state is initialized with the corresponding power sum \(P_k(x)\). The implementation then processes primes in increasing order. For a given prime \(p\), it computes the factor \(t_p=p^k(p^k-1)\), and for every state \(x\ge p^2\) it subtracts the indexed states at \(\left\lfloor x/p^2\right\rfloor,\left\lfloor x/p^3\right\rfloor,\dots\).
The states are updated from larger quotients to smaller quotients. That ordering is important: when the implementation needs the value at \(\left\lfloor x/p^a\right\rfloor\), it still reads the pre-update value for the current prime, which is exactly what the recurrence requires.
After all relevant primes have been processed, the state corresponding to \(x=N\) is \(S_k(N)\). The implementation adds that value to the running total and continues with the next \(k\).
Sieving primes up to \(\sqrt{N}\) costs \(O(\sqrt{N}\log\log N)\) time and \(O(\sqrt{N})\) memory. The Faulhaber preprocessing for \(k\le 50\) is tiny compared with the rest of the computation.
The main optimization is that the recurrence is evaluated only on compressed quotient states of the form \(\left\lfloor N/m\right\rfloor\) with powerful \(m\). That state set is far smaller than \(\{1,2,\dots,N\}\), which is what makes \(N=10^{12}\) feasible.
For each \(k\), the runtime is driven by the number of reachable pairs consisting of a quotient state and a prime power \(p^a\) with \(a\ge 2\). In other words, the algorithm scales with the compressed transition graph rather than with all integers up to \(N\). Memory is proportional to the number of stored quotient states together with the prime list and a few auxiliary tables.
Fur jede positive ganze Zahl \(k\) sei
$$S_k(N)=\sum_{n\le N}\operatorname{rad}(n)^k,$$
wobei \(\operatorname{rad}(n)\) das Produkt der verschiedenen Primteiler von \(n\) ist. Gesucht wird
$$\sum_{k=1}^{50} S_k(10^{12}) \pmod{10^9+7}.$$
Ein direkter Ansatz ware unmoglich, weil man fur alle \(n\le 10^{12}\) die Radikalfunktion auswerten musste. Die Implementierung beginnt deshalb mit der leichteren Potenzsumme \(\sum n^k\) und korrigiert diese dann Primzahl fur Primzahl, bis das lokale Verhalten genau \(\operatorname{rad}(n)^k\) entspricht.
Fixiere ein \(k\). Dann muss zuerst \(S_k(N)\) berechnet und dieser Vorgang anschliessend fur \(k=1,2,\dots,50\) wiederholt werden.
Die multiplikative Funktion \(n^k\) hat die lokalen Werte
$$n^k=\prod_{p^a\parallel n} p^{ak}.$$
Dagegen gilt
$$\operatorname{rad}(n)^k=\prod_{p^a\parallel n} p^k.$$
Der Unterschied tritt also nur dann auf, wenn ein Primteiler mit Exponent \(a\ge 2\) vorkommt. Bei Exponent \(1\) liefern beide Funktionen denselben Faktor \(p^k\).
Darum startet die Methode mit der Summenfunktion
$$P_k(x)=\sum_{n\le x} n^k,$$
die sich modulo \(10^9+7\) mithilfe eines Faulhaber-Polynoms schnell fur viele verschiedene Werte von \(x\) auswerten lasst.
Angenommen, wir besitzen vorubergehend eine multiplikative Gewichtung, bei der die bereits behandelten Primzahlen sich wie \(\operatorname{rad}(n)^k\) verhalten, wahrend die ubrigen Primzahlen noch wie \(n^k\) wirken. Fur die nachste Primzahl \(p\) mussen dann nur die Primzahlpotenzen \(p^a\) mit \(a\ge 2\) angepasst werden.
Lokalerweise soll
$$p^{ak}\quad\text{durch}\quad p^k \qquad (a\ge 2)$$
ersetzt werden. Die Differenz lautet
$$p^{ak}-p^k=p^k\left(p^{(a-1)k}-1\right).$$
Mit der geometrischen Reihe wird daraus
$$p^{ak}-p^k=p^k(p^k-1)\left(1+p^k+p^{2k}+\cdots+p^{(a-2)k}\right).$$
Der zugehorige Korrekturfaktor fur die Primzahl \(p\) ist also
$$t_p=p^k(p^k-1).$$
Sei \(F(x)\) die aktuelle Summenfunktion vor der Korrektur der Primzahl \(p\). Dann verandert das Subtrahieren von
$$t_p\sum_{a\ge 2} F\!\left(\left\lfloor\frac{x}{p^a}\right\rfloor\right)$$
genau jene lokalen Faktoren \(p^{ak}\) mit \(a\ge 2\) in den Zielwert \(p^k\).
Betrachte namlich \(n=p^e m\) mit \(p\nmid m\). Der alte lokale Beitrag bei \(p\) ist \(p^{ek}\). Die Subtraktion liefert
$$t_p\sum_{a=2}^{e} p^{(e-a)k}=p^k(p^k-1)\sum_{j=0}^{e-2} p^{jk}=p^{ek}-p^k.$$
Danach bleibt also
$$p^{ek}-(p^{ek}-p^k)=p^k,$$
genau der richtige lokale Faktor fur \(\operatorname{rad}(p^e)^k\). Fur \(e=0\) oder \(e=1\) gibt es keine Subtraktion; dort stimmt der Wert bereits.
Damit lautet der Primzahl-Ubergang
$$F_{\text{neu}}(x)=F_{\text{alt}}(x)-t_p\sum_{a\ge 2} F_{\text{alt}}\!\left(\left\lfloor\frac{x}{p^a}\right\rfloor\right).$$
Ist \(p>\sqrt{N}\), dann gilt \(p^2>N\). Eine solche Primzahl kann in einem \(n\le N\) nur mit Exponent \(0\) oder \(1\) auftreten. Fur diese Exponenten stimmen \(n^k\) und \(\operatorname{rad}(n)^k\) aber ohnehin uberein. Also mussen nur Primzahlen mit
$$p\le \sqrt{N}$$
korrigiert werden. Deshalb reicht ein Sieb bis \(\lfloor\sqrt{N}\rfloor\) vollstandig aus.
Die Rekursion fragt niemals beliebige Werte von \(x\) ab. Ausgehend von \(N\) wird in jedem Schritt durch \(p^a\) mit \(a\ge 2\) dividiert. Nach mehreren Schritten hat jedes abgefragte Argument die Form
$$\left\lfloor\frac{N}{m}\right\rfloor,$$
wobei \(m\) eine powerful number ist, also eine Zahl, in deren Primfaktorzerlegung jeder Exponent entweder \(0\) oder mindestens \(2\) ist.
Dies ist die zentrale Kompression. Anstatt Werte fur alle \(1\le x\le N\) zu speichern, halt die Implementierung nur die verschiedenen Quotienten fest, die aus solchen Nennern entstehen. Diese Zustande werden einmal sortiert, einmal indiziert und danach fur alle funfzig Werte von \(k\) wiederverwendet.
Fur jedes feste \(k\) startet das Verfahren bei \(P_k(x)\), verarbeitet die Primzahlen in aufsteigender Reihenfolge und erhalt schliesslich
$$S_k(N)=\sum_{n\le N}\operatorname{rad}(n)^k \pmod{10^9+7}.$$
Die gesuchte Endgroesse ist dann
$$\sum_{k=1}^{50} S_k(N)\pmod{10^9+7}.$$
Fur \(k=1\) ist die Ausgangssumme
$$P_1(10)=1+2+\cdots+10=55.$$
Nur Primzahlen bis \(\sqrt{10}\) mussen korrigiert werden, also nur \(p=2\) und \(p=3\).
Fur \(p=2\) gilt
$$t_2=2(2-1)=2.$$
Die Korrektur verwendet \(2^2\) und \(2^3\):
$$2\left(P_1\!\left(\left\lfloor\frac{10}{4}\right\rfloor\right)+P_1\!\left(\left\lfloor\frac{10}{8}\right\rfloor\right)\right)=2\left(P_1(2)+P_1(1)\right)=2(3+1)=8.$$
Danach steht die Summe bei \(55-8=47\).
Fur \(p=3\) ist
$$t_3=3(3-1)=6,$$
und nur \(3^2\) tragt bei:
$$6\,P_1\!\left(\left\lfloor\frac{10}{9}\right\rfloor\right)=6\,P_1(1)=6.$$
Somit folgt
$$S_1(10)=55-8-6=41.$$
Das stimmt mit der direkten Kontrolle uberein:
$$\operatorname{rad}(1),\dots,\operatorname{rad}(10)=1,2,3,2,5,6,7,2,3,10,$$
deren Summe ebenfalls \(41\) ist.
Die C++, Python- und Java-Implementierungen haben dieselbe Struktur. Zuerst erzeugen sie alle Primzahlen bis \(\lfloor\sqrt{N}\rfloor\). Danach werden Faulhaber-Koeffizienten modulo \(10^9+7\) vorbereitet, damit \(\sum_{n\le x} n^k\) fur jeden benotigten Quotienten \(x\) schnell ausgewertet werden kann.
Anschliessend enumeriert die Implementierung alle powerful Nenner, die in der Rekursion auftreten konnen. Daraus baut sie die verschiedenen Quotientenzustande \(\left\lfloor N/m\right\rfloor\), sortiert sie absteigend und legt Indextabellen an, sodass jede spatere Division sofort auf einen bereits bekannten Zustand zeigt.
Fur jedes \(k\in\{1,\dots,50\}\) wird jeder gespeicherte Zustand zunachst mit der Potenzsumme \(P_k(x)\) initialisiert. Danach werden die Primzahlen in aufsteigender Reihenfolge abgearbeitet. Fur eine feste Primzahl \(p\) berechnet die Implementierung den Faktor \(t_p=p^k(p^k-1)\) und subtrahiert fur jeden Zustand \(x\ge p^2\) die indizierten Zustande bei \(\left\lfloor x/p^2\right\rfloor,\left\lfloor x/p^3\right\rfloor,\dots\).
Die Aktualisierung lauft von grossen zu kleinen Quotienten. Genau diese Reihenfolge ist notwendig: Wenn der Wert bei \(\left\lfloor x/p^a\right\rfloor\) gebraucht wird, muss noch der alte Wert vor der aktuellen Primzahlkorrektur gelesen werden, und das wird durch die absteigende Ordnung sichergestellt.
Nachdem alle relevanten Primzahlen verarbeitet sind, ist der Zustand zu \(x=N\) gleich \(S_k(N)\). Dieser Wert wird zur Gesamtsumme addiert, und danach folgt das nachste \(k\).
Das Sieben der Primzahlen bis \(\sqrt{N}\) kostet \(O(\sqrt{N}\log\log N)\) Zeit und \(O(\sqrt{N})\) Speicher. Die Faulhaber-Vorbereitung fur \(k\le 50\) ist im Vergleich dazu sehr klein.
Die entscheidende Optimierung ist, dass die Rekursion nur auf komprimierten Quotientenzustanden der Form \(\left\lfloor N/m\right\rfloor\) mit powerful \(m\) ausgewertet wird. Diese Zustandsmenge ist viel kleiner als \(\{1,2,\dots,N\}\), und genau das macht \(N=10^{12}\) praktisch behandelbar.
Fur jedes \(k\) wird die Laufzeit durch die Zahl der erreichbaren Paare aus Quotientenzustand und Primzahlpotenz \(p^a\) mit \(a\ge 2\) bestimmt. Anders gesagt skaliert das Verfahren mit dem komprimierten Ubergangsgraphen und nicht mit allen ganzen Zahlen bis \(N\). Der Speicherverbrauch ist proportional zur Anzahl der gespeicherten Quotienten sowie zur Primzahlliste und einigen kleinen Hilfstabellen.
Her pozitif tamsayi \(k\) icin
$$S_k(N)=\sum_{n\le N}\operatorname{rad}(n)^k$$
tanimini yapalim. Burada \(\operatorname{rad}(n)\), \(n\)'in farkli asal carpanlarinin carpimidir. Istek su degerdir:
$$\sum_{k=1}^{50} S_k(10^{12}) \pmod{10^9+7}.$$
Dogrudan cozum, \(10^{12}\)'ye kadar her sayi icin \(\operatorname{rad}(n)\) hesaplamayi gerektirir ve bu uygulanabilir degildir. Uygulama bunun yerine once daha kolay olan \(\sum n^k\) toplamini kurar, sonra asal sayi sayi ilerleyerek yerel davranisi \(\operatorname{rad}(n)^k\) olacak sekilde duzeltir.
Bir \(k\) sabitleyelim. Once \(S_k(N)\) hesaplanir, sonra ayni surec \(k=1,2,\dots,50\) icin tekrarlanir.
\(n^k\) carpimsal fonksiyonunun asal kuvvetlerdeki yerel degeri
$$n^k=\prod_{p^a\parallel n} p^{ak}$$
seklindedir. Oysa
$$\operatorname{rad}(n)^k=\prod_{p^a\parallel n} p^k.$$
Yani fark sadece bir asalin ussu \(a\ge 2\) oldugunda ortaya cikar. Ussu \(1\) olan asal carpanlarda iki fonksiyon zaten ayni katsayiyi verir.
Bu nedenle uygulama once
$$P_k(x)=\sum_{n\le x} n^k$$
toplamini kullanir. Bu toplam, \(10^9+7\) modunda Faulhaber polinomlari yardimiyla cok sayida farkli \(x\) icin hizli bicimde hesaplanir.
Gecici olarak, islenmis asallarin \(\operatorname{rad}(n)^k\) gibi davrandigi, henuz islenmeyen asallarin ise \(n^k\) gibi kaldigi bir carpimsal agirlik dusunelim. Siradaki asal \(p\) icin sadece \(p^a\) ve \(a\ge 2\) durumlari duzeltilmelidir.
Tek bir yerel faktor icin amacimiz
$$p^{ak}\quad\text{yerine}\quad p^k \qquad (a\ge 2)$$
kullanmaktir. Aradaki fark
$$p^{ak}-p^k=p^k\left(p^{(a-1)k}-1\right)$$
olur. Geometrik seri ile bu ifade
$$p^{ak}-p^k=p^k(p^k-1)\left(1+p^k+p^{2k}+\cdots+p^{(a-2)k}\right)$$
haline gelir. Dolayisiyla asal \(p\) icin duzeltme katsayisi
$$t_p=p^k(p^k-1)$$
olur.
\(p\) asali icin duzeltme yapilmadan onceki toplami \(F(x)\) ile gosterelim. O zaman
$$t_p\sum_{a\ge 2} F\!\left(\left\lfloor\frac{x}{p^a}\right\rfloor\right)$$
ifadesini cikarmak, tam olarak \(a\ge 2\) olan yerel \(p^{ak}\) katkilarini \(p^k\)'ye indirger.
Bunu gormek icin \(n=p^e m\) ve \(p\nmid m\) alin. Eski yerel katkisi \(p^{ek}\)'dir. Cikartilan kisim
$$t_p\sum_{a=2}^{e} p^{(e-a)k}=p^k(p^k-1)\sum_{j=0}^{e-2} p^{jk}=p^{ek}-p^k$$
olur. Dolayisiyla yeni yerel katkisi
$$p^{ek}-(p^{ek}-p^k)=p^k$$
olur; bu da \(\operatorname{rad}(p^e)^k\) icin gereken dogru degerdir. \(e=0\) veya \(e=1\) oldugunda hic cikarma olmaz; zaten yerel deger dogrudur.
Boylece asal bazli gecis
$$F_{\text{yeni}}(x)=F_{\text{eski}}(x)-t_p\sum_{a\ge 2} F_{\text{eski}}\!\left(\left\lfloor\frac{x}{p^a}\right\rfloor\right)$$
seklini alir.
Eger \(p>\sqrt{N}\) ise \(p^2>N\) olur. Bu durumda \(p\), \(n\le N\) olan bir sayida sadece \(0\) veya \(1\) ussu ile bulunabilir. Bu uslerde ise \(n^k\) ile \(\operatorname{rad}(n)^k\) zaten ayni davranir. Dolayisiyla duzeltme yalnizca
$$p\le \sqrt{N}$$
olan asallar icin gereklidir. Uygulamanin \(\lfloor\sqrt{N}\rfloor\)'ye kadar elek kurmasinin nedeni budur.
Rekursiyon keyfi \(x\) degerleri istemez. Baslangic \(N\)'den itibaren her adimda \(a\ge 2\) olan bir \(p^a\) ile bolme yapilir. Birkac adim sonra sorgulanan her arguman su bicimdedir:
$$\left\lfloor\frac{N}{m}\right\rfloor,$$
burada \(m\) bir powerful number'dir; yani asal uslerinin her biri ya \(0\) ya da en az \(2\) olur.
Temel hizlandirma budur. \(1\) ile \(N\) arasindaki tum degerleri saklamak yerine uygulama, sadece bu powerful paydalardan gelen ayri bolumleri saklar. Bu bolumler bir kez siralanir, bir kez indekslenir ve sonra elli farkli \(k\) degeri icin tekrar tekrar kullanilir.
Her sabit \(k\) icin algoritma \(P_k(x)\) ile baslar, asallari kucukten buyuge dogru duzeltir ve sonunda
$$S_k(N)=\sum_{n\le N}\operatorname{rad}(n)^k \pmod{10^9+7}$$
degerini elde eder. Nihai cevap ise
$$\sum_{k=1}^{50} S_k(N)\pmod{10^9+7}$$
olur.
\(k=1\) icin baslangic toplami
$$P_1(10)=1+2+\cdots+10=55$$
olur. Sadece \(\sqrt{10}\)'a kadar olan asallari duzeltmek gerektiginden yalnizca \(p=2\) ve \(p=3\) gerekir.
\(p=2\) icin
$$t_2=2(2-1)=2$$
olur. Duzeltme, \(2^2\) ve \(2^3\) uzerinden gelir:
$$2\left(P_1\!\left(\left\lfloor\frac{10}{4}\right\rfloor\right)+P_1\!\left(\left\lfloor\frac{10}{8}\right\rfloor\right)\right)=2\left(P_1(2)+P_1(1)\right)=2(3+1)=8.$$
Toplam boylece \(55-8=47\) olur.
\(p=3\) icin
$$t_3=3(3-1)=6$$
ve sadece \(3^2\) katkisi vardir:
$$6\,P_1\!\left(\left\lfloor\frac{10}{9}\right\rfloor\right)=6\,P_1(1)=6.$$
Sonuc olarak
$$S_1(10)=55-8-6=41$$
elde edilir. Bu, dogrudan
$$\operatorname{rad}(1),\dots,\operatorname{rad}(10)=1,2,3,2,5,6,7,2,3,10$$
degerlerinin toplami ile de aynidir.
C++, Python ve Java uygulamalari ayni yapidandir. Once \(\lfloor\sqrt{N}\rfloor\)'ye kadar tum asallar uretilir. Ardindan Faulhaber tipinde polinom katsayilari \(10^9+7\) modunda hazirlanir; boylece \(\sum_{n\le x} n^k\) ifadesi gereken her bolum durumu icin hizla hesaplanabilir.
Daha sonra rekursiyonda ortaya cikabilecek powerful paydalar uretilir. Bunlardan \(\left\lfloor N/m\right\rfloor\) bicimindeki ayri durumlar kurulur, buyukten kucuge siralanir ve her biri icin indeks tablolar olusturulur. Boylece daha sonraki her bolme islemi daha once bilinen bir duruma baglanir.
Her \(k\in\{1,\dots,50\}\) icin butun durumlar once \(P_k(x)\) ile baslatilir. Sonra asallar artan sirayla islenir. Sabit bir asal \(p\) icin \(t_p=p^k(p^k-1)\) hesaplanir ve \(x\ge p^2\) olan her durumda \(\left\lfloor x/p^2\right\rfloor,\left\lfloor x/p^3\right\rfloor,\dots\) durumlarinin indekslenmis degerleri cikartilir.
Durumlar buyuk bolumlerden kucuk bolumlere dogru guncellenir. Bu sira zorunludur; cunku \(\left\lfloor x/p^a\right\rfloor\) degerine bakildiginda, o degerin mevcut asal icin henuz guncellenmemis eski surumunun okunmasi gerekir ve azalan sira bunu saglar.
Butun ilgili asallar bittiginde \(x=N\) durumundaki deger \(S_k(N)\) olur. Bu deger genel toplama eklenir ve sonraki \(k\) icin ayni surec tekrarlanir.
\(\sqrt{N}\)'ye kadar asal elemesi \(O(\sqrt{N}\log\log N)\) zaman ve \(O(\sqrt{N})\) bellek gerektirir. \(k\le 50\) icin Faulhaber hazirligi bunun yaninda cok kucuktur.
Asil optimizasyon, rekursiyonun yalnizca powerful \(m\) icin \(\left\lfloor N/m\right\rfloor\) bicimindeki sikistirilmis durumlarda calismasidir. Bu durum kumesi \(\{1,2,\dots,N\}\) kadar buyuk olmadigi icin \(N=10^{12}\) pratik hale gelir.
Her \(k\) icin calisma suresi, ulasilabilen bolum durumlari ile \(a\ge 2\) olan asal kuvvet gecislerinin sayisina baglidir. Yani algoritma tum sayilar uzerinden degil, sikistirilmis gecis grafigi uzerinden olceklenir. Bellek kullanimi da saklanan bolum sayisi, asal listesi ve kucuk yardimci tablolara orantilidir.
Para cada entero positivo \(k\), definimos
$$S_k(N)=\sum_{n\le N}\operatorname{rad}(n)^k,$$
donde \(\operatorname{rad}(n)\) es el producto de los divisores primos distintos de \(n\). La cantidad pedida es
$$\sum_{k=1}^{50} S_k(10^{12}) \pmod{10^9+7}.$$
Un metodo directo obligaria a evaluar \(\operatorname{rad}(n)\) para todos los \(n\le 10^{12}\), algo inviable. La implementacion parte en cambio de la suma mas sencilla \(\sum n^k\) y luego la corrige primo por primo hasta que el comportamiento local coincide con \(\operatorname{rad}(n)^k\).
Fijemos un valor de \(k\). Primero se calcula \(S_k(N)\), y despues se repite el mismo procedimiento para \(k=1,2,\dots,50\).
La funcion multiplicativa \(n^k\) tiene valores locales
$$n^k=\prod_{p^a\parallel n} p^{ak}.$$
En cambio,
$$\operatorname{rad}(n)^k=\prod_{p^a\parallel n} p^k.$$
La diferencia aparece solo cuando un primo tiene exponente \(a\ge 2\). Si un primo aparece una sola vez, ambas funciones aportan el mismo factor \(p^k\).
Por eso el algoritmo comienza con
$$P_k(x)=\sum_{n\le x} n^k,$$
que puede evaluarse rapidamente para muchos valores de \(x\) mediante un polinomio de Faulhaber modulo \(10^9+7\).
Supongamos que tenemos un peso multiplicativo temporal en el que los primos ya procesados se comportan como \(\operatorname{rad}(n)^k\), mientras que los primos aun no procesados siguen comportandose como \(n^k\). Para el siguiente primo \(p\), solo hay que corregir las potencias \(p^a\) con \(a\ge 2\).
Queremos reemplazar localmente
$$p^{ak}\quad\text{por}\quad p^k \qquad (a\ge 2).$$
La diferencia es
$$p^{ak}-p^k=p^k\left(p^{(a-1)k}-1\right).$$
Usando una serie geometrica, esto se reescribe como
$$p^{ak}-p^k=p^k(p^k-1)\left(1+p^k+p^{2k}+\cdots+p^{(a-2)k}\right).$$
Por tanto, el coeficiente de correccion asociado al primo \(p\) es
$$t_p=p^k(p^k-1).$$
Sea \(F(x)\) la suma acumulada actual antes de corregir el primo \(p\). Entonces restar
$$t_p\sum_{a\ge 2} F\!\left(\left\lfloor\frac{x}{p^a}\right\rfloor\right)$$
convierte exactamente cada factor local \(p^{ak}\) con \(a\ge 2\) en el valor deseado \(p^k\).
En efecto, si \(n=p^e m\) con \(p\nmid m\), la contribucion local antigua es \(p^{ek}\). La parte restada vale
$$t_p\sum_{a=2}^{e} p^{(e-a)k}=p^k(p^k-1)\sum_{j=0}^{e-2} p^{jk}=p^{ek}-p^k.$$
Despues de la correccion queda
$$p^{ek}-(p^{ek}-p^k)=p^k,$$
que es exactamente el factor local correcto para \(\operatorname{rad}(p^e)^k\). Si \(e=0\) o \(e=1\), no hay nada que restar y el valor ya era correcto.
La transicion primo a primo es entonces
$$F_{\text{nuevo}}(x)=F_{\text{viejo}}(x)-t_p\sum_{a\ge 2} F_{\text{viejo}}\!\left(\left\lfloor\frac{x}{p^a}\right\rfloor\right).$$
Si \(p>\sqrt{N}\), entonces \(p^2>N\). Un primo asi solo puede aparecer en un numero \(n\le N\) con exponente \(0\) o \(1\). Pero en esos casos \(n^k\) y \(\operatorname{rad}(n)^k\) ya coinciden. Por eso solo hay que corregir los primos que satisfacen
$$p\le \sqrt{N}.$$
Esa es la razon por la que la implementacion construye el cribado solo hasta \(\lfloor\sqrt{N}\rfloor\).
La recurrencia nunca pide valores arbitrarios de \(x\). Partiendo de \(N\), cada paso divide por algun \(p^a\) con \(a\ge 2\). Tras varias correcciones, cada argumento consultado tiene la forma
$$\left\lfloor\frac{N}{m}\right\rfloor,$$
donde \(m\) es un powerful number, es decir, un entero en cuya factorizacion cada exponente primo es \(0\) o al menos \(2\).
Aqui esta la compresion esencial. En vez de almacenar valores para todos los \(1\le x\le N\), la implementacion guarda solo los cocientes distintos producidos por esos denominadores poderosos. Esos estados se ordenan una vez, se indexan una vez y despues se reutilizan para los cincuenta valores de \(k\).
Para cada \(k\) fijo, el algoritmo comienza con \(P_k(x)\), aplica las correcciones primo por primo en orden creciente y obtiene al final
$$S_k(N)=\sum_{n\le N}\operatorname{rad}(n)^k \pmod{10^9+7}.$$
La respuesta final es
$$\sum_{k=1}^{50} S_k(N)\pmod{10^9+7}.$$
Para \(k=1\), la suma base es
$$P_1(10)=1+2+\cdots+10=55.$$
Solo hay que corregir los primos hasta \(\sqrt{10}\), asi que intervienen \(p=2\) y \(p=3\).
Para \(p=2\),
$$t_2=2(2-1)=2.$$
La correccion usa \(2^2\) y \(2^3\):
$$2\left(P_1\!\left(\left\lfloor\frac{10}{4}\right\rfloor\right)+P_1\!\left(\left\lfloor\frac{10}{8}\right\rfloor\right)\right)=2\left(P_1(2)+P_1(1)\right)=2(3+1)=8.$$
Por tanto, el total pasa a ser \(55-8=47\).
Para \(p=3\),
$$t_3=3(3-1)=6,$$
y solo contribuye \(3^2\):
$$6\,P_1\!\left(\left\lfloor\frac{10}{9}\right\rfloor\right)=6\,P_1(1)=6.$$
Asi obtenemos
$$S_1(10)=55-8-6=41.$$
Esto coincide con la comprobacion directa
$$\operatorname{rad}(1),\dots,\operatorname{rad}(10)=1,2,3,2,5,6,7,2,3,10,$$
cuya suma tambien es \(41\).
Las implementaciones en C++, Python y Java siguen exactamente la misma estructura. Primero generan todos los primos hasta \(\lfloor\sqrt{N}\rfloor\). Despues precalculan coeficientes tipo Faulhaber modulo \(10^9+7\), de modo que \(\sum_{n\le x} n^k\) pueda evaluarse rapidamente para cada cociente \(x\) que vaya apareciendo.
A continuacion enumeran los denominadores poderosos que pueden surgir en la recurrencia. A partir de ellos construyen los estados de cociente \(\left\lfloor N/m\right\rfloor\), los ordenan de mayor a menor y crean tablas de indices para que toda division posterior apunte directamente a un estado ya existente.
Para cada \(k\in\{1,\dots,50\}\), todos los estados almacenados se inicializan con la suma de potencias \(P_k(x)\). Luego se recorren los primos en orden creciente. Para un primo fijo \(p\), la implementacion calcula \(t_p=p^k(p^k-1)\) y, para cada estado \(x\ge p^2\), resta los estados indexados \(\left\lfloor x/p^2\right\rfloor,\left\lfloor x/p^3\right\rfloor,\dots\).
La actualizacion se hace desde cocientes grandes hacia cocientes pequenos. Esa direccion importa: cuando se necesita el valor en \(\left\lfloor x/p^a\right\rfloor\), todavia debe leerse la version anterior a la correccion actual del primo, y el orden descendente garantiza precisamente eso.
Cuando ya se han procesado todos los primos relevantes, el estado correspondiente a \(x=N\) es \(S_k(N)\). Ese valor se acumula en la suma global y el mismo procedimiento se repite para el siguiente \(k\).
El cribado de primos hasta \(\sqrt{N}\) cuesta \(O(\sqrt{N}\log\log N)\) tiempo y \(O(\sqrt{N})\) memoria. El preprocesamiento de Faulhaber para \(k\le 50\) es pequeno comparado con el resto.
La optimizacion decisiva es que la recurrencia se evalua solo sobre estados comprimidos de la forma \(\left\lfloor N/m\right\rfloor\) con \(m\) poderoso. Ese conjunto es muchisimo menor que \(\{1,2,\dots,N\}\), y por eso \(N=10^{12}\) se vuelve manejable.
Para cada \(k\), el tiempo de ejecucion depende del numero de pares alcanzables formados por un estado de cociente y una potencia prima \(p^a\) con \(a\ge 2\). En otras palabras, el algoritmo escala con el grafo de transiciones comprimido, no con todos los enteros hasta \(N\). La memoria es proporcional al numero de cocientes almacenados, junto con la lista de primos y unas pocas tablas auxiliares.
对每个正整数 \(k\)，定义
$$S_k(N)=\sum_{n\le N}\operatorname{rad}(n)^k,$$
其中 \(\operatorname{rad}(n)\) 表示 \(n\) 的所有不同素因子的乘积。题目要求计算
$$\sum_{k=1}^{50} S_k(10^{12}) \pmod{10^9+7}.$$
如果直接做，就必须对所有 \(n\le 10^{12}\) 逐个求 \(\operatorname{rad}(n)\)，这显然不可行。实现采用的思路是先从更容易的幂和 \(\sum n^k\) 出发，再按素数逐步修正，直到每个素数幂的局部贡献都变成 \(\operatorname{rad}(n)^k\) 所需要的形式。
先固定一个 \(k\)。目标是求出 \(S_k(N)\)，然后对 \(k=1,2,\dots,50\) 重复同样的流程并把结果累加。
函数 \(n^k\) 是乘法性的，而且在素数幂上满足
$$n^k=\prod_{p^a\parallel n} p^{ak}.$$
另一方面，
$$\operatorname{rad}(n)^k=\prod_{p^a\parallel n} p^k.$$
这说明两者的差别只发生在某个素数的指数 \(a\ge 2\) 时。若某个素数只出现一次，那么两种权重都贡献同一个因子 \(p^k\)。
因此算法首先建立
$$P_k(x)=\sum_{n\le x} n^k,$$
并利用 Faulhaber 多项式在模 \(10^9+7\) 下对大量不同的 \(x\) 快速求值。
设想当前有一个临时的乘法性权重：已经处理过的素数按照 \(\operatorname{rad}(n)^k\) 的规则计算，还没处理的素数仍然按照 \(n^k\) 的规则计算。对于下一个素数 \(p\)，真正需要修正的只有 \(p^a\) 且 \(a\ge 2\) 的情形。
我们希望把局部因子
$$p^{ak}\quad\text{替换为}\quad p^k \qquad (a\ge 2).$$
两者的差是
$$p^{ak}-p^k=p^k\left(p^{(a-1)k}-1\right).$$
把它写成等比级数后得到
$$p^{ak}-p^k=p^k(p^k-1)\left(1+p^k+p^{2k}+\cdots+p^{(a-2)k}\right).$$
因此，对素数 \(p\) 的修正常数正好是
$$t_p=p^k(p^k-1).$$
设 \(F(x)\) 表示修正素数 \(p\) 之前的当前前缀和函数。那么减去
$$t_p\sum_{a\ge 2} F\!\left(\left\lfloor\frac{x}{p^a}\right\rfloor\right)$$
就会把所有 \(a\ge 2\) 的局部贡献 \(p^{ak}\) 精确地改成目标值 \(p^k\)。
证明很直接。若 \(n=p^e m\) 且 \(p\nmid m\)，那么旧的局部贡献是 \(p^{ek}\)。减去的部分为
$$t_p\sum_{a=2}^{e} p^{(e-a)k}=p^k(p^k-1)\sum_{j=0}^{e-2} p^{jk}=p^{ek}-p^k.$$
于是修正后的局部贡献变成
$$p^{ek}-(p^{ek}-p^k)=p^k,$$
这恰好就是 \(\operatorname{rad}(p^e)^k\) 需要的值。如果 \(e=0\) 或 \(e=1\)，则根本不会发生减法，而原本的局部值已经是正确的。
因此，按素数进行更新的总公式是
$$F_{\text{new}}(x)=F_{\text{old}}(x)-t_p\sum_{a\ge 2} F_{\text{old}}\!\left(\left\lfloor\frac{x}{p^a}\right\rfloor\right).$$
如果 \(p>\sqrt{N}\)，那么 \(p^2>N\)。这样的素数在任何 \(n\le N\) 中只能以指数 \(0\) 或 \(1\) 出现。而在这两种指数下，\(n^k\) 与 \(\operatorname{rad}(n)^k\) 本来就完全一致。因此真正需要修正的只有
$$p\le \sqrt{N}$$
的素数。这也解释了为什么实现只筛到 \(\lfloor\sqrt{N}\rfloor\)。
递推过程中并不会访问任意的 \(x\)。从 \(N\) 出发，每一步都只会除以某个 \(p^a\)，其中 \(a\ge 2\)。经过若干步之后，所有被访问到的参数都具有形式
$$\left\lfloor\frac{N}{m}\right\rfloor,$$
其中 \(m\) 是一个 powerful number，也就是在素因子分解中，每个指数不是 \(0\)，就是至少 \(2\)。
这一步是整个算法能跑起来的关键。实现并不保存从 \(1\) 到 \(N\) 的所有状态，而只保存由这些 powerful 分母产生的不同商值。所有这些商值只需排序一次、编号一次，之后就可以在全部五十个 \(k\) 上重复使用。
对于固定的 \(k\)，算法从 \(P_k(x)\) 开始，按照素数从小到大的顺序依次修正，最后得到
$$S_k(N)=\sum_{n\le N}\operatorname{rad}(n)^k \pmod{10^9+7}.$$
最终答案就是
$$\sum_{k=1}^{50} S_k(N)\pmod{10^9+7}.$$
当 \(k=1\) 时，基础幂和为
$$P_1(10)=1+2+\cdots+10=55.$$
由于只需要处理不超过 \(\sqrt{10}\) 的素数，所以这里只涉及 \(p=2\) 和 \(p=3\)。
对 \(p=2\)，有
$$t_2=2(2-1)=2.$$
对应的修正来自 \(2^2\) 和 \(2^3\)：
$$2\left(P_1\!\left(\left\lfloor\frac{10}{4}\right\rfloor\right)+P_1\!\left(\left\lfloor\frac{10}{8}\right\rfloor\right)\right)=2\left(P_1(2)+P_1(1)\right)=2(3+1)=8.$$
因此总和先变成 \(55-8=47\)。
对 \(p=3\)，
$$t_3=3(3-1)=6,$$
只需要处理 \(3^2\)：
$$6\,P_1\!\left(\left\lfloor\frac{10}{9}\right\rfloor\right)=6\,P_1(1)=6.$$
于是
$$S_1(10)=55-8-6=41.$$
这和直接计算
$$\operatorname{rad}(1),\dots,\operatorname{rad}(10)=1,2,3,2,5,6,7,2,3,10$$
再求和得到的 \(41\) 完全一致。
C++、Python 和 Java 三个实现的结构完全一致。首先筛出所有不超过 \(\lfloor\sqrt{N}\rfloor\) 的素数。然后预处理 Faulhaber 型多项式系数，使得 \(\sum_{n\le x} n^k\) 能够在模 \(10^9+7\) 下对任意需要的商值 \(x\) 快速求出。
接着，程序枚举递推中可能出现的 powerful 分母，并据此构造所有不同的商状态 \(\left\lfloor N/m\right\rfloor\)。这些状态会按从大到小排序，再建立索引表，从而保证后续每一次除法都能直接落到某个已经存在的状态上。
对于每个 \(k\in\{1,\dots,50\}\)，所有保存的状态先用幂和 \(P_k(x)\) 初始化。之后程序按素数递增的顺序处理。对某个固定素数 \(p\)，先计算 \(t_p=p^k(p^k-1)\)，然后对每个满足 \(x\ge p^2\) 的状态，减去 \(\left\lfloor x/p^2\right\rfloor,\left\lfloor x/p^3\right\rfloor,\dots\) 这些已编号状态所对应的值。
状态更新的顺序是从较大的商到较小的商。这一点非常关键：当程序需要读取 \(\left\lfloor x/p^a\right\rfloor\) 的值时，必须读到当前素数修正之前的旧值，而按降序更新正好保证了这一点。
当所有相关素数都处理完以后，与 \(x=N\) 对应的状态就是 \(S_k(N)\)。实现把这个值加入总答案，然后继续处理下一个 \(k\)。
筛到 \(\sqrt{N}\) 的素数需要 \(O(\sqrt{N}\log\log N)\) 时间和 \(O(\sqrt{N})\) 空间。对于 \(k\le 50\) 的 Faulhaber 预处理，相比主过程可以看作很小。
真正的优化来自状态压缩：递推只在 \(\left\lfloor N/m\right\rfloor\) 这类状态上进行，其中 \(m\) 是 powerful 数。这样的状态集合远远小于 \(\{1,2,\dots,N\}\)，这正是 \(N=10^{12}\) 仍然可计算的原因。
对每个 \(k\) 而言，运行时间由可达的状态数以及所有 \(a\ge 2\) 的素数幂转移共同决定。换句话说，算法的规模取决于压缩后的转移图，而不是取决于 \(1\) 到 \(N\) 的全部整数。内存使用量则与保存的商状态数量、素数表以及少量辅助表成正比。
Для каждого положительного целого \(k\) определим
$$S_k(N)=\sum_{n\le N}\operatorname{rad}(n)^k,$$
где \(\operatorname{rad}(n)\) есть произведение различных простых делителей числа \(n\). Требуется найти
$$\sum_{k=1}^{50} S_k(10^{12}) \pmod{10^9+7}.$$
Прямой перебор невозможен: пришлось бы вычислять \(\operatorname{rad}(n)\) для всех \(n\le 10^{12}\). Поэтому реализация начинает с гораздо более простой суммы \(\sum n^k\), а затем исправляет ее по простым числам, пока локальное поведение не станет ровно таким, как у \(\operatorname{rad}(n)^k\).
Зафиксируем одно значение \(k\). Нужно вычислить \(S_k(N)\), а затем повторить тот же процесс для \(k=1,2,\dots,50\).
Мультипликативная функция \(n^k\) имеет локальные значения
$$n^k=\prod_{p^a\parallel n} p^{ak}.$$
С другой стороны,
$$\operatorname{rad}(n)^k=\prod_{p^a\parallel n} p^k.$$
Значит, различие появляется только тогда, когда некоторый простой делитель входит с показателем \(a\ge 2\). Если показатель равен \(1\), обе функции дают один и тот же множитель \(p^k\).
Поэтому алгоритм стартует с сумматорной функции
$$P_k(x)=\sum_{n\le x} n^k,$$
которую можно быстро вычислять для многих значений \(x\) с помощью полинома Фаульхабера по модулю \(10^9+7\).
Предположим, что у нас есть временный мультипликативный вес, в котором уже обработанные простые ведут себя как \(\operatorname{rad}(n)^k\), а остальные простые пока еще ведут себя как \(n^k\). Для следующего простого \(p\) нужно исправить только случаи \(p^a\) при \(a\ge 2\).
Локально мы хотим заменить
$$p^{ak}\quad\text{на}\quad p^k \qquad (a\ge 2).$$
Разность равна
$$p^{ak}-p^k=p^k\left(p^{(a-1)k}-1\right).$$
Через геометрическую прогрессию получаем
$$p^{ak}-p^k=p^k(p^k-1)\left(1+p^k+p^{2k}+\cdots+p^{(a-2)k}\right).$$
Значит, коэффициент поправки для простого \(p\) равен
$$t_p=p^k(p^k-1).$$
Пусть \(F(x)\) обозначает текущую сумматорную функцию до обработки простого \(p\). Тогда вычитание
$$t_p\sum_{a\ge 2} F\!\left(\left\lfloor\frac{x}{p^a}\right\rfloor\right)$$
ровно и переводит каждый локальный вклад \(p^{ak}\) при \(a\ge 2\) в нужное значение \(p^k\).
Действительно, для числа \(n=p^e m\) при \(p\nmid m\) старый локальный вклад равен \(p^{ek}\). Вычитаемая часть равна
$$t_p\sum_{a=2}^{e} p^{(e-a)k}=p^k(p^k-1)\sum_{j=0}^{e-2} p^{jk}=p^{ek}-p^k.$$
После вычитания остается
$$p^{ek}-(p^{ek}-p^k)=p^k,$$
то есть именно тот локальный множитель, который нужен для \(\operatorname{rad}(p^e)^k\). Если \(e=0\) или \(e=1\), вычитания нет, и значение уже было правильным.
Итак, переход по простому имеет вид
$$F_{\text{new}}(x)=F_{\text{old}}(x)-t_p\sum_{a\ge 2} F_{\text{old}}\!\left(\left\lfloor\frac{x}{p^a}\right\rfloor\right).$$
Если \(p>\sqrt{N}\), то \(p^2>N\). Такой простой может входить в любое число \(n\le N\) только с показателем \(0\) или \(1\). Но в этих случаях веса \(n^k\) и \(\operatorname{rad}(n)^k\) уже совпадают. Следовательно, исправлять нужно только простые
$$p\le \sqrt{N}.$$
Именно поэтому в реализации простые просеиваются только до \(\lfloor\sqrt{N}\rfloor\).
Рекурсия никогда не обращается к произвольным значениям \(x\). Начиная с \(N\), каждый шаг делит на некоторое \(p^a\) с \(a\ge 2\). После нескольких шагов любой запрашиваемый аргумент имеет вид
$$\left\lfloor\frac{N}{m}\right\rfloor,$$
где \(m\) является powerful number, то есть числом, в разложении которого каждый простой показатель либо равен \(0\), либо не меньше \(2\).
Здесь и происходит ключевое сжатие. Вместо хранения значений для всех \(1\le x\le N\) реализация хранит только различные частные, возникающие из таких знаменателей. Эти состояния один раз сортируются, один раз индексируются, а затем повторно используются для всех пятидесяти значений \(k\).
Для каждого фиксированного \(k\) алгоритм стартует с \(P_k(x)\), обрабатывает простые в возрастающем порядке и в итоге получает
$$S_k(N)=\sum_{n\le N}\operatorname{rad}(n)^k \pmod{10^9+7}.$$
Окончательный ответ равен
$$\sum_{k=1}^{50} S_k(N)\pmod{10^9+7}.$$
При \(k=1\) базовая сумма равна
$$P_1(10)=1+2+\cdots+10=55.$$
Исправлять нужно только простые до \(\sqrt{10}\), то есть \(p=2\) и \(p=3\).
Для \(p=2\)
$$t_2=2(2-1)=2.$$
Поправка использует \(2^2\) и \(2^3\):
$$2\left(P_1\!\left(\left\lfloor\frac{10}{4}\right\rfloor\right)+P_1\!\left(\left\lfloor\frac{10}{8}\right\rfloor\right)\right)=2\left(P_1(2)+P_1(1)\right)=2(3+1)=8.$$
После этого получаем \(55-8=47\).
Для \(p=3\)
$$t_3=3(3-1)=6,$$
и вклад дает только \(3^2\):
$$6\,P_1\!\left(\left\lfloor\frac{10}{9}\right\rfloor\right)=6\,P_1(1)=6.$$
Следовательно,
$$S_1(10)=55-8-6=41.$$
Это совпадает с прямой проверкой
$$\operatorname{rad}(1),\dots,\operatorname{rad}(10)=1,2,3,2,5,6,7,2,3,10,$$
сумма которых действительно равна \(41\).
Реализации на C++, Python и Java устроены одинаково. Сначала они строят список всех простых чисел до \(\lfloor\sqrt{N}\rfloor\). Затем заранее вычисляют коэффициенты типа Фаульхабера по модулю \(10^9+7\), чтобы быстро получать \(\sum_{n\le x} n^k\) для любого нужного частного \(x\).
После этого программа перечисляет все powerful знаменатели, которые могут возникнуть в рекурсии. Из них строятся различные состояния вида \(\left\lfloor N/m\right\rfloor\); они сортируются по убыванию и снабжаются индексами, чтобы каждое последующее деление сразу переходило к уже известному состоянию.
Для каждого \(k\in\{1,\dots,50\}\) все сохраненные состояния сначала инициализируются значениями \(P_k(x)\). Затем простые числа обрабатываются в возрастающем порядке. Для данного простого \(p\) вычисляется \(t_p=p^k(p^k-1)\), и для каждого состояния \(x\ge p^2\) вычитаются уже индексированные состояния \(\left\lfloor x/p^2\right\rfloor,\left\lfloor x/p^3\right\rfloor,\dots\).
Обновление идет от больших частных к меньшим. Это принципиально важно: когда требуется значение в \(\left\lfloor x/p^a\right\rfloor\), нужно читать старую версию до текущей поправки по простому, и убывающий порядок как раз это гарантирует.
После обработки всех relevant простых состояние, соответствующее \(x=N\), уже равно \(S_k(N)\). Это значение добавляется к общей сумме, и затем выполняется следующий проход для нового \(k\).
Просеивание простых до \(\sqrt{N}\) требует \(O(\sqrt{N}\log\log N)\) времени и \(O(\sqrt{N})\) памяти. Предварительная подготовка Faulhaber для \(k\le 50\) по сравнению с основным этапом очень мала.
Главная оптимизация состоит в том, что рекурсия считается только на сжатых состояниях вида \(\left\lfloor N/m\right\rfloor\), где \(m\) является powerful числом. Это множество состояний намного меньше, чем \(\{1,2,\dots,N\}\), и именно поэтому \(N=10^{12}\) становится достижимым.
Для каждого \(k\) время работы определяется числом достижимых пар: состояние частного плюс простая степень \(p^a\) при \(a\ge 2\). Иначе говоря, масштаб алгоритма задается сжатым графом переходов, а не всеми целыми числами до \(N\). Память пропорциональна числу сохраненных частных, списку простых и нескольким небольшим вспомогательным таблицам.
لكل عدد صحيح موجب \(k\)، نعرّف
$$S_k(N)=\sum_{n\le N}\operatorname{rad}(n)^k,$$
حيث تمثل \(\operatorname{rad}(n)\) حاصل ضرب العوامل الاولية المميزة للعدد \(n\). المطلوب هو حساب
$$\sum_{k=1}^{50} S_k(10^{12}) \pmod{10^9+7}.$$
الحل المباشر غير ممكن عمليا، لانه يتطلب حساب \(\operatorname{rad}(n)\) لكل \(n\le 10^{12}\). لذلك تبدا الخوارزمية من مجموع اسهل هو \(\sum n^k\)، ثم تصححه اوليا بعد اولي حتى يصبح السلوك المحلي مطابقا تماما لـ \(\operatorname{rad}(n)^k\).
لنثبت قيمة واحدة لـ \(k\). الهدف اولا هو حساب \(S_k(N)\)، ثم تكرار الفكرة نفسها لكل \(k=1,2,\dots,50\).
الدالة الضربية \(n^k\) تملك القيم المحلية
$$n^k=\prod_{p^a\parallel n} p^{ak}.$$
بينما
$$\operatorname{rad}(n)^k=\prod_{p^a\parallel n} p^k.$$
لذلك لا يظهر الفرق الا عندما يكون اس العوامل الاولية \(a\ge 2\). اما اذا ظهر العامل الاولي مرة واحدة فقط، فان الدالتين تعطيان العامل نفسه \(p^k\).
لهذا السبب تبدا الطريقة من
$$P_k(x)=\sum_{n\le x} n^k,$$
وهو مجموع يمكن تقييمه بسرعة لكثير من قيم \(x\) باستخدام كثير حدود فاولهابر بت合同 \(10^9+7\).
افترض ان لدينا وزنا ضربيًا مؤقتا: الاوليات التي عولجت سابقا تتصرف كما في \(\operatorname{rad}(n)^k\)، بينما الاوليات التي لم تعالج بعد ما زالت تتصرف كما في \(n^k\). عند الانتقال الى اولي جديد \(p\)، لا نحتاج الى تصحيح الا الحالات التي يظهر فيها \(p^a\) مع \(a\ge 2\).
محليا نريد استبدال
$$p^{ak}\quad\to\quad p^k \qquad (a\ge 2).$$
والفرق بينهما هو
$$p^{ak}-p^k=p^k\left(p^{(a-1)k}-1\right).$$
وباستخدام متسلسلة هندسية نحصل على
$$p^{ak}-p^k=p^k(p^k-1)\left(1+p^k+p^{2k}+\cdots+p^{(a-2)k}\right).$$
ومن ثم يكون معامل التصحيح المرتبط بالاولي \(p\) هو
$$t_p=p^k(p^k-1).$$
لنرمز الى دالة المجموع الحالية قبل تصحيح الاولي \(p\) بالرمز \(F(x)\). عندئذ فان طرح
$$t_p\sum_{a\ge 2} F\!\left(\left\lfloor\frac{x}{p^a}\right\rfloor\right)$$
يحوّل بالضبط كل مساهمة محلية من الشكل \(p^{ak}\) مع \(a\ge 2\) الى القيمة المطلوبة \(p^k\).
لرؤية ذلك، خذ \(n=p^e m\) حيث \(p\nmid m\). المساهمة المحلية القديمة هي \(p^{ek}\). اما الجزء المطروح فهو
$$t_p\sum_{a=2}^{e} p^{(e-a)k}=p^k(p^k-1)\sum_{j=0}^{e-2} p^{jk}=p^{ek}-p^k.$$
وبالتالي تصبح المساهمة الجديدة
$$p^{ek}-(p^{ek}-p^k)=p^k,$$
وهي بالضبط القيمة المحلية الصحيحة لـ \(\operatorname{rad}(p^e)^k\). اما اذا كان \(e=0\) او \(e=1\)، فلا يحدث اي طرح اصلا، لان القيمة تكون صحيحة من البداية.
اذن تكون خطوة الانتقال الخاصة بكل اولي
$$F_{\text{new}}(x)=F_{\text{old}}(x)-t_p\sum_{a\ge 2} F_{\text{old}}\!\left(\left\lfloor\frac{x}{p^a}\right\rfloor\right).$$
اذا كان \(p>\sqrt{N}\)، فحينها \(p^2>N\). وهذا يعني ان الاولي \(p\) لا يمكن ان يظهر في اي \(n\le N\) الا باس \(0\) او \(1\). وفي هاتين الحالتين تتطابق قيمتا \(n^k\) و \(\operatorname{rad}(n)^k\) مسبقا. لذلك لا حاجة الى التصحيح الا للاوليات التي تحقق
$$p\le \sqrt{N}.$$
ولهذا تكتفي الخوارزمية بانتاج الاوليات حتى \(\lfloor\sqrt{N}\rfloor\).
العلاقة التكرارية لا تحتاج الى قيم عشوائية لـ \(x\). انطلاقا من \(N\)، كل انتقال يقسم على بعض \(p^a\) مع \(a\ge 2\). وبعد عدة خطوات، يصبح كل وسيط مطلوب من الشكل
$$\left\lfloor\frac{N}{m}\right\rfloor,$$
حيث يكون \(m\) من نوع powerful number، اي ان كل اس اولي فيه اما \(0\) او على الاقل \(2\).
هذه هي فكرة الضغط الحاسمة. بدلا من تخزين القيم لكل \(1\le x\le N\)، تخزن الخوارزمية فقط القيم المختلفة الناتجة من هذه المقامات القوية. ثم ترتب هذه الحالات مرة واحدة، وتفهرس مرة واحدة، وتستعمل من جديد لكل القيم الخمسين لـ \(k\).
لكل \(k\) ثابت، تبدا الخوارزمية من \(P_k(x)\)، ثم تطبق التصحيحات على الاوليات بترتيب تصاعدي، وفي النهاية نحصل على
$$S_k(N)=\sum_{n\le N}\operatorname{rad}(n)^k \pmod{10^9+7}.$$
والجواب النهائي هو
$$\sum_{k=1}^{50} S_k(N)\pmod{10^9+7}.$$
عند \(k=1\)، يكون مجموع الاساس هو
$$P_1(10)=1+2+\cdots+10=55.$$
وبما اننا لا نحتاج الا الى الاوليات حتى \(\sqrt{10}\)، فالاوليان الوحيدان المؤثران هنا هما \(2\) و\(3\).
بالنسبة الى \(p=2\)،
$$t_2=2(2-1)=2.$$
وياتي التصحيح من \(2^2\) و\(2^3\):
$$2\left(P_1\!\left(\left\lfloor\frac{10}{4}\right\rfloor\right)+P_1\!\left(\left\lfloor\frac{10}{8}\right\rfloor\right)\right)=2\left(P_1(2)+P_1(1)\right)=2(3+1)=8.$$
ومن ثم يصبح المجموع \(55-8=47\).
اما بالنسبة الى \(p=3\)،
$$t_3=3(3-1)=6,$$
ولا يساهم الا \(3^2\):
$$6\,P_1\!\left(\left\lfloor\frac{10}{9}\right\rfloor\right)=6\,P_1(1)=6.$$
اذن
$$S_1(10)=55-8-6=41.$$
وهذا يطابق الفحص المباشر للقيم
$$\operatorname{rad}(1),\dots,\operatorname{rad}(10)=1,2,3,2,5,6,7,2,3,10,$$
ومجموعها بالفعل \(41\).
تتبع تطبيقات C++ وPython وJava البنية نفسها. اولا يتم توليد جميع الاوليات حتى \(\lfloor\sqrt{N}\rfloor\). بعد ذلك تُحضَّر معاملات على نمط فاولهابر بت合同 \(10^9+7\)، بحيث يمكن حساب \(\sum_{n\le x} n^k\) بسرعة لكل خارج قسمة مطلوب.
ثم تقوم الخوارزمية بتعداد المقامات من نوع powerful التي يمكن ان تظهر داخل العلاقة التكرارية. ومن هذه المقامات تبني الحالات المختلفة من الشكل \(\left\lfloor N/m\right\rfloor\)، وترتبها تنازليا، وتجهز جداول فهرسة بحيث تشير كل عملية قسمة لاحقة مباشرة الى حالة معروفة مسبقا.
لكل \(k\in\{1,\dots,50\}\)، تهيأ جميع الحالات المحفوظة بقيم \(P_k(x)\). بعد ذلك تُعالج الاوليات بترتيب تصاعدي. بالنسبة الى اولي معين \(p\)، تحسب الخوارزمية \(t_p=p^k(p^k-1)\)، ثم تطرح من كل حالة تحقق \(x\ge p^2\) القيم المفهرسة المقابلة لـ \(\left\lfloor x/p^2\right\rfloor,\left\lfloor x/p^3\right\rfloor,\dots\).
ويتم التحديث من الخارجات الكبيرة الى الخارجات الصغيرة. هذا الترتيب ضروري، لانه عندما تحتاج الخوارزمية الى القيمة عند \(\left\lfloor x/p^a\right\rfloor\)، يجب ان تقرأ النسخة القديمة السابقة لتصحيح الاولي الحالي، والترتيب التنازلي يحقق ذلك تماما.
بعد الانتهاء من جميع الاوليات اللازمة، تصبح الحالة الموافقة لـ \(x=N\) مساوية لـ \(S_k(N)\). وتُضاف هذه القيمة الى المجموع الكلي، ثم يعاد المسار كله للقيمة التالية من \(k\).
توليد الاوليات حتى \(\sqrt{N}\) يكلف \(O(\sqrt{N}\log\log N)\) زمنا و\(O(\sqrt{N})\) ذاكرة. اما التحضير الخاص بصيغ فاولهابر حتى \(k\le 50\) فهو صغير جدا مقارنة بباقي العمل.
التحسين الجوهري هو ان العلاقة التكرارية لا تُقيَّم الا على حالات مضغوطة من الشكل \(\left\lfloor N/m\right\rfloor\) حيث \(m\) عدد powerful. هذه المجموعة اصغر بكثير من \(\{1,2,\dots,N\}\)، وهذا بالتحديد ما يجعل \(N=10^{12}\) قابلا للحساب.
بالنسبة الى كل \(k\)، يحدد زمن التنفيذ عدد الازواج الممكنة من حالة خارج قسمة مع قوة اولية \(p^a\) حيث \(a\ge 2\). بعبارة اخرى، الخوارزمية تتوسع وفق رسم الانتقالات المضغوط، لا وفق جميع الاعداد حتى \(N\). والذاكرة تتناسب مع عدد الحالات المخزنة، بالاضافة الى قائمة الاوليات وبعض الجداول الصغيرة المساعدة.