Problem Summary

Define

$$f(x,M)=\sum_{\substack{p\le M\\ p\mid x^{15}+1}} p,$$

and let

$$S(N,M)=\sum_{x=1}^{N} f(x,M).$$

We must evaluate this exactly for very large parameters, so checking every \(x\) separately is impossible. The key observation is that divisibility by a fixed prime depends only on residue classes modulo that prime.

Mathematical Approach

1. Reverse the order of summation

Instead of iterating over \(x\), sum prime by prime:

$$S(N,M)=\sum_{p\le M} p\,A_p(N),$$

where

$$A_p(N)=\left|\left\{1\le x\le N : x^{15}\equiv -1 \pmod p\right\}\right|.$$

So for each prime \(p\), the entire problem reduces to counting how many residue classes solve \(x^{15}\equiv -1 \pmod p\), and then counting how often those classes occur up to \(N\).

2. Count the solutions of \(x^{15}\equiv -1 \pmod p\)

For odd primes, the multiplicative group \((\mathbb{Z}/p\mathbb{Z})^\times\) is cyclic of order \(p-1\). If \(g\) is a generator and \(x=g^k\), then \(-1=g^{(p-1)/2}\), so the congruence becomes

$$15k\equiv \frac{p-1}{2}\pmod{p-1}.$$

A linear congruence \(ak\equiv b\pmod m\) has solutions exactly when \(\gcd(a,m)\mid b\), and in that case it has \(\gcd(a,m)\) distinct residue classes modulo \(m\). Therefore the number of solutions is

$$d_p=\gcd(15,p-1).$$

This always works for odd \(p\): the number \(d_p\) is an odd divisor of \(15\), hence an odd divisor of \(p-1\), and every odd divisor of \(p-1\) also divides \((p-1)/2\). The special case \(p=2\) also fits the same formula, because \(d_2=\gcd(15,1)=1\) and the unique solution is \(x\equiv 1\pmod 2\).

3. Describe all solution classes explicitly

Because \((-1)^{15}=-1\), one solution is already known: \(x\equiv -1\pmod p\). Every other solution differs from it by a 15th root of unity. The subgroup

$$U_p=\left\{u\in (\mathbb{Z}/p\mathbb{Z})^\times : u^{15}\equiv 1\pmod p\right\}$$

has exactly \(d_p\) elements. If \(r\) has exact order \(d_p\), then

$$U_p=\{r^0,r^1,\dots,r^{d_p-1}\},$$

so the full solution set is

$$h_j\equiv -r^j\pmod p,\qquad 0\le j<d_p.$$

This is exactly what the implementation uses: find one element of order \(d_p\), then walk once around the resulting geometric cycle. Since \(d_p\in\{1,3,5,15\}\), each prime produces at most 15 relevant residue classes.

4. Count integers up to \(N\) in one residue class

If \(1\le h\le p-1\), the integers \(x\le N\) with \(x\equiv h\pmod p\) are

$$h,\ h+p,\ h+2p,\ \dots,$$

so their count is

$$\left\lfloor\frac{N-h}{p}\right\rfloor+1=\left\lfloor\frac{N+p-h}{p}\right\rfloor.$$

This formula automatically becomes \(0\) when the first representative \(h\) is already larger than \(N\). Therefore the contribution of one prime is

$$C_p(N)=p\sum_{j=0}^{d_p-1}\left\lfloor\frac{N+p-h_j}{p}\right\rfloor,$$

and the final answer is

$$\boxed{S(N,M)=\sum_{p\le M} p\sum_{j=0}^{d_p-1}\left\lfloor\frac{N+p-h_j}{p}\right\rfloor.}$$

5. Worked example: \(N=10,\ M=20\)

For \(p=2\), the only class is \(h=1\), so there are \(5\) odd numbers up to \(10\), giving \(C_2=2\cdot 5=10\).

For \(p=7\), we have \(d_7=\gcd(15,6)=3\). One element of order \(3\) is \(r=4\), hence the solution classes are \(h\in\{3,5,6\}\). Up to \(10\), these occur \(2,1,1\) times respectively, so

$$C_7=7(2+1+1)=28.$$

For \(p=11\), \(d_{11}=5\), and the five classes are \(h\in\{7,6,2,8,10\}\). Each appears exactly once in \(1,\dots,10\), so

$$C_{11}=11\cdot 5=55.$$

The remaining contributing primes up to \(20\) are \(3,5,13,19\), with contributions \(9,10,26,19\). Prime \(17\) contributes \(0\), because its only solution class is \(16\pmod{17}\), which does not appear below \(10\). Thus

$$S(10,20)=10+9+10+28+55+26+19=157.$$

How the Code Works

The C++, Python, and Java implementations all follow the same number-theoretic plan. First they sieve all primes up to \(M\). For each prime, they compute \(d=\gcd(15,p-1)\), construct one element of exact order \(d\), generate the \(d\) solution classes for \(x^{15}\equiv -1\pmod p\), and add \(p\) times the number of integers up to \(N\) in each class. No search over all \(x\le N\) ever occurs.

The C++ implementation also exploits the fact that prime contributions are independent: once the prime list has been built, different ranges of primes can be processed in parallel and their partial sums added at the end.

Complexity Analysis

Generating all primes up to \(M\) with a sieve costs \(O(M\log\log M)\) time and \(O(M)\) memory in the direct form used here. After that, the work per prime is small: one gcd, a short search for an element of order \(d\), and then processing only \(d\le 15\) residue classes. So the runtime is dominated by the prime scan, not by \(N\), which is the decisive improvement over brute force.

References

1. Problem page: https://projecteuler.net/problem=421
2. Multiplicative groups modulo a prime: Wikipedia — Multiplicative group of integers modulo \(n\)
3. Primitive roots: Wikipedia — Primitive root modulo \(n\)
4. Linear congruences: Wikipedia — Linear congruence theorem
5. Sieve of Eratosthenes: Wikipedia — Sieve of Eratosthenes