Problem 421: Prime Factors of $n^{15}+1$

Project Euler Problem 421 Solution

EulerSolve provides an optimized solution for Project Euler Problem 421, Prime Factors of $n^{15}+1$, with C++, Python, Java, and a step-by-step mathematical explanation.

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$....

Detailed mathematical approach

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

Problem page: https://projecteuler.net/problem=421
Multiplicative groups modulo a prime: Wikipedia — Multiplicative group of integers modulo $n$
Primitive roots: Wikipedia — Primitive root modulo $n$
Linear congruences: Wikipedia — Linear congruence theorem
Sieve of Eratosthenes: Wikipedia — Sieve of Eratosthenes

Problem 421 source code

C++

#include <algorithm>
#include <cstdint>
#include <iostream>
#include <limits>
#include <pthread.h>
#include <string>
#include <unistd.h>
#include <vector>

namespace {

using u64 = std::uint64_t;
using u128 = unsigned __int128;
using i64 = std::int64_t;

constexpr u64 kDefaultN = 100'000'000'000ULL;
constexpr int kDefaultM = 100'000'000;

struct Options {
    u64 n = kDefaultN;
    int m = kDefaultM;
    bool run_checkpoints = true;
    bool allow_multithreading = true;
    unsigned requested_threads = 0U;
};

std::string to_string_u128(u128 v) {
    if (v == 0) return "0";
    std::string s;
    while (v > 0) {
        s.push_back(static_cast<char>('0' + static_cast<int>(v % 10)));
        v /= 10;
    }
    std::reverse(s.begin(), s.end());
    return s;
}

bool parse_u64_after_prefix(const std::string& arg, const char* prefix, u64& value) {
    const std::string p(prefix);
    if (arg.rfind(p, 0) != 0) {
        return false;
    }
    const std::string tail = arg.substr(p.size());
    if (tail.empty()) {
        return false;
    }
    u64 parsed = 0ULL;
    for (const char ch : tail) {
        if (ch < '0' || ch > '9') {
            return false;
        }
        const u64 digit = static_cast<u64>(ch - '0');
        if (parsed > (std::numeric_limits<u64>::max() - digit) / 10ULL) {
            return false;
        }
        parsed = parsed * 10ULL + digit;
    }
    value = parsed;
    return true;
}

bool parse_int_after_prefix(const std::string& arg, const char* prefix, int& value) {
    u64 parsed = 0ULL;
    if (!parse_u64_after_prefix(arg, prefix, parsed)) {
        return false;
    }
    if (parsed > static_cast<u64>(std::numeric_limits<int>::max())) {
        return false;
    }
    value = static_cast<int>(parsed);
    return true;
}

bool parse_unsigned_after_prefix(const std::string& arg,
                                 const char* prefix,
                                 unsigned& value) {
    u64 parsed = 0ULL;
    if (!parse_u64_after_prefix(arg, prefix, parsed)) {
        return false;
    }
    if (parsed > static_cast<u64>(std::numeric_limits<unsigned>::max())) {
        return false;
    }
    value = static_cast<unsigned>(parsed);
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (arg == "--single-thread") {
            options.allow_multithreading = false;
            continue;
        }
        if (parse_u64_after_prefix(arg, "--n=", options.n)) {
            continue;
        }
        if (parse_int_after_prefix(arg, "--m=", options.m)) {
            continue;
        }
        if (parse_unsigned_after_prefix(arg, "--threads=", options.requested_threads)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return true;
}

std::vector<int> sieve_primes(const int n) {
    std::vector<bool> sieve(static_cast<std::size_t>(n) + 1, true);
    if (n >= 0) sieve[0] = false;
    if (n >= 1) sieve[1] = false;
    for (int p = 2; static_cast<i64>(p) * p <= n; ++p) {
        if (!sieve[static_cast<std::size_t>(p)]) continue;
        for (int q = p * p; q <= n; q += p) {
            sieve[static_cast<std::size_t>(q)] = false;
        }
    }
    std::vector<int> primes;
    primes.reserve(n / 10);
    for (int i = 2; i <= n; ++i) {
        if (sieve[static_cast<std::size_t>(i)]) primes.push_back(i);
    }
    return primes;
}

u64 gcd_u64(u64 a, u64 b) {
    while (b != 0) {
        u64 t = a % b;
        a = b;
        b = t;
    }
    return a;
}

u64 mod_pow(u64 a, u64 e, u64 mod) {
    u64 r = 1 % mod;
    u64 x = a % mod;
    while (e > 0) {
        if (e & 1ULL) r = (r * x) % mod;
        x = (x * x) % mod;
        e >>= 1ULL;
    }
    return r;
}

u64 primitive_root_of_unity(u64 p, u64 d) {
    if (d == 1) return 1;
    for (u64 a = 2; ; ++a) {
        u64 h = mod_pow(a, (p - 1) / d, p);
        u64 z = 1;
        bool ok = true;
        for (u64 i = 1; i < d; ++i) {
            z = (z * h) % p;
            if (z == 1) {
                ok = false;
                break;
            }
        }
        if (ok) return h;
    }
}

struct SumTask {
    const std::vector<int>* primes = nullptr;
    std::size_t start = 0;
    std::size_t end = 0;
    u64 n = 0;
    u128 sum = 0;
};

static void* sum_worker(void* arg) {
    auto* task = static_cast<SumTask*>(arg);
    const auto& primes = *task->primes;
    const u64 n = task->n;
    u128 local = 0;
    for (std::size_t i = task->start; i < task->end; ++i) {
        const u64 p = static_cast<u64>(primes[i]);
        const u64 d = gcd_u64(15, p - 1);
        const u64 r = primitive_root_of_unity(p, d);
        u64 h = p - 1;
        for (u64 j = 0; j < d; ++j) {
            h = (h * r) % p;
            local += static_cast<u128>((n + p - h) / p) * static_cast<u128>(p);
        }
    }
    task->sum = local;
    return nullptr;
}

unsigned thread_count(bool allow_multithreading, unsigned requested_threads) {
    if (!allow_multithreading) {
        return 1U;
    }
    unsigned threads = requested_threads;
    if (threads == 0U) {
        long n = sysconf(_SC_NPROCESSORS_ONLN);
        if (n < 1) n = 1;
        threads = static_cast<unsigned>(n);
    }
    return std::max(1U, threads);
}

u128 sum_for(u64 n, int m, bool allow_multithreading, unsigned requested_threads) {
    const auto primes = sieve_primes(m);
    if (primes.empty()) {
        return 0;
    }

    const unsigned threads = thread_count(allow_multithreading, requested_threads);
    if (threads < 2U || primes.size() < 20000U) {
        u128 total = 0;
        for (int p_int : primes) {
            const u64 p = static_cast<u64>(p_int);
            const u64 d = gcd_u64(15, p - 1);
            const u64 r = primitive_root_of_unity(p, d);
            u64 h = p - 1;
            for (u64 j = 0; j < d; ++j) {
                h = (h * r) % p;
                total += static_cast<u128>((n + p - h) / p) * static_cast<u128>(p);
            }
        }
        return total;
    }

    const std::size_t count = primes.size();
    const unsigned active = std::min<unsigned>(threads, static_cast<unsigned>(count));
    const std::size_t block = (count + active - 1) / active;
    std::vector<SumTask> tasks(active);
    std::vector<pthread_t> workers(active);

    for (unsigned t = 0; t < active; ++t) {
        const std::size_t start = static_cast<std::size_t>(t) * block;
        const std::size_t end = std::min(count, start + block);
        tasks[t] = {&primes, start, end, n, 0};
        pthread_create(&workers[t], nullptr, sum_worker, &tasks[t]);
    }

    u128 total = 0;
    for (unsigned t = 0; t < active; ++t) {
        pthread_join(workers[t], nullptr);
        total += tasks[t].sum;
    }
    return total;
}

u64 s_single_bruteforce(u64 n, int m) {
    const auto primes = sieve_primes(m);
    u64 sum = 0;
    for (int p_int : primes) {
        const u64 p = static_cast<u64>(p_int);
        const u64 r = mod_pow(n % p, 15ULL, p);
        if ((r + 1ULL) % p == 0ULL) sum += p;
    }
    return sum;
}

u128 brute_small(u64 n_limit, int m) {
    u128 s = 0;
    for (u64 n = 1; n <= n_limit; ++n) {
        s += s_single_bruteforce(n, m);
    }
    return s;
}

bool run_checkpoints() {
    if (s_single_bruteforce(2ULL, 10) != 3ULL) return false;
    if (s_single_bruteforce(2ULL, 1000) != 345ULL) return false;
    if (s_single_bruteforce(10ULL, 100) != 31ULL) return false;
    if (s_single_bruteforce(10ULL, 1000) != 483ULL) return false;
    if (sum_for(200ULL, 1000, false, 1) != brute_small(200ULL, 1000)) return false;
    if (sum_for(1000ULL, 5000, false, 1) != brute_small(1000ULL, 5000)) return false;
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }

    if (options.run_checkpoints && !run_checkpoints()) {
        std::cerr << "Checkpoint failed\n";
        return 2;
    }

    const u128 ans = sum_for(options.n, options.m, options.allow_multithreading, options.requested_threads);
    std::cout << to_string_u128(ans) << '\n';
    return 0;
}

Python

import math

def solve():
    N = 100_000_000_000
    M = 100_000_000

    def sieve_primes(n):
        sieve = bytearray([1]) * (n + 1)
        sieve[0] = sieve[1] = 0
        for p in range(2, int(n**0.5) + 1):
            if sieve[p]:
                for q in range(p*p, n+1, p):
                    sieve[q] = 0
        return [i for i in range(2, n+1) if sieve[i]]

    def mod_pow(a, e, mod):
        r = 1 % mod
        x = a % mod
        while e > 0:
            if e & 1: r = r * x % mod
            x = x * x % mod
            e >>= 1
        return r

    def primitive_root_of_unity(p, d):
        if d == 1: return 1
        for a in range(2, p):
            h = mod_pow(a, (p-1)//d, p)
            z = 1
            ok = True
            for _ in range(1, d):
                z = z * h % p
                if z == 1:
                    ok = False
                    break
            if ok: return h

    primes = sieve_primes(M)
    total = 0
    for p in primes:
        d = math.gcd(15, p - 1)
        r = primitive_root_of_unity(p, d)
        h = p - 1
        for _ in range(d):
            h = h * r % p
            total += ((N + p - h) // p) * p

    return str(total)

if __name__ == '__main__':
    print(solve())

Java

public class Euler421 {
    public static String solve() {
        long N = 100_000_000_000L;
        int M = 100_000_000;
        boolean[] sieve = new boolean[M + 1];
        java.util.Arrays.fill(sieve, true);
        sieve[0] = sieve[1] = false;
        for (int p = 2; (long) p * p <= M; p++)
            if (sieve[p])
                for (int q = p * p; q <= M; q += p)
                    sieve[q] = false;
        long total = 0;
        for (int p = 2; p <= M; p++) {
            if (!sieve[p])
                continue;
            int d = gcd(15, p - 1);
            long r = primRootOfUnity(p, d);
            long h = p - 1;
            for (int i = 0; i < d; i++) {
                h = h * r % p;
                total += ((N + p - h) / p) * p;
            }
        }
        return String.valueOf(total);
    }

    static int gcd(int a, int b) {
        while (b != 0) {
            int t = b;
            b = a % b;
            a = t;
        }
        return a;
    }

    static long primRootOfUnity(int p, int d) {
        if (d == 1)
            return 1;
        for (int a = 2; a < p; a++) {
            long h = modPow(a, (p - 1) / d, p);
            long z = 1;
            boolean ok = true;
            for (int i = 1; i < d; i++) {
                z = z * h % p;
                if (z == 1) {
                    ok = false;
                    break;
                }
            }
            if (ok)
                return h;
        }
        return 1;
    }

    static long modPow(long b, long e, long m) {
        long r = 1;
        b %= m;
        while (e > 0) {
            if ((e & 1) != 0)
                r = r * b % m;
            b = b * b % m;
            e >>= 1;
        }
        return r;
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}