Detailed mathematical approach

Problem Summary

The goal is to evaluate

$$S(n)=\sum_{k=1}^{n}\sum_{p=1}^{n}(1-k^2)^p \pmod{M},\qquad M=10^9+7.$$

For the actual input \(n=10^6\), a direct double loop would require about \(10^{12}\) term updates, so the computation must be reorganized. The key observation is that for each fixed \(k\), the inner sum is a finite geometric series. Once that series is collapsed to a closed form, the problem becomes a single loop over \(k\), with every contribution evaluated by modular exponentiation.

A useful checkpoint is

$$S(4)=51160,$$

which the implementation uses to confirm that the closed form matches direct evaluation on a small case.

Mathematical Approach

The entire method comes from simplifying the inner sum for one fixed value of \(k\).

Step 1: Isolate the Inner Sum for One \(k\)

Define

$$u_k=1-k^2.$$

Then the contribution of this \(k\) to the outer sum is

$$T_k=\sum_{p=1}^{n}u_k^p.$$

Therefore

$$S(n)=\sum_{k=1}^{n}T_k.$$

This separates the problem into \(n\) independent geometric sums.

Step 2: Collapse the Geometric Series

For \(u_k\neq 1\), the usual finite geometric-series identity gives

$$T_k=\frac{u_k^{n+1}-u_k}{u_k-1}=u_k\cdot\frac{u_k^n-1}{u_k-1}.$$

This is exactly the form used by the implementation: one fast power computes \(u_k^n\), and the remaining work is a small number of modular multiplications.

Step 3: Simplify the Denominator

Because \(u_k=1-k^2\), we have

$$u_k-1=-k^2.$$

So the same quantity can also be written as

$$T_k=\frac{u_k-u_k^{n+1}}{k^2}.$$

This makes the arithmetic easier to interpret: the denominator is just \(k^2\). Under the Euler modulus \(M=10^9+7\), every \(k\) in the range \(1\le k\le 10^6\) is nonzero modulo \(M\), so the denominator is invertible for the actual problem input.

Step 4: Replace Division by a Modular Inverse

All calculations are performed modulo the prime \(M\). For any nonzero residue \(a\), Fermat's little theorem gives

$$a^{-1}\equiv a^{M-2}\pmod{M}.$$

Hence, for \(u_k\not\equiv 1\pmod{M}\),

$$T_k\equiv u_k\,(u_k^n-1)\,(u_k-1)^{-1}\pmod{M}.$$

This turns the whole computation into repeated modular exponentiation and multiplication, which is far cheaper than summing \(n\) powers for each \(k\).

Step 5: Handle the Degenerate Case Safely

If \(u_k\equiv 1\pmod{M}\), then every term in the inner sum equals \(1\), so

$$T_k\equiv n\pmod{M}.$$

For the actual Project Euler parameters this case does not occur, because \(u_k\equiv 1\) would imply \(k^2\equiv 0\pmod{M}\), impossible for \(1\le k<M\). Still, the implementation keeps this branch so the formula never attempts to divide by zero modulo \(M\).

Worked Example: \(n=4\)

Now compute the checkpoint directly.

For \(k=1\), \(u_1=1-1^2=0\), so

$$T_1=0+0+0+0=0.$$

For \(k=2\), \(u_2=-3\), hence

$$T_2=(-3)+(-3)^2+(-3)^3+(-3)^4=-3+9-27+81=60.$$

For \(k=3\), \(u_3=-8\), so

$$T_3=-8+64-512+4096=3640.$$

For \(k=4\), \(u_4=-15\), therefore

$$T_4=-15+225-3375+50625=47460.$$

Adding the four contributions gives

$$S(4)=0+60+3640+47460=51160,$$

which matches the checkpoint used by the implementation.

How the Code Works

The C++, Python, and Java implementations all follow the same mathematical structure. They iterate once over \(k=1,2,\dots,n\), reduce \(1-k^2\) into the range \(0\) to \(M-1\), and then evaluate the closed form for \(T_k\) modulo \(M\).

For each \(k\), the implementation computes the power \(u_k^n \bmod M\) with binary exponentiation. If the denominator is nonzero, it computes the modular inverse by raising that denominator to the power \(M-2\), again with binary exponentiation. The resulting contribution is added to the running total modulo \(M\).

One implementation also includes two small validation steps: the known checkpoint \(S(4)=51160\), and a comparison between the optimized formula and direct summation for a small input. Those checks confirm that the algebraic reduction is correct before the large target value is evaluated.

Complexity Analysis

There is one outer loop over \(k=1\) to \(n\). For each \(k\), the implementation performs a constant amount of arithmetic plus modular exponentiation for \(u_k^n\), and modular exponentiation for the inverse. Binary exponentiation costs \(O(\log n)\) for the first power and \(O(\log M)\) for the inverse, so the full running time is

$$O\bigl(n(\log n+\log M)\bigr).$$

With the fixed Euler modulus \(M=10^9+7\), this is effectively \(O(n\log n)\). The memory usage is

$$O(1),$$

since the program stores only a handful of current residues. This is a major improvement over the naive \(O(n^2)\) direct summation.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=479
2. Geometric series: Wikipedia — Geometric series
3. Fermat's little theorem: Wikipedia — Fermat's little theorem
4. Modular exponentiation: Wikipedia — Modular exponentiation