For an integer \(k\), define the ordinary power sum and its cumulative sum by
$$F_k(n)=\sum_{j=1}^{n} j^k,\qquad S_k(n)=\sum_{i=1}^{n}F_k(i).$$
The concrete task is to evaluate
$$\sum_{p\in \mathcal{P}} \left(S_{10000}(10^{12}) \bmod p\right),\qquad \mathcal{P}=\{p\text{ prime}: 2\cdot 10^9 \le p \le 2\cdot 10^9+2000\}.$$
A direct computation is hopeless: \(n=10^{12}\) is far too large for nested summation, and even evaluating \(F_k(n)\) term by term would be impossible. The solution therefore transforms the double sum into polynomial evaluation modulo each prime.
The implementations solve the problem prime by prime. For each prime \(p\) in the interval, the goal is to obtain \(S_k(n)\bmod p\) without summing up to \(n\).
Start from the definition
$$S_k(n)=\sum_{i=1}^{n}\sum_{j=1}^{i} j^k.$$
Swap the order of summation. A fixed value \(j\) appears in every inner sum with \(i\ge j\), so it is counted exactly \(n+1-j\) times. Hence
$$S_k(n)=\sum_{j=1}^{n}(n+1-j)j^k.$$
Now separate the two terms:
$$S_k(n)=(n+1)\sum_{j=1}^{n}j^k-\sum_{j=1}^{n}j^{k+1}=(n+1)F_k(n)-F_{k+1}(n).$$
So the whole problem reduces to evaluating two ordinary power sums modulo \(p\): \(F_k(n)\) and \(F_{k+1}(n)\).
For a fixed exponent \(e\), the function
$$F_e(x)=\sum_{t=1}^{x} t^e$$
is a polynomial in \(x\) of degree \(e+1\). One way to see this is through the forward difference
$$F_e(x)-F_e(x-1)=x^e,$$
whose right-hand side has degree \(e\), so \(F_e\) must have degree \(e+1\).
That means \(F_e(x)\) is uniquely determined by its values at any \(e+2\) distinct points. Over the field \(\mathbb{F}_p\), the implementations use the points
$$x=0,1,2,\dots,e+1.$$
This is valid because every prime in the target interval is much larger than \(e+1\), so these nodes remain distinct modulo \(p\).
Let \(m=e+1\). Define
$$y_i=F_e(i)\bmod p\qquad (0\le i\le m).$$
These values are easy to generate iteratively:
$$y_0=0,\qquad y_i=y_{i-1}+i^e \pmod p.$$
So for each prime \(p\), the implementation computes \(i^e\bmod p\) for \(i=1,2,\dots,m\), accumulates those values, and obtains the full interpolation table \((0,y_0),(1,y_1),\dots,(m,y_m)\).
No Bernoulli numbers or closed forms are needed. The only ingredients are modular exponentiation and cumulative addition.
Given the samples \(y_0,\dots,y_m\), Lagrange interpolation reconstructs the value at any \(x\in\mathbb{F}_p\):
$$F_e(x)=\sum_{i=0}^{m} y_i \prod_{\substack{0\le j\le m\\j\ne i}}\frac{x-j}{i-j}\pmod p.$$
The denominator can be simplified explicitly:
$$\prod_{\substack{0\le j\le m\\j\ne i}}(i-j)=i!\,(-1)^{m-i}(m-i)!.$$
This identity is what makes the implementation fast. Once factorials and inverse factorials are available modulo \(p\), each basis denominator is obtained in constant time.
A naive interpolation would recompute
$$\prod_{\substack{0\le j\le m\\j\ne i}}(x-j)$$
from scratch for every \(i\), costing \(O(m^2)\). The implementations avoid that by storing prefix and suffix products:
$$P_r=\prod_{j=0}^{r-1}(x-j),\qquad Q_r=\prod_{j=r}^{m}(x-j).$$
Then the numerator for node \(i\) is simply
$$\prod_{\substack{0\le j\le m\\j\ne i}}(x-j)=P_i\,Q_{i+1}.$$
After one forward pass and one backward pass, every interpolation term is available in \(O(1)\), so evaluating \(F_e(x)\) costs only \(O(m)\).
For each prime \(p\), the implementation evaluates the two power sums at
$$x=n\bmod p,$$
which is sufficient because the polynomial is being evaluated in \(\mathbb{F}_p\). It obtains
$$F_k(n)\bmod p\qquad \text{and}\qquad F_{k+1}(n)\bmod p,$$
then applies
$$S_k(n)\equiv (n+1)F_k(n)-F_{k+1}(n)\pmod p.$$
Primes in the interval are detected with deterministic Miller-Rabin for 64-bit integers, preceded by a few small trial divisions. Each prime contributes one residue, and the final answer is the ordinary sum of those residues.
Here we want \(S_2(10)\bmod 17\). First use the identity
$$S_2(10)\equiv 11\,F_2(10)-F_3(10)\pmod{17}.$$
For \(F_2\), the degree is \(3\), so the samples at \(x=0,1,2,3\) are
$$0,\ 1,\ 1+2^2=5,\ 1+2^2+3^2=14.$$
Thus the interpolation data modulo \(17\) is
$$y^{(2)}=(0,1,5,14).$$
For \(F_3\), the degree is \(4\), so the samples at \(x=0,1,2,3,4\) are
$$0,\ 1,\ 9,\ 36,\ 100 \equiv 0,\ 1,\ 9,\ 2,\ 15 \pmod{17}.$$
Hence
$$y^{(3)}=(0,1,9,2,15).$$
Lagrange interpolation at \(x=10\) gives
$$F_2(10)\equiv 11 \pmod{17},\qquad F_3(10)\equiv 16 \pmod{17}.$$
Therefore
$$S_2(10)\equiv 11\cdot 11-16\equiv 3 \pmod{17}.$$
A direct check confirms this: \(S_2(10)=1210\), and \(1210\equiv 3\pmod{17}\).
The C++, Python, and Java implementations all follow the same structure. They scan the interval from \(2\cdot 10^9\) to \(2\cdot 10^9+2000\), keep only primes, and process each surviving modulus independently.
For one prime \(p\), the implementation first builds factorials and inverse factorials up to \(k+2\). Only one modular inverse is needed: after inverting the largest factorial, the remaining inverse factorials are recovered by a backward pass. Because every target prime is far larger than \(k+2\), all factorial values involved are invertible modulo \(p\).
Next, it evaluates the degree-\(k+1\) and degree-\(k+2\) power-sum polynomials. Each sample table is generated by modular exponentiation plus prefix summation. If \(n\bmod p\) already lies among the interpolation nodes, the value is returned immediately from the table; otherwise the prefix/suffix Lagrange formula is used.
Finally, the implementation combines the two polynomial values with
$$S_k(n)\equiv (n+1)F_k(n)-F_{k+1}(n)\pmod p,$$
and adds that residue to the running total. No global modulus is applied to the final accumulation.
Let \(k\) be the exponent parameter and let \(\pi\) be the number of primes in the target interval. For one prime \(p\), building factorial and inverse-factorial tables costs \(O(k)\) modular operations and \(O(k)\) memory. Constructing the sample values for \(F_k\) and \(F_{k+1}\) requires \(O(k)\) modular exponentiation calls; counting modular multiplications explicitly, this is \(O(k\log k)\). The two Lagrange evaluations themselves are linear, so their additional cost is \(O(k)\).
Therefore the overall running time is \(O(\pi k\log k)\) modular multiplications, with \(O(k)\) memory. Since the prime interval is short, the per-prime interpolation work dominates the total runtime.
Für ein Integer \(k\) definieren wir die gewöhnliche Potenzsumme und ihre kumulierte Variante durch
$$F_k(n)=\sum_{j=1}^{n} j^k,\qquad S_k(n)=\sum_{i=1}^{n}F_k(i).$$
Die konkrete Aufgabe ist
$$\sum_{p\in \mathcal{P}} \left(S_{10000}(10^{12}) \bmod p\right),\qquad \mathcal{P}=\{p\text{ prim}: 2\cdot 10^9 \le p \le 2\cdot 10^9+2000\}.$$
Direktes Summieren ist ausgeschlossen: Schon \(n=10^{12}\) ist zu groß, und die verschachtelte Definition von \(S_k(n)\) wäre astronomisch teuer. Deshalb formt die Lösung das Problem in eine Polynom-Auswertung modulo jeder Primzahl um.
Die Implementierungen behandeln jede Primzahl \(p\) im Intervall getrennt. Für jedes \(p\) muss \(S_k(n)\bmod p\) berechnet werden, ohne bis \(n\) zu summieren.
Ausgehend von
$$S_k(n)=\sum_{i=1}^{n}\sum_{j=1}^{i} j^k$$
vertauschen wir die Summationsreihenfolge. Ein fester Wert \(j\) erscheint in allen inneren Summen mit \(i\ge j\), also genau \(n+1-j\)-mal. Damit folgt
$$S_k(n)=\sum_{j=1}^{n}(n+1-j)j^k.$$
Nun trennen wir die beiden Teile:
$$S_k(n)=(n+1)\sum_{j=1}^{n}j^k-\sum_{j=1}^{n}j^{k+1}=(n+1)F_k(n)-F_{k+1}(n).$$
Somit reduziert sich die Aufgabe auf zwei gewöhnliche Potenzsummen modulo \(p\): \(F_k(n)\) und \(F_{k+1}(n)\).
Für festes \(e\) ist
$$F_e(x)=\sum_{t=1}^{x} t^e$$
ein Polynom in \(x\) vom Grad \(e+1\). Das sieht man zum Beispiel an der Vorwärtsdifferenz
$$F_e(x)-F_e(x-1)=x^e,$$
deren rechte Seite Grad \(e\) hat. Also muss \(F_e\) selbst Grad \(e+1\) besitzen.
Ein Polynom dieses Grades ist durch \(e+2\) verschiedene Stützstellen eindeutig bestimmt. Über dem Körper \(\mathbb{F}_p\) verwenden die Implementierungen die Punkte
$$x=0,1,2,\dots,e+1.$$
Das ist zulässig, weil jede Ziel-Primzahl viel größer als \(e+1\) ist und diese Knoten daher modulo \(p\) verschieden bleiben.
Setze \(m=e+1\) und definiere
$$y_i=F_e(i)\bmod p\qquad (0\le i\le m).$$
Diese Werte lassen sich iterativ aufbauen:
$$y_0=0,\qquad y_i=y_{i-1}+i^e \pmod p.$$
Für jede Primzahl \(p\) berechnet die Implementierung also nacheinander \(i^e\bmod p\) für \(i=1,2,\dots,m\), summiert die Werte auf und erhält damit die komplette Interpolationstabelle \((0,y_0),(1,y_1),\dots,(m,y_m)\).
Bernoulli-Zahlen oder explizite Faulhaber-Formeln werden nicht benötigt. Es genügen modulare Potenzierung und Präfixsummen.
Mit den Stützwerten \(y_0,\dots,y_m\) liefert die Lagrange-Formel den Wert für jedes \(x\in\mathbb{F}_p\):
$$F_e(x)=\sum_{i=0}^{m} y_i \prod_{\substack{0\le j\le m\\j\ne i}}\frac{x-j}{i-j}\pmod p.$$
Der Nenner vereinfacht sich zu
$$\prod_{\substack{0\le j\le m\\j\ne i}}(i-j)=i!\,(-1)^{m-i}(m-i)!.$$
Genau diese Identität macht die Auswertung schnell. Sobald Fakultäten und inverse Fakultäten modulo \(p\) vorliegen, ist jeder Basisnenner in konstanter Zeit verfügbar.
Naiv müsste man für jedes \(i\) das Produkt
$$\prod_{\substack{0\le j\le m\\j\ne i}}(x-j)$$
neu berechnen, was \(O(m^2)\) kosten würde. Die Implementierungen vermeiden das mit Präfix- und Suffixprodukten:
$$P_r=\prod_{j=0}^{r-1}(x-j),\qquad Q_r=\prod_{j=r}^{m}(x-j).$$
Dann gilt für den Zähler des Knotens \(i\)
$$\prod_{\substack{0\le j\le m\\j\ne i}}(x-j)=P_i\,Q_{i+1}.$$
Nach einem Vorwärts- und einem Rückwärtsdurchlauf stehen alle Interpolationsterme in \(O(1)\) bereit, also kostet die gesamte Auswertung nur \(O(m)\).
Für jede Primzahl \(p\) werden die beiden Potenzsummen an der Stelle
$$x=n\bmod p$$
ausgewertet. Das genügt, weil das Polynom in \(\mathbb{F}_p\) betrachtet wird. Man erhält also
$$F_k(n)\bmod p\qquad \text{und}\qquad F_{k+1}(n)\bmod p,$$
und setzt dann
$$S_k(n)\equiv (n+1)F_k(n)-F_{k+1}(n)\pmod p.$$
Die Primzahlen im Intervall werden mit deterministischem Miller-Rabin für 64-Bit-Zahlen erkannt, ergänzt durch einige kleine Teilbarkeitstests. Jede Primzahl liefert genau einen Restwert; die Endantwort ist die gewöhnliche Summe dieser Restwerte.
Gesucht ist \(S_2(10)\bmod 17\). Zuerst verwenden wir
$$S_2(10)\equiv 11\,F_2(10)-F_3(10)\pmod{17}.$$
Für \(F_2\) ist der Grad \(3\), also nehmen wir die Werte bei \(x=0,1,2,3\):
$$0,\ 1,\ 1+2^2=5,\ 1+2^2+3^2=14.$$
Damit lautet die Interpolationstabelle modulo \(17\)
$$y^{(2)}=(0,1,5,14).$$
Für \(F_3\) ist der Grad \(4\), also betrachten wir \(x=0,1,2,3,4\):
$$0,\ 1,\ 9,\ 36,\ 100 \equiv 0,\ 1,\ 9,\ 2,\ 15 \pmod{17}.$$
Somit gilt
$$y^{(3)}=(0,1,9,2,15).$$
Die Lagrange-Auswertung bei \(x=10\) ergibt
$$F_2(10)\equiv 11 \pmod{17},\qquad F_3(10)\equiv 16 \pmod{17}.$$
Daher
$$S_2(10)\equiv 11\cdot 11-16\equiv 3 \pmod{17}.$$
Die Direktkontrolle passt: \(S_2(10)=1210\), und \(1210\equiv 3\pmod{17}\).
Die C++-, Python- und Java-Implementierungen folgen derselben Struktur. Sie durchlaufen das Intervall von \(2\cdot 10^9\) bis \(2\cdot 10^9+2000\), verwerfen Nicht-Primzahlen und behandeln jeden verbleibenden Modulus unabhängig.
Für eine Primzahl \(p\) werden zunächst Fakultäten und inverse Fakultäten bis \(k+2\) aufgebaut. Dafür ist nur ein einziges modulares Inverses nötig: Nachdem die größte Fakultät invertiert wurde, entstehen die übrigen inversen Fakultäten in einem Rückwärtslauf. Weil jede Ziel-Primzahl deutlich größer als \(k+2\) ist, sind alle dabei vorkommenden Fakultäten modulo \(p\) invertierbar.
Anschließend wertet die Implementierung die Potenzsummenpolynome vom Grad \(k+1\) und \(k+2\) aus. Jede Stützwerttabelle entsteht durch modulare Potenzierung und Präfixsummen. Falls \(n\bmod p\) bereits unter den Interpolationsknoten liegt, wird der Wert direkt aus der Tabelle gelesen; sonst kommt die Präfix/Suffix-Variante der Lagrange-Formel zum Einsatz.
Zum Schluss werden die beiden Polynomwerte über
$$S_k(n)\equiv (n+1)F_k(n)-F_{k+1}(n)\pmod p$$
kombiniert und zum laufenden Gesamtwert addiert. Eine globale Modulo-Reduktion des Endergebnisses erfolgt nicht.
Sei \(k\) der Exponentenparameter und \(\pi\) die Anzahl der Primzahlen im Zielintervall. Für eine feste Primzahl \(p\) kosten Fakultäten und inverse Fakultäten \(O(k)\) modulare Operationen und \(O(k)\) Speicher. Das Erzeugen der Stützwerte für \(F_k\) und \(F_{k+1}\) benötigt \(O(k)\) Aufrufe modularer Potenzierung; zählt man modulare Multiplikationen explizit, ergibt das \(O(k\log k)\). Die beiden Lagrange-Auswertungen sind linear und tragen also zusätzlich nur \(O(k)\) bei.
Insgesamt ergibt sich damit \(O(\pi k\log k)\) modulare Multiplikationen bei \(O(k)\) Speicher. Weil das Primzahlintervall kurz ist, dominiert die Interpolationsarbeit pro Primzahl die Laufzeit.
Bir \(k\) tamsayısı için sıradan kuvvet toplamını ve onun kümülatif sürümünü
$$F_k(n)=\sum_{j=1}^{n} j^k,\qquad S_k(n)=\sum_{i=1}^{n}F_k(i)$$
şeklinde tanımlayalım. Somut görev şu toplamı hesaplamaktır:
$$\sum_{p\in \mathcal{P}} \left(S_{10000}(10^{12}) \bmod p\right),\qquad \mathcal{P}=\{p\text{ asal}: 2\cdot 10^9 \le p \le 2\cdot 10^9+2000\}.$$
Doğrudan hesaplama mümkün değildir. \(n=10^{12}\) olduğu için tek başına \(F_k(n)\)'yi bile terim terim toplamak imkansızdır; \(S_k(n)\)'nin iç içe toplamı ise çok daha pahalıdır. Bu yüzden çözüm, problemi her asal için modüler polinom değerlendirmesine çevirir.
Uygulamalar, aralıktaki her asal \(p\) için \(S_k(n)\bmod p\) değerini ayrı ayrı hesaplar. Amaç, \(n\)'e kadar toplam almadan bu kalanı elde etmektir.
Tanımdan başlayalım:
$$S_k(n)=\sum_{i=1}^{n}\sum_{j=1}^{i} j^k.$$
Toplam sırasını değiştirirsek, sabit bir \(j\) değeri \(i\ge j\) olan tüm iç toplamların içinde yer alır; yani tam olarak \(n+1-j\) kez sayılır. Böylece
$$S_k(n)=\sum_{j=1}^{n}(n+1-j)j^k$$
olur. İki parçayı ayırınca
$$S_k(n)=(n+1)\sum_{j=1}^{n}j^k-\sum_{j=1}^{n}j^{k+1}=(n+1)F_k(n)-F_{k+1}(n)$$
elde edilir. Yani problem, mod \(p\) altında iki kuvvet toplamını bulmaya iner: \(F_k(n)\) ve \(F_{k+1}(n)\).
Sabit \(e\) için
$$F_e(x)=\sum_{t=1}^{x} t^e$$
ifadesi \(x\)'e göre derece \(e+1\) bir polinomdur. Bunu ileri farktan görebiliriz:
$$F_e(x)-F_e(x-1)=x^e.$$
Sağ tarafın derecesi \(e\) olduğundan, \(F_e\)'in derecesi bir fazlası olan \(e+1\) olur.
Bu nedenle derece \(e+1\) olan polinom, birbirinden farklı \(e+2\) noktadaki değerleriyle tamamen belirlenir. \(\mathbb{F}_p\) sonlu cismi üzerinde uygulama şu düğümleri kullanır:
$$x=0,1,2,\dots,e+1.$$
Hedef asallar \(e+1\)'den çok büyük olduğu için bu noktalar mod \(p\) altında çakışmaz.
\(m=e+1\) olsun ve
$$y_i=F_e(i)\bmod p\qquad (0\le i\le m)$$
tanımını yapalım. Bu değerler iteratif olarak çok kolay üretilir:
$$y_0=0,\qquad y_i=y_{i-1}+i^e \pmod p.$$
Dolayısıyla her asal \(p\) için uygulama, \(i=1,2,\dots,m\) boyunca \(i^e\bmod p\) değerlerini hesaplar, bunları birikimli toplar ve tam enterpolasyon tablosunu elde eder.
Bernoulli sayıları ya da kapalı formüller kullanılmaz. Yalnızca modüler üs alma ve önek toplamlar yeterlidir.
\(y_0,\dots,y_m\) örnekleri bilindiğinde, herhangi bir \(x\in\mathbb{F}_p\) için
$$F_e(x)=\sum_{i=0}^{m} y_i \prod_{\substack{0\le j\le m\\j\ne i}}\frac{x-j}{i-j}\pmod p$$
yazılabilir. Payda kısmı açık biçimde sadeleşir:
$$\prod_{\substack{0\le j\le m\\j\ne i}}(i-j)=i!\,(-1)^{m-i}(m-i)!.$$
Bu özdeşlik kritik önemdedir. Faktöriyel ve ters faktöriyel değerleri mod \(p\) altında hazır olduğunda, her düğümün paydası sabit zamanda elde edilir.
Naif yaklaşımda her \(i\) için
$$\prod_{\substack{0\le j\le m\\j\ne i}}(x-j)$$
ürünü baştan hesaplanır ve maliyet \(O(m^2)\) olur. Uygulamalar bunu önek ve sonek çarpımlarıyla önler:
$$P_r=\prod_{j=0}^{r-1}(x-j),\qquad Q_r=\prod_{j=r}^{m}(x-j).$$
Böylece \(i\). düğümün payı
$$\prod_{\substack{0\le j\le m\\j\ne i}}(x-j)=P_i\,Q_{i+1}$$
şeklinde hemen bulunur. Bir ileri ve bir geri geçişten sonra tüm Lagrange terimleri \(O(1)\) sürede hazırdır; dolayısıyla toplam değerlendirme maliyeti \(O(m)\) olur.
Her asal \(p\) için iki kuvvet toplamı şu noktada değerlendirilir:
$$x=n\bmod p.$$
Bu yeterlidir çünkü polinom \(\mathbb{F}_p\) içinde değerlendirilmektedir. Sonuçta
$$F_k(n)\bmod p\qquad \text{ve}\qquad F_{k+1}(n)\bmod p$$
elde edilir ve ardından
$$S_k(n)\equiv (n+1)F_k(n)-F_{k+1}(n)\pmod p$$
uygulanır. Aralıktaki asallar, 64-bit için deterministik Miller-Rabin ve birkaç küçük bölünebilme kontrolüyle seçilir. Her asal tek bir kalan katkısı verir; nihai cevap bu kalanların normal toplamıdır.
\(S_2(10)\bmod 17\) isteniyor. Önce
$$S_2(10)\equiv 11\,F_2(10)-F_3(10)\pmod{17}$$
özdeşliğini kullanırız. \(F_2\)'nin derecesi \(3\) olduğu için \(x=0,1,2,3\) noktalarındaki değerler
$$0,\ 1,\ 1+2^2=5,\ 1+2^2+3^2=14$$
olur; yani
$$y^{(2)}=(0,1,5,14).$$
\(F_3\)'ün derecesi \(4\) olduğu için \(x=0,1,2,3,4\) noktalarındaki değerler
$$0,\ 1,\ 9,\ 36,\ 100 \equiv 0,\ 1,\ 9,\ 2,\ 15 \pmod{17}$$
ve dolayısıyla
$$y^{(3)}=(0,1,9,2,15)$$
olur. Lagrange enterpolasyonu \(x=10\) için
$$F_2(10)\equiv 11 \pmod{17},\qquad F_3(10)\equiv 16 \pmod{17}$$
sonucunu verir. O halde
$$S_2(10)\equiv 11\cdot 11-16\equiv 3 \pmod{17}.$$
Doğrudan kontrol de aynıdır: \(S_2(10)=1210\) ve \(1210\equiv 3\pmod{17}\).
C++, Python ve Java uygulamalarının yapısı aynıdır. Önce \(2\cdot 10^9\) ile \(2\cdot 10^9+2000\) arasındaki sayılar dolaşılır, asal olmayanlar elenir ve kalan her modül bağımsız olarak işlenir.
Tek bir asal \(p\) için önce \(k+2\)'ye kadar faktöriyel ve ters faktöriyel dizileri hazırlanır. Bunun için tek bir modüler ters yeterlidir: en büyük faktöriyel terslendikten sonra diğer ters faktöriyeller geri yönde elde edilir. Hedef asalların tamamı \(k+2\)'den çok büyük olduğu için kullanılan tüm faktöriyeller mod \(p\) altında terslenebilir.
Daha sonra derece \(k+1\) ve \(k+2\) olan kuvvet-toplam polinomları değerlendirilir. Her örnek tablo modüler üs alma ve önek toplama ile oluşturulur. Eğer \(n\bmod p\) zaten düğümlerden biriyse değer doğrudan tablodan alınır; değilse önek/sonek temelli Lagrange formülü uygulanır.
Son adımda iki polinom değeri
$$S_k(n)\equiv (n+1)F_k(n)-F_{k+1}(n)\pmod p$$
ilişkisiyle birleştirilir ve çalışan toplama eklenir. Son toplam için ayrıca küresel bir mod uygulanmaz.
\(k\) üs parametresi, \(\pi\) ise hedef aralıktaki asal sayısı olsun. Sabit bir asal \(p\) için faktöriyel ve ters faktöriyel tabloları \(O(k)\) modüler işlem ve \(O(k)\) bellek ister. \(F_k\) ile \(F_{k+1}\) için örnek değerleri üretmek \(O(k)\) adet modüler üs alma çağrısı gerektirir; modüler çarpım sayısı üzerinden sayarsak bu kısım \(O(k\log k)\) olur. İki Lagrange değerlendirmesi doğrusal olduğundan ek maliyet \(O(k)\)'dir.
Böylece toplam süre \(O(\pi k\log k)\) modüler çarpım, bellek ise \(O(k)\) olur. Asal aralığı kısa olduğu için çalışma süresini esas olarak her asal için yapılan enterpolasyon belirler.
Para un entero \(k\), definimos la suma de potencias ordinaria y su versión acumulada por
$$F_k(n)=\sum_{j=1}^{n} j^k,\qquad S_k(n)=\sum_{i=1}^{n}F_k(i).$$
La tarea concreta es evaluar
$$\sum_{p\in \mathcal{P}} \left(S_{10000}(10^{12}) \bmod p\right),\qquad \mathcal{P}=\{p\text{ primo}: 2\cdot 10^9 \le p \le 2\cdot 10^9+2000\}.$$
Una evaluación directa no es viable. Ya \(F_k(10^{12})\) sería imposible de sumar término a término, y la definición anidada de \(S_k(n)\) es todavía más costosa. La idea central consiste en convertir el problema, para cada primo \(p\), en una evaluación de polinomios módulo \(p\).
Las implementaciones resuelven el problema primo por primo. Para cada \(p\) del intervalo, hay que obtener \(S_k(n)\bmod p\) sin recorrer todos los enteros hasta \(n\).
Partimos de
$$S_k(n)=\sum_{i=1}^{n}\sum_{j=1}^{i} j^k.$$
Si intercambiamos el orden de suma, un valor fijo \(j\) aparece en todos los sumandos interiores con \(i\ge j\), es decir, exactamente \(n+1-j\) veces. Por tanto,
$$S_k(n)=\sum_{j=1}^{n}(n+1-j)j^k.$$
Separando los dos términos obtenemos
$$S_k(n)=(n+1)\sum_{j=1}^{n}j^k-\sum_{j=1}^{n}j^{k+1}=(n+1)F_k(n)-F_{k+1}(n).$$
Así, todo se reduce a calcular dos sumas de potencias ordinarias módulo \(p\): \(F_k(n)\) y \(F_{k+1}(n)\).
Para un exponente fijo \(e\), la función
$$F_e(x)=\sum_{t=1}^{x} t^e$$
es un polinomio en \(x\) de grado \(e+1\). Una forma de verlo es mediante la diferencia hacia adelante
$$F_e(x)-F_e(x-1)=x^e,$$
cuyo lado derecho tiene grado \(e\); por eso \(F_e\) debe tener grado \(e+1\).
Un polinomio de ese grado queda determinado por \(e+2\) valores en puntos distintos. Sobre \(\mathbb{F}_p\), las implementaciones usan los nodos
$$x=0,1,2,\dots,e+1.$$
Esto funciona porque todos los primos objetivo son muchísimo mayores que \(e+1\), de modo que esos nodos siguen siendo distintos módulo \(p\).
Sea \(m=e+1\) y definamos
$$y_i=F_e(i)\bmod p\qquad (0\le i\le m).$$
Estos valores se generan de forma iterativa:
$$y_0=0,\qquad y_i=y_{i-1}+i^e \pmod p.$$
Por lo tanto, para cada primo \(p\), la implementación calcula \(i^e\bmod p\) para \(i=1,2,\dots,m\), va acumulando la suma y obtiene así toda la tabla de interpolación \((0,y_0),(1,y_1),\dots,(m,y_m)\).
No hacen falta números de Bernoulli ni fórmulas cerradas. Basta con exponenciación modular y sumas prefijo.
Con las muestras \(y_0,\dots,y_m\), la fórmula de Lagrange permite reconstruir el valor en cualquier \(x\in\mathbb{F}_p\):
$$F_e(x)=\sum_{i=0}^{m} y_i \prod_{\substack{0\le j\le m\\j\ne i}}\frac{x-j}{i-j}\pmod p.$$
El denominador se simplifica a
$$\prod_{\substack{0\le j\le m\\j\ne i}}(i-j)=i!\,(-1)^{m-i}(m-i)!.$$
Esa identidad es la clave de eficiencia. Una vez conocidas las factoriales y las inversas factoriales módulo \(p\), cada denominador de base se obtiene en tiempo constante.
Si para cada \(i\) se recalculara desde cero
$$\prod_{\substack{0\le j\le m\\j\ne i}}(x-j),$$
el coste sería \(O(m^2)\). Las implementaciones lo evitan guardando productos prefijo y sufijo:
$$P_r=\prod_{j=0}^{r-1}(x-j),\qquad Q_r=\prod_{j=r}^{m}(x-j).$$
Entonces el numerador del nodo \(i\) es simplemente
$$\prod_{\substack{0\le j\le m\\j\ne i}}(x-j)=P_i\,Q_{i+1}.$$
Tras un recorrido hacia adelante y otro hacia atrás, todos los términos de Lagrange quedan disponibles en \(O(1)\), así que la evaluación completa cuesta solo \(O(m)\).
Para cada primo \(p\), las dos sumas de potencias se evalúan en
$$x=n\bmod p.$$
Eso basta porque el polinomio se está evaluando dentro de \(\mathbb{F}_p\). Se obtienen
$$F_k(n)\bmod p\qquad \text{y}\qquad F_{k+1}(n)\bmod p,$$
y después se aplica
$$S_k(n)\equiv (n+1)F_k(n)-F_{k+1}(n)\pmod p.$$
Los primos del intervalo se detectan con Miller-Rabin determinista para enteros de 64 bits, precedido por unas pocas divisiones pequeñas. Cada primo aporta un residuo y la respuesta final es la suma ordinaria de todos ellos.
Queremos \(S_2(10)\bmod 17\). Primero usamos
$$S_2(10)\equiv 11\,F_2(10)-F_3(10)\pmod{17}.$$
Para \(F_2\), el grado es \(3\), así que los valores en \(x=0,1,2,3\) son
$$0,\ 1,\ 1+2^2=5,\ 1+2^2+3^2=14.$$
La tabla queda
$$y^{(2)}=(0,1,5,14).$$
Para \(F_3\), el grado es \(4\), y en \(x=0,1,2,3,4\) obtenemos
$$0,\ 1,\ 9,\ 36,\ 100 \equiv 0,\ 1,\ 9,\ 2,\ 15 \pmod{17}.$$
Por tanto,
$$y^{(3)}=(0,1,9,2,15).$$
La interpolación de Lagrange en \(x=10\) produce
$$F_2(10)\equiv 11 \pmod{17},\qquad F_3(10)\equiv 16 \pmod{17}.$$
Entonces
$$S_2(10)\equiv 11\cdot 11-16\equiv 3 \pmod{17}.$$
La verificación directa coincide: \(S_2(10)=1210\), y \(1210\equiv 3\pmod{17}\).
Las implementaciones en C++, Python y Java siguen la misma organización. Recorren el intervalo de \(2\cdot 10^9\) a \(2\cdot 10^9+2000\), descartan los no primos y procesan cada módulo primo por separado.
Para un primo \(p\), la implementación construye primero factoriales e inversas factoriales hasta \(k+2\). Solo hace falta una inversión modular: al invertir la factorial más grande, el resto de inversas factoriales se obtiene con un recorrido hacia atrás. Como todos los primos del problema son mucho mayores que \(k+2\), todas esas factoriales son invertibles módulo \(p\).
Después se evalúan los polinomios de sumas de potencias de grado \(k+1\) y \(k+2\). Cada tabla de muestras se forma con exponenciación modular y sumas prefijo. Si \(n\bmod p\) ya coincide con uno de los nodos de interpolación, el valor se devuelve directamente; en caso contrario se usa la versión con prefijos y sufijos de Lagrange.
Por último, los dos valores polinomiales se combinan con
$$S_k(n)\equiv (n+1)F_k(n)-F_{k+1}(n)\pmod p,$$
y ese residuo se añade al acumulador. No se aplica un módulo global a la suma final.
Sea \(k\) el parámetro del exponente y \(\pi\) el número de primos del intervalo objetivo. Para un primo fijo \(p\), construir factoriales e inversas factoriales cuesta \(O(k)\) operaciones modulares y \(O(k)\) memoria. Generar las muestras para \(F_k\) y \(F_{k+1}\) requiere \(O(k)\) llamadas a exponenciación modular; si contamos multiplicaciones modulares de forma explícita, eso es \(O(k\log k)\). Las dos evaluaciones de Lagrange son lineales, así que añaden solo \(O(k)\).
En conjunto, el tiempo total es \(O(\pi k\log k)\) en multiplicaciones modulares, con \(O(k)\) memoria. Como el intervalo de primos es corto, domina el trabajo de interpolación realizado para cada primo.
对一个整数 \(k\)，定义普通幂和与它的累积和：
$$F_k(n)=\sum_{j=1}^{n} j^k,\qquad S_k(n)=\sum_{i=1}^{n}F_k(i).$$
本题实际要求计算
$$\sum_{p\in \mathcal{P}} \left(S_{10000}(10^{12}) \bmod p\right),\qquad \mathcal{P}=\{p\text{ 是素数}: 2\cdot 10^9 \le p \le 2\cdot 10^9+2000\}.$$
直接按定义计算完全不可行。仅仅 \(F_k(10^{12})\) 就无法逐项累加，而 \(S_k(n)\) 还是一个双层求和。真正可行的做法，是把问题化成“对每个素数 \(p\)，在 \(\mathbb{F}_p\) 中求某个多项式在一点的值”。
三份实现的数学核心完全一致：对区间中的每个素数 \(p\)，高效求出 \(S_k(n)\bmod p\)，最后把这些余数作为普通整数相加。
从定义出发：
$$S_k(n)=\sum_{i=1}^{n}\sum_{j=1}^{i} j^k.$$
交换求和顺序。固定某个 \(j\) 后，它会出现在所有满足 \(i\ge j\) 的内层和中，因此总共被计算了 \(n+1-j\) 次。所以
$$S_k(n)=\sum_{j=1}^{n}(n+1-j)j^k.$$
再把它拆开：
$$S_k(n)=(n+1)\sum_{j=1}^{n}j^k-\sum_{j=1}^{n}j^{k+1}=(n+1)F_k(n)-F_{k+1}(n).$$
这一步非常关键，因为它把原题化成了两个更标准的对象：只要能快速求出 \(F_k(n)\bmod p\) 和 \(F_{k+1}(n)\bmod p\)，就能得到 \(S_k(n)\bmod p\)。
对固定指数 \(e\)，定义
$$F_e(x)=\sum_{t=1}^{x} t^e.$$
这是一个关于 \(x\) 的次数 \(e+1\) 多项式。一个直接的理由是前向差分满足
$$F_e(x)-F_e(x-1)=x^e,$$
右边是次数 \(e\) 的多项式，因此左边对应的原函数 \(F_e(x)\) 必然是次数 \(e+1\)。
一旦知道这一点，就可以用多项式插值。次数 \(e+1\) 的多项式，只要知道 \(e+2\) 个不同点上的值，就能唯一确定它。在本题中，最方便的节点就是
$$x=0,1,2,\dots,e+1.$$
因为目标素数都在 \(2\cdot 10^9\) 附近，而 \(e\) 最多只是 \(10001\)，所以这些节点在模 \(p\) 意义下彼此不同，插值完全合法。
设 \(m=e+1\)，并令
$$y_i=F_e(i)\bmod p\qquad (0\le i\le m).$$
这些值不需要任何复杂公式，直接递推即可：
$$y_0=0,\qquad y_i=y_{i-1}+i^e \pmod p.$$
也就是说，只要依次计算 \(1^e,2^e,\dots,m^e\) 在模 \(p\) 下的值，再做前缀累加，就能得到完整的插值数据
$$\bigl(0,y_0\bigr),\bigl(1,y_1\bigr),\dots,\bigl(m,y_m\bigr).$$
这正是实现中使用的方法。它没有依赖伯努利数，也没有调用显式的 Faulhaber 闭式，而是完全从定义出发构造样本。
已知 \(m+1\) 个样本值 \(y_0,\dots,y_m\) 后，可以用拉格朗日公式恢复任意 \(x\in\mathbb{F}_p\) 处的函数值：
$$F_e(x)=\sum_{i=0}^{m} y_i \prod_{\substack{0\le j\le m\\j\ne i}}\frac{x-j}{i-j}\pmod p.$$
其中最重要的简化是分母：
$$\prod_{\substack{0\le j\le m\\j\ne i}}(i-j)=i!\,(-1)^{m-i}(m-i)!.$$
这个公式把看似复杂的分母变成了阶乘与逆阶乘，因此每个基函数的分母都可以 \(O(1)\) 得到。由于题目中的素数远大于 \(m\)，这些阶乘在模 \(p\) 下都可逆，这一点在实现中非常重要。
如果对每个 \(i\) 都单独重算
$$\prod_{\substack{0\le j\le m\\j\ne i}}(x-j),$$
就会产生 \(O(m^2)\) 的代价。实现采用标准优化：预处理前缀积和后缀积
$$P_r=\prod_{j=0}^{r-1}(x-j),\qquad Q_r=\prod_{j=r}^{m}(x-j).$$
这样一来，去掉第 \(i\) 个节点后的分子就能写成
$$\prod_{\substack{0\le j\le m\\j\ne i}}(x-j)=P_i\,Q_{i+1}.$$
也就是说，只需要一次从左到右和一次从右到左的扫描，就能为全部节点准备好分子。于是一次插值求值的总代价从平方时间降为线性时间 \(O(m)\)。
对每个素数 \(p\)，实现并不是在整数意义下把 \(n=10^{12}\) 直接塞进插值公式，而是先取
$$x=n\bmod p.$$
因为整个计算发生在 \(\mathbb{F}_p\) 中，最终只需要这个剩余类即可。先求得
$$F_k(n)\bmod p\qquad \text{和}\qquad F_{k+1}(n)\bmod p,$$
再代回
$$S_k(n)\equiv (n+1)F_k(n)-F_{k+1}(n)\pmod p.$$
区间中的素数通过确定性的 64 位 Miller-Rabin 判定法筛出，前面再配合少量小素数试除。每个素数给出一个余数，最后答案就是这些余数的普通和，而不是再对某个全局模数取模。
我们来完整走一遍小例子，看看整个流程如何落地。目标是计算 \(S_2(10)\bmod 17\)。首先根据步骤 1，
$$S_2(10)\equiv 11\,F_2(10)-F_3(10)\pmod{17}.$$
对于 \(F_2\)，它是次数 \(3\) 的多项式，所以只要 \(x=0,1,2,3\) 四个点的值：
$$F_2(0)=0,\quad F_2(1)=1,\quad F_2(2)=1+2^2=5,\quad F_2(3)=1+2^2+3^2=14.$$
于是有
$$y^{(2)}=(0,1,5,14).$$
对于 \(F_3\)，它是次数 \(4\) 的多项式，所以取 \(x=0,1,2,3,4\)：
$$F_3(0)=0,\quad F_3(1)=1,\quad F_3(2)=9,\quad F_3(3)=36,\quad F_3(4)=100.$$
把它们模 \(17\) 化简后得到
$$y^{(3)}=(0,1,9,2,15).$$
接着用拉格朗日公式在 \(x=10\) 处求值，可得
$$F_2(10)\equiv 11 \pmod{17},\qquad F_3(10)\equiv 16 \pmod{17}.$$
因此
$$S_2(10)\equiv 11\cdot 11-16\equiv 3 \pmod{17}.$$
直接验证也完全一致：按定义算出 \(S_2(10)=1210\)，而 \(1210\equiv 3\pmod{17}\)。这个例子虽然很小，但流程与大数据情形完全相同。
C++、Python 和 Java 三个版本的整体结构一致。它们先遍历 \(2\cdot 10^9\) 到 \(2\cdot 10^9+2000\) 的所有整数，只保留其中的素数，然后对每个素数模数独立处理。
对某个素数 \(p\)，实现首先建立到 \(k+2\) 为止的阶乘表和逆阶乘表。这里不需要为每个数单独求逆，只要先求出最大阶乘的逆元，再通过一次反向递推恢复其余逆阶乘即可。因为所有目标素数都远大于 \(k+2\)，所以这些阶乘值模 \(p\) 后都可逆。
接下来，分别为指数 \(k\) 和 \(k+1\) 构造样本表，并用线性时间的拉格朗日插值求出目标值。如果 \(n\bmod p\) 恰好落在节点 \(0,1,\dots,e+1\) 之中，就直接返回样本值；否则才走前缀积/后缀积的插值路径。
最后，程序用
$$S_k(n)\equiv (n+1)F_k(n)-F_{k+1}(n)\pmod p$$
把两个幂和组合成所需余数，并把这个余数加到总和中。最终累计的是普通整数和，没有额外的全局取模。
记指数参数为 \(k\)，目标区间中的素数个数为 \(\pi\)。对单个素数 \(p\) 而言，构造阶乘和逆阶乘需要 \(O(k)\) 次模运算与 \(O(k)\) 空间。为 \(F_k\) 与 \(F_{k+1}\) 构造样本表时，需要 \(O(k)\) 次模幂调用；如果按模乘次数精确计量，这一部分是 \(O(k\log k)\)。两次拉格朗日求值本身都是线性的，因此附加代价是 \(O(k)\)。
综合起来，总时间复杂度为 \(O(\pi k\log k)\) 次模乘，空间复杂度为 \(O(k)\)。由于素数区间本身很短，实际运行时间的主导部分就是“每个素数上的两次插值求值”。
Для целого \(k\) определим обычную сумму степеней и её накопленную версию:
$$F_k(n)=\sum_{j=1}^{n} j^k,\qquad S_k(n)=\sum_{i=1}^{n}F_k(i).$$
В задаче требуется вычислить
$$\sum_{p\in \mathcal{P}} \left(S_{10000}(10^{12}) \bmod p\right),\qquad \mathcal{P}=\{p\text{ — простое}: 2\cdot 10^9 \le p \le 2\cdot 10^9+2000\}.$$
Прямой подсчёт невозможен: даже \(F_k(10^{12})\) нельзя получить суммированием по определению, а для \(S_k(n)\) пришлось бы выполнять ещё и внешнюю сумму. Поэтому решение переводит задачу в вычисление значений полиномов по модулю каждого простого \(p\).
Все реализации работают одинаково: для каждого простого \(p\) из заданного интервала нужно быстро получить \(S_k(n)\bmod p\), не итерируясь до \(n\).
Начинаем с определения
$$S_k(n)=\sum_{i=1}^{n}\sum_{j=1}^{i} j^k.$$
Поменяем порядок суммирования. Фиксированное значение \(j\) входит во все внутренние суммы с \(i\ge j\), то есть встречается ровно \(n+1-j\) раз. Следовательно,
$$S_k(n)=\sum_{j=1}^{n}(n+1-j)j^k.$$
После разделения членов получаем
$$S_k(n)=(n+1)\sum_{j=1}^{n}j^k-\sum_{j=1}^{n}j^{k+1}=(n+1)F_k(n)-F_{k+1}(n).$$
Итак, задача сводится к вычислению двух сумм степеней по модулю \(p\): \(F_k(n)\) и \(F_{k+1}(n)\).
При фиксированном \(e\) функция
$$F_e(x)=\sum_{t=1}^{x} t^e$$
является полиномом по \(x\) степени \(e+1\). Это видно, например, из разности
$$F_e(x)-F_e(x-1)=x^e,$$
где правая часть имеет степень \(e\). Значит, сама функция должна иметь степень на единицу выше.
Полином степени \(e+1\) однозначно задаётся значениями в \(e+2\) различных точках. В поле \(\mathbb{F}_p\) реализации используют узлы
$$x=0,1,2,\dots,e+1.$$
Это допустимо, потому что все простые из задачи намного больше \(e+1\), а значит эти точки остаются различными по модулю \(p\).
Пусть \(m=e+1\) и
$$y_i=F_e(i)\bmod p\qquad (0\le i\le m).$$
Эти значения строятся итеративно:
$$y_0=0,\qquad y_i=y_{i-1}+i^e \pmod p.$$
То есть для каждого простого \(p\) реализация последовательно вычисляет \(i^e\bmod p\) при \(i=1,2,\dots,m\), накапливает сумму и получает полную таблицу интерполяции \((0,y_0),(1,y_1),\dots,(m,y_m)\).
Никаких чисел Бернулли и готовых формул Фаульхабера здесь не требуется. Используются только быстрое возведение в степень по модулю и префиксные суммы.
Зная значения \(y_0,\dots,y_m\), можно восстановить \(F_e(x)\) для любого \(x\in\mathbb{F}_p\):
$$F_e(x)=\sum_{i=0}^{m} y_i \prod_{\substack{0\le j\le m\\j\ne i}}\frac{x-j}{i-j}\pmod p.$$
Ключевое упрощение относится к знаменателю:
$$\prod_{\substack{0\le j\le m\\j\ne i}}(i-j)=i!\,(-1)^{m-i}(m-i)!.$$
Благодаря этой формуле знаменатели выражаются через факториалы и обратные факториалы, а значит каждый базисный множитель можно получать за \(O(1)\). Поскольку \(p\) намного больше \(m\), все нужные факториалы обратимы по модулю \(p\).
Если для каждого \(i\) отдельно пересчитывать
$$\prod_{\substack{0\le j\le m\\j\ne i}}(x-j),$$
то получится квадратичная сложность \(O(m^2)\). Реализации используют стандартную оптимизацию через префиксные и суффиксные произведения:
$$P_r=\prod_{j=0}^{r-1}(x-j),\qquad Q_r=\prod_{j=r}^{m}(x-j).$$
Тогда числитель для узла \(i\) равен
$$\prod_{\substack{0\le j\le m\\j\ne i}}(x-j)=P_i\,Q_{i+1}.$$
После одного прохода слева направо и одного справа налево все числители становятся доступны за \(O(1)\), и итоговая оценка \(F_e(x)\) стоит только \(O(m)\).
Для каждого простого \(p\) обе суммы степеней вычисляются в точке
$$x=n\bmod p.$$
Этого достаточно, потому что вычисление происходит в поле \(\mathbb{F}_p\). В результате получаются значения
$$F_k(n)\bmod p\qquad \text{и}\qquad F_{k+1}(n)\bmod p,$$
после чего применяется формула
$$S_k(n)\equiv (n+1)F_k(n)-F_{k+1}(n)\pmod p.$$
Простые числа в заданном интервале находятся детерминированным тестом Миллера-Рабина для 64-битных целых, дополненным несколькими малыми проверками делимости. Каждый простой даёт один остаток, а окончательный ответ — это обычная сумма этих остатков.
Найдём \(S_2(10)\bmod 17\). Сначала используем тождество
$$S_2(10)\equiv 11\,F_2(10)-F_3(10)\pmod{17}.$$
Для \(F_2\) степень равна \(3\), поэтому берём значения в точках \(0,1,2,3\):
$$0,\ 1,\ 1+2^2=5,\ 1+2^2+3^2=14.$$
Значит,
$$y^{(2)}=(0,1,5,14).$$
Для \(F_3\) степень равна \(4\), и в точках \(0,1,2,3,4\) получаем
$$0,\ 1,\ 9,\ 36,\ 100 \equiv 0,\ 1,\ 9,\ 2,\ 15 \pmod{17}.$$
Следовательно,
$$y^{(3)}=(0,1,9,2,15).$$
Интерполяция Лагранжа в точке \(x=10\) даёт
$$F_2(10)\equiv 11 \pmod{17},\qquad F_3(10)\equiv 16 \pmod{17}.$$
Отсюда
$$S_2(10)\equiv 11\cdot 11-16\equiv 3 \pmod{17}.$$
Прямая проверка совпадает: \(S_2(10)=1210\), а \(1210\equiv 3\pmod{17}\).
Реализации на C++, Python и Java устроены одинаково. Они просматривают числа от \(2\cdot 10^9\) до \(2\cdot 10^9+2000\), отбрасывают составные и независимо обрабатывают каждый простой модуль.
Для одного простого \(p\) сначала строятся массивы факториалов и обратных факториалов до \(k+2\). При этом нужен только один модульный обратный элемент: после обращения наибольшего факториала остальные обратные факториалы восстанавливаются обратным проходом. Поскольку каждый целевой простой намного больше \(k+2\), все эти значения факториалов обратимы по модулю \(p\).
Далее строятся таблицы значений для степеней \(k\) и \(k+1\), а затем выполняются две линейные по размеру интерполяции. Если \(n\bmod p\) уже совпадает с одним из узлов \(0,1,\dots,e+1\), нужное значение возвращается сразу из таблицы; иначе используется схема с префиксными и суффиксными произведениями.
В конце два найденных значения объединяются по формуле
$$S_k(n)\equiv (n+1)F_k(n)-F_{k+1}(n)\pmod p,$$
и полученный остаток добавляется к общему результату. Глобальное взятие по модулю не выполняется.
Обозначим через \(k\) параметр степени, а через \(\pi\) число простых в целевом интервале. Для одного простого \(p\) построение факториалов и обратных факториалов требует \(O(k)\) модульных операций и \(O(k)\) памяти. Формирование таблиц для \(F_k\) и \(F_{k+1}\) использует \(O(k)\) вызовов быстрого возведения в степень; если считать модульные умножения явно, это даёт \(O(k\log k)\). Обе интерполяции Лагранжа линейны, так что добавляют ещё \(O(k)\).
Итоговая сложность составляет \(O(\pi k\log k)\) модульных умножений при памяти \(O(k)\). Так как интервал простых короткий, доминирует именно работа по интерполяции для каждого простого \(p\).
لعدد صحيح \(k\)، نعرّف مجموع القوى العادي ومجموعه التراكمي كما يلي:
$$F_k(n)=\sum_{j=1}^{n} j^k,\qquad S_k(n)=\sum_{i=1}^{n}F_k(i).$$
والمطلوب فعليًا هو حساب
$$\sum_{p\in \mathcal{P}} \left(S_{10000}(10^{12}) \bmod p\right).$$
وترمز \(\mathcal{P}\) إلى مجموعة الأعداد الأولية الواقعة في المجال \(2\cdot 10^9 \le p \le 2\cdot 10^9+2000\).
الحساب المباشر مستحيل عمليًا. فحتى \(F_k(10^{12})\) لا يمكن جمعه حدًا حدًا، أما \(S_k(n)\) ففيه جمعان متداخلان. لذلك تعتمد الفكرة على تحويل المسألة، لكل عدد أولي \(p\)، إلى تقييم كثير حدود داخل الحقل \(\mathbb{F}_p\).
النسخ الثلاث في C++ وPython وJava تعتمد الفكرة نفسها: نعالج كل عدد أولي \(p\) على حدة، ونحسب \(S_k(n)\bmod p\) من دون المرور على جميع القيم حتى \(n\).
نبدأ من التعريف
$$S_k(n)=\sum_{i=1}^{n}\sum_{j=1}^{i} j^k.$$
إذا بدّلنا ترتيب الجمع، فإن القيمة الثابتة \(j\) تظهر في كل مجموع داخلي يحقق \(i\ge j\)، أي أنها تُحسب بالضبط \(n+1-j\) مرة. ومن ثم
$$S_k(n)=\sum_{j=1}^{n}(n+1-j)j^k.$$
وبفصل الحدود نحصل على
$$S_k(n)=(n+1)\sum_{j=1}^{n}j^k-\sum_{j=1}^{n}j^{k+1}=(n+1)F_k(n)-F_{k+1}(n).$$
إذًا لبّ المسألة هو حساب قيمتين فقط: \(F_k(n)\) و\(F_{k+1}(n)\) بترديد \(p\).
لأس ثابت \(e\)، تكون الدالة
$$F_e(x)=\sum_{t=1}^{x} t^e$$
كثير حدود في \(x\) درجته \(e+1\). ويمكن رؤية ذلك من فرقها الأمامي:
$$F_e(x)-F_e(x-1)=x^e.$$
وبما أن الطرف الأيمن كثير حدود درجته \(e\)، فإن \(F_e(x)\) نفسه يجب أن تكون درجته أعلى بواحد.
ولهذا فإن كثير الحدود من الدرجة \(e+1\) يتحدد تمامًا بمعرفة قيمه في \(e+2\) نقاط مختلفة. وتستخدم الخوارزمية داخل \(\mathbb{F}_p\) العقد التالية:
$$x=0,1,2,\dots,e+1.$$
وهذا جائز لأن جميع الأوليات المستهدفة أكبر بكثير من \(e+1\)، فلا يحدث تطابق بين هذه العقد تحت الترديد \(p\).
لنضع \(m=e+1\)، ولنعرف
$$y_i=F_e(i)\bmod p\qquad (0\le i\le m).$$
هذه القيم تُبنى تراكميًا بسهولة:
$$y_0=0,\qquad y_i=y_{i-1}+i^e \pmod p.$$
أي أن التنفيذ يحسب \(i^e\bmod p\) للقيم \(i=1,2,\dots,m\)، ثم يجمعها تجميعيًا ليحصل على جدول الاستيفاء الكامل
$$\bigl(0,y_0\bigr),\bigl(1,y_1\bigr),\dots,\bigl(m,y_m\bigr).$$
ولا يحتاج هذا إلى أعداد برنولي أو إلى صيغة مغلقة جاهزة؛ فالمعادلات الأساسية هنا تقوم فقط على الأس المعياري والجمع التراكمي.
بعد معرفة \(y_0,\dots,y_m\)، يمكن استرجاع قيمة \(F_e(x)\) لأي \(x\in\mathbb{F}_p\) عبر
$$F_e(x)=\sum_{i=0}^{m} y_i \prod_{\substack{0\le j\le m\\j\ne i}}\frac{x-j}{i-j}\pmod p.$$
أما العامل الأساسي في السرعة فهو تبسيط المقام إلى
$$\prod_{\substack{0\le j\le m\\j\ne i}}(i-j)=i!\,(-1)^{m-i}(m-i)!.$$
وبهذا يتحول المقام إلى جداء من عامليات ومعكوساتها. وبما أن \(p\) أكبر بكثير من \(m\)، فإن هذه القيم كلها قابلة للعكس بترديد \(p\)، وهو ما تستغله الخوارزمية مباشرة.
لو أعدنا حساب
$$\prod_{\substack{0\le j\le m\\j\ne i}}(x-j)$$
من الصفر لكل \(i\)، فستصبح الكلفة \(O(m^2)\). لذلك تستخدم التطبيقات جداءات بادئة ولاحقة:
$$P_r=\prod_{j=0}^{r-1}(x-j),\qquad Q_r=\prod_{j=r}^{m}(x-j).$$
وعندها يصبح بسط العقدة \(i\) مساويًا لـ
$$\prod_{\substack{0\le j\le m\\j\ne i}}(x-j)=P_i\,Q_{i+1}.$$
أي إن مرورًا واحدًا من اليسار إلى اليمين ومرورًا واحدًا من اليمين إلى اليسار يكفيان لتجهيز جميع البسوط. وعندها تصبح عملية الاستيفاء كلها خطية \(O(m)\) بدلًا من تربيعية.
لكل عدد أولي \(p\)، لا تُقيَّم كثيرات الحدود عند \(n\) كعدد صحيح ضخم مباشرة، بل عند
$$x=n\bmod p.$$
وهذا يكفي تمامًا لأن كل الحساب يجري في \(\mathbb{F}_p\). فنحصل أولًا على القيمتين
$$F_k(n)\bmod p,\qquad F_{k+1}(n)\bmod p.$$
ثم نطبق
$$S_k(n)\equiv (n+1)F_k(n)-F_{k+1}(n)\pmod p.$$
وتُستخرج الأوليات داخل المجال بواسطة اختبار Miller-Rabin الحتمي المناسب للأعداد ذات 64 بت، مع بعض اختبارات القسمة الصغيرة قبل ذلك. وكل عدد أولي يضيف باقيًا واحدًا إلى المجموع النهائي، وهذا المجموع هو مجموع صحيح عادي لا يخضع لمعيار عالمي إضافي.
لنحسب \(S_2(10)\bmod 17\). بدايةً نستخدم الهوية
$$S_2(10)\equiv 11\,F_2(10)-F_3(10)\pmod{17}.$$
بالنسبة إلى \(F_2\)، درجته \(3\)، لذا تكفينا القيم عند \(x=0,1,2,3\):
$$0,\ 1,\ 1+2^2=5,\ 1+2^2+3^2=14.$$
ومن ثم
$$y^{(2)}=(0,1,5,14).$$
أما \(F_3\)، فدرجته \(4\)، والقيم عند \(x=0,1,2,3,4\) هي
$$0,\ 1,\ 9,\ 36,\ 100 \equiv 0,\ 1,\ 9,\ 2,\ 15 \pmod{17}.$$
إذًا
$$y^{(3)}=(0,1,9,2,15).$$
وعند تطبيق استيفاء لاغرانج عند \(x=10\) نحصل على
$$F_2(10)\equiv 11 \pmod{17},\qquad F_3(10)\equiv 16 \pmod{17}.$$
وبالتالي
$$S_2(10)\equiv 11\cdot 11-16\equiv 3 \pmod{17}.$$
والتحقق المباشر يؤكد النتيجة: \(S_2(10)=1210\)، و\(1210\equiv 3\pmod{17}\).
تسير تطبيقات C++ وPython وJava جميعها وفق البنية نفسها. فهي تمر على الأعداد من \(2\cdot 10^9\) إلى \(2\cdot 10^9+2000\)، وتحتفظ فقط بالأعداد الأولية، ثم تعالج كل معيار أولي على نحو مستقل.
لكل أولي \(p\)، يبني التنفيذ أولًا جداول العامليات ومعكوسات العامليات حتى \(k+2\). ولا يحتاج إلا إلى معكوس معياري واحد: بعد قلب أكبر عاملي، تُستنتج بقية المعكوسات في مرور عكسي. وبما أن جميع الأوليات المستهدفة أكبر بكثير من \(k+2\)، فإن كل هذه القيم قابلة للعكس بترديد \(p\).
بعد ذلك تُنشأ جداول العينات لدرجتي \(k\) و\(k+1\)، ثم تُجرى عمليتا استيفاء خطيتان. وإذا كان \(n\bmod p\) يقع أصلًا بين العقد \(0,1,\dots,e+1\)، فيُعاد الجواب مباشرة من الجدول؛ وإلا تُستخدم صيغة لاغرانج المعتمدة على الجداءات البادئة واللاحقة.
وفي النهاية يدمج التنفيذ القيمتين وفق
$$S_k(n)\equiv (n+1)F_k(n)-F_{k+1}(n)\pmod p,$$
ثم يضيف الباقي الناتج إلى المجموع الجاري. ولا يوجد معيار نهائي إضافي على المجموع الكلي.
إذا رمزنا إلى وسيط الأس بـ \(k\)، وإلى عدد الأوليات في المجال المستهدف بـ \(\pi\)، فإن بناء جداول العامليات ومعكوساتها لأولية واحدة \(p\) يحتاج إلى \(O(k)\) عملية معيارية وإلى \(O(k)\) ذاكرة. أما بناء العينات لكل من \(F_k\) و\(F_{k+1}\) فيتطلب \(O(k)\) استدعاءً للأس المعياري؛ وإذا قسنا التعقيد بعدد الضربات المعيارية، فهذه المرحلة تكلف \(O(k\log k)\). وتبقى عمليتا استيفاء لاغرانج خطيتين، أي \(O(k)\) إضافية.
إذًا يكون التعقيد الكلي \(O(\pi k\log k)\) من حيث الضربات المعيارية، مع ذاكرة \(O(k)\). ونظرًا إلى أن مجال الأوليات قصير، فإن الجزء المسيطر فعليًا هو عمل الاستيفاء المنفذ لكل أولي.