Detailed mathematical approach

Problem Summary

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.

Mathematical Approach

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\).

Step 1: Collapse the Double Sum

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)\).

Step 2: View \(F_e(x)\) as a Polynomial

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\).

Step 3: Build the Sample Values

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.

Step 4: Use Lagrange Interpolation over \(\mathbb{F}_p\)

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.

Step 5: Evaluate All Numerators in Linear 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)\).

Step 6: Assemble \(S_k(n)\bmod p\) and Sum over Primes

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.

Worked Example: \(k=2\), \(n=10\), \(p=17\)

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}\).

How the Code Works

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.

Complexity Analysis

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.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=487
2. Power sums and Faulhaber-type identities: Wikipedia — Faulhaber's formula
3. Polynomial interpolation: Wikipedia — Lagrange polynomial
4. Primality testing: Wikipedia — Miller-Rabin primality test