Problem 659: Largest Prime

Project Euler Problem 659 Solution

EulerSolve provides an optimized solution for Project Euler Problem 659, Largest Prime, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For each positive integer $k$, define $$N_k=4k^2+1,$$ and let $P(k)$ denote the largest prime factor of $N_k$. The goal is to compute $$S(K)=\sum_{k=1}^{K} P(k)\pmod{10^{18}},\qquad K=10^7.$$ Factoring each $N_k$ independently would be far too slow, so the efficient approach is to sieve prime divisors across all $k$ at once. Mathematical Approach The key observation is that a prime divisor of $4k^2+1$ forces $k$ into very specific residue classes. That turns the problem into a structured sieve on arithmetic progressions rather than $10^7$ separate factorizations. Step 1: Restrict the Shape of the Prime Divisors If a prime $p$ divides $4k^2+1$, then $$4k^2\equiv -1\pmod p,$$ so $-1$ must be a quadratic residue modulo $p$. Since $4k^2+1$ is always odd, $p=2$ never occurs. For an odd prime, the existence of a solution to $$x^2\equiv -1\pmod p$$ implies $$p\equiv 1\pmod 4.$$ Therefore only primes congruent to $1 \bmod 4$ can appear as prime divisors of the numbers in this family. Step 2: Convert a Prime Divisor into Residue Classes of $k$ Fix an odd prime $p\equiv 1\pmod 4$....

Detailed mathematical approach

Problem Summary

For each positive integer $k$, define

$$N_k=4k^2+1,$$

and let $P(k)$ denote the largest prime factor of $N_k$. The goal is to compute

$$S(K)=\sum_{k=1}^{K} P(k)\pmod{10^{18}},\qquad K=10^7.$$

Factoring each $N_k$ independently would be far too slow, so the efficient approach is to sieve prime divisors across all $k$ at once.

Mathematical Approach

The key observation is that a prime divisor of $4k^2+1$ forces $k$ into very specific residue classes. That turns the problem into a structured sieve on arithmetic progressions rather than $10^7$ separate factorizations.

Step 1: Restrict the Shape of the Prime Divisors

If a prime $p$ divides $4k^2+1$, then

$$4k^2\equiv -1\pmod p,$$

so $-1$ must be a quadratic residue modulo $p$. Since $4k^2+1$ is always odd, $p=2$ never occurs. For an odd prime, the existence of a solution to

$$x^2\equiv -1\pmod p$$

implies

$$p\equiv 1\pmod 4.$$

Therefore only primes congruent to $1 \bmod 4$ can appear as prime divisors of the numbers in this family.

Step 2: Convert a Prime Divisor into Residue Classes of $k$

Fix an odd prime $p\equiv 1\pmod 4$. Once we find a square root $x$ of $-1$ modulo $p$, we have

$$x^2\equiv -1\pmod p.$$

Because $4k^2=(2k)^2$, the divisibility condition $p\mid 4k^2+1$ is equivalent to

$$2k\equiv \pm x \pmod p.$$

Since $2$ is invertible modulo every odd prime, this becomes

$$k\equiv \pm x\cdot 2^{-1}\pmod p.$$

So each eligible prime contributes at most two arithmetic progressions in $k$. The implementations use Tonelli-Shanks to obtain the modular square root and then sweep exactly those two progressions.

Step 3: Remove the Full Prime Power

When a progression reaches an index $k$, the current value attached to that index is checked for divisibility by $p$. If $p$ divides it, the algorithm divides by $p$ repeatedly until the factor is gone:

$$N_k=p^{v_p(N_k)}\cdot M_k,\qquad p\nmid M_k.$$

That repeated division matters because values such as $325=5^2\cdot 13$ contain higher powers. Processing the complete $p$-adic valuation ensures the remaining cofactor is always correct after the sweep for $p$.

Step 4: Track the Largest Processed Prime Factor

The prime sweep runs in increasing order. For each $k$, the implementation stores two pieces of information:

$$\text{current remainder of }N_k,\qquad \text{largest processed prime that has divided }N_k.$$

Because the primes are handled from small to large, whenever a new prime divides the current remainder it automatically becomes the largest processed prime seen so far for that $k$. After all eligible small primes have been removed, the largest prime factor is simply the larger of those two stored quantities.

Step 5: Why Scanning Primes Only up to $2K$ Is Enough

The implementations generate primes only up to $2K$. This is sufficient because every prime factor $p\le 2K$ of any $N_k$ is removed during the sieve. After that, the remaining cofactor has no prime divisor $\le 2K$.

Suppose the leftover cofactor were composite and greater than $1$. Then it would have at least two prime factors, both strictly larger than $2K$, so it would be at least

$$(2K+1)^2=4K^2+4K+1,$$

which is larger than every value in the family because

$$4k^2+1\le 4K^2+1.$$

This contradiction shows that the leftover is either $1$ or one prime. Hence the final answer for each $k$ is

$$P(k)=\max(\text{largest processed prime for }k,\ \text{leftover cofactor for }k).$$

Worked Example: The First Ten Terms

For $K=10$, the values are small enough to inspect directly:

$$\begin{aligned} 4(1)^2+1&=5 &&\Rightarrow P(1)=5,\\ 4(2)^2+1&=17 &&\Rightarrow P(2)=17,\\ 4(3)^2+1&=37 &&\Rightarrow P(3)=37,\\ 4(4)^2+1&=65=5\cdot 13 &&\Rightarrow P(4)=13,\\ 4(5)^2+1&=101 &&\Rightarrow P(5)=101,\\ 4(6)^2+1&=145=5\cdot 29 &&\Rightarrow P(6)=29,\\ 4(7)^2+1&=197 &&\Rightarrow P(7)=197,\\ 4(8)^2+1&=257 &&\Rightarrow P(8)=257,\\ 4(9)^2+1&=325=5^2\cdot 13 &&\Rightarrow P(9)=13,\\ 4(10)^2+1&=401 &&\Rightarrow P(10)=401. \end{aligned}$$

Therefore

$$S(10)=5+17+37+13+101+29+197+257+13+401=1070.$$

This is also a useful checkpoint for the implementation. The sieve interpretation is visible here: the prime $5$ hits the residue classes $k\equiv 1,4\pmod 5$, so among the first ten indices it strips factors from $k=1,4,6,9$.

How the Code Works

The C++, Python, and Java implementations all follow the same pipeline. First they build all primes up to $2K$ with a linear sieve. Then they initialize an array containing the values $4k^2+1$ and a second array that stores the largest prime factor already confirmed during the sweep.

For each prime $p\equiv 1\pmod 4$, the implementation computes a square root of $-1$ modulo $p$, converts that root into the two valid residue classes of $k$, and walks those arithmetic progressions with step $p$. At each visited index, if the current remainder is divisible by $p$, the code divides out all powers of $p$ and updates the stored largest processed prime.

After every eligible prime up to $2K$ has been handled, each index $k$ has a remaining cofactor that is either $1$ or prime. The implementation takes the maximum of that leftover value and the largest processed prime, adds it to the running total, and reduces the sum modulo $10^{18}$.

Complexity Analysis

Generating all primes up to $2K$ with a linear sieve costs $O(K)$ time and $O(K)$ memory. The residue-class sweeps contribute about

$$\sum_{\substack{p\le 2K\\ p\equiv 1\!\!\!\pmod 4}} \frac{2K}{p},$$

which has the usual near-harmonic sieve growth, so the overall running time is $O(K\log\log K)$ in practice. The repeated divisions by prime powers do not change that asymptotic bound, and the dominant storage remains the two arrays of length $K+1$, so the memory usage is $O(K)$.

Footnotes and References

Problem page: https://projecteuler.net/problem=659
Tonelli-Shanks algorithm: Wikipedia — Tonelli-Shanks algorithm
Quadratic residue: Wikipedia — Quadratic residue
Euler's criterion: Wikipedia — Euler's criterion
Modular square root: Wikipedia — Modular square root

Problem 659 source code

C++

#include <algorithm>
#include <cassert>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <vector>

namespace {

using u32 = std::uint32_t;
using u64 = std::uint64_t;
using u128 = __uint128_t;

constexpr u64 kMod18 = 1'000'000'000'000'000'000ULL;

u64 mod_mul(u64 a, u64 b, u64 mod) {
    return static_cast<u64>((static_cast<u128>(a) * static_cast<u128>(b)) % mod);
}

u64 mod_pow(u64 base, u64 exp, u64 mod) {
    u64 result = 1 % mod;
    base %= mod;
    while (exp > 0) {
        if (exp & 1ULL) result = mod_mul(result, base, mod);
        base = mod_mul(base, base, mod);
        exp >>= 1ULL;
    }
    return result;
}

u64 tonelli_shanks(u64 n, u64 p) {
    if (n == 0) return 0;
    if (p == 2) return n;
    if (mod_pow(n, (p - 1) / 2, p) != 1) return 0;
    if (p % 4 == 3) return mod_pow(n, (p + 1) / 4, p);

    u64 q = p - 1;
    int s = 0;
    while ((q & 1ULL) == 0) {
        q >>= 1ULL;
        ++s;
    }

    u64 z = 2;
    while (mod_pow(z, (p - 1) / 2, p) != p - 1) ++z;

    u64 m = static_cast<u64>(s);
    u64 c = mod_pow(z, q, p);
    u64 t = mod_pow(n, q, p);
    u64 r = mod_pow(n, (q + 1) / 2, p);

    while (t != 1) {
        u64 tt = t;
        u64 i = 0;
        while (tt != 1) {
            tt = mod_mul(tt, tt, p);
            ++i;
        }
        u64 b = mod_pow(c, 1ULL << (m - i - 1), p);
        r = mod_mul(r, b, p);
        c = mod_mul(b, b, p);
        t = mod_mul(t, c, p);
        m = i;
    }

    return r;
}

std::vector<int> sieve_primes(int n) {
    std::vector<int> primes;
    std::vector<bool> composite(static_cast<std::size_t>(n + 1), false);
    for (int i = 2; i <= n; ++i) {
        if (!composite[static_cast<std::size_t>(i)]) primes.push_back(i);
        for (int p : primes) {
            const long long v = 1LL * p * i;
            if (v > n) break;
            composite[static_cast<std::size_t>(v)] = true;
            if (i % p == 0) break;
        }
    }
    return primes;
}

u64 solve_case(int kmax) {
    const int limit = 2 * kmax;
    const std::vector<int> primes = sieve_primes(limit);

    std::vector<u64> rem(static_cast<std::size_t>(kmax + 1), 0);
    std::vector<u32> last_small(static_cast<std::size_t>(kmax + 1), 1);
    for (int k = 1; k <= kmax; ++k) {
        rem[static_cast<std::size_t>(k)] = 4ULL * static_cast<u64>(k) * static_cast<u64>(k) + 1ULL;
    }

    for (int p : primes) {
        if (p == 2 || (p & 3) != 1) continue;

        const u64 root = tonelli_shanks(static_cast<u64>(p - 1), static_cast<u64>(p));
        if (root == 0) continue;

        const u64 inv2 = static_cast<u64>(p + 1) / 2ULL;
        const int k1 = static_cast<int>(mod_mul(root, inv2, static_cast<u64>(p)));
        const int k2 = (k1 == 0) ? 0 : (p - k1);

        auto process_residue = [&](int residue) {
            int start = residue;
            if (start <= 0) start += p;
            for (int k = start; k <= kmax; k += p) {
                u64& v = rem[static_cast<std::size_t>(k)];
                if (v % static_cast<u64>(p) != 0) continue;
                do {
                    v /= static_cast<u64>(p);
                } while (v % static_cast<u64>(p) == 0);
                last_small[static_cast<std::size_t>(k)] = static_cast<u32>(p);
            }
        };

        process_residue(k1);
        if (k2 != k1) process_residue(k2);
    }

    u64 ans = 0;
    for (int k = 1; k <= kmax; ++k) {
        const u64 tail = rem[static_cast<std::size_t>(k)];
        u64 lpf = static_cast<u64>(last_small[static_cast<std::size_t>(k)]);
        if (tail > lpf) lpf = tail;
        ans += lpf % kMod18;
        if (ans >= kMod18) ans %= kMod18;
    }
    return ans % kMod18;
}

}  // namespace

int main() {
    assert(solve_case(1) == 5ULL);
    assert(solve_case(2) == 22ULL);
    assert(solve_case(10) == 1070ULL);
    assert(solve_case(100) == 433752ULL);

    const u64 ans = solve_case(10'000'000);
    std::cout << std::setw(18) << std::setfill('0') << ans << "\n";
    return 0;
}

Python

MOD18 = 1000000000000000000

def tonelli_shanks(n, p):
    if n == 0: return 0
    if p == 2: return n
    if pow(n, (p - 1) // 2, p) != 1: return 0
    if p % 4 == 3: return pow(n, (p + 1) // 4, p)

    q = p - 1
    s = 0
    while (q & 1) == 0:
        q >>= 1
        s += 1

    z = 2
    while pow(z, (p - 1) // 2, p) != p - 1:
        z += 1

    m = s
    c = pow(z, q, p)
    t = pow(n, q, p)
    r = pow(n, (q + 1) // 2, p)

    while t != 1:
        tt = t
        i = 0
        while tt != 1:
            tt = (tt * tt) % p
            i += 1
        
        b = pow(c, 1 << (m - i - 1), p)
        r = (r * b) % p
        c = (b * b) % p
        t = (t * c) % p
        m = i

    return r

def sieve_primes(n):
    composite = bytearray(n + 1)
    primes = []
    for i in range(2, n + 1):
        if not composite[i]:
            primes.append(i)
        for p in primes:
            v = p * i
            if v > n: break
            composite[v] = 1
            if i % p == 0: break
    return primes

def solve_case(kmax):
    limit = 2 * kmax
    primes = sieve_primes(limit)

    rem = [0] * (kmax + 1)
    last_small = [1] * (kmax + 1)
    for k in range(1, kmax + 1):
        rem[k] = 4 * k * k + 1

    for p in primes:
        if p == 2 or (p & 3) != 1: continue

        root = tonelli_shanks(p - 1, p)
        if root == 0: continue

        inv2 = (p + 1) // 2
        k1 = (root * inv2) % p
        k2 = 0 if k1 == 0 else p - k1

        def process_residue(residue):
            start = residue
            if start <= 0: start += p
            for k in range(start, kmax + 1, p):
                v = rem[k]
                if v % p != 0: continue
                while v % p == 0:
                    v //= p
                rem[k] = v
                last_small[k] = p

        process_residue(k1)
        if k2 != k1: process_residue(k2)

    ans = 0
    for k in range(1, kmax + 1):
        tail = rem[k]
        lpf = last_small[k]
        if tail > lpf: lpf = tail
        ans = (ans + lpf) % MOD18

    return "{:018d}".format(ans)

def solve():
    ans_str = solve_case(10000000)
    return ans_str.lstrip('0') or '0'

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

Java

import java.util.ArrayList;

public class Euler659 {
    static final long kMod18 = 1000000000000000000L;

    static long modPow(long base, long exp, long mod) {
        long result = 1 % mod;
        base %= mod;
        while (exp > 0) {
            if ((exp & 1) != 0) {
                // To avoid overflow in long * long % mod, use BigInteger or Russian Peasant
                // Here mod is at most p <= 2*10^7, so base*base easily fits in long.
                result = (result * base) % mod;
            }
            base = (base * base) % mod;
            exp >>= 1;
        }
        return result;
    }

    static long tonelliShanks(long n, long p) {
        if (n == 0)
            return 0;
        if (p == 2)
            return n;
        if (modPow(n, (p - 1) / 2, p) != 1)
            return 0;
        if (p % 4 == 3)
            return modPow(n, (p + 1) / 4, p);

        long q = p - 1;
        int s = 0;
        while ((q & 1) == 0) {
            q >>= 1;
            ++s;
        }

        long z = 2;
        while (modPow(z, (p - 1) / 2, p) != p - 1)
            ++z;

        long m = s;
        long c = modPow(z, q, p);
        long t = modPow(n, q, p);
        long r = modPow(n, (q + 1) / 2, p);

        while (t != 1) {
            long tt = t;
            long i = 0;
            while (tt != 1) {
                tt = (tt * tt) % p;
                ++i;
            }
            long b = modPow(c, 1L << (m - i - 1), p);
            r = (r * b) % p;
            c = (b * b) % p;
            t = (t * c) % p;
            m = i;
        }

        return r;
    }

    static ArrayList<Integer> sievePrimes(int n) {
        ArrayList<Integer> primes = new ArrayList<>();
        byte[] composite = new byte[n + 1];
        for (int i = 2; i <= n; ++i) {
            if (composite[i] == 0)
                primes.add(i);
            for (int p : primes) {
                long v = 1L * p * i;
                if (v > n)
                    break;
                composite[(int) v] = 1;
                if (i % p == 0)
                    break;
            }
        }
        return primes;
    }

    static long solveCase(int kmax) {
        int limit = 2 * kmax;
        ArrayList<Integer> primes = sievePrimes(limit);

        long[] rem = new long[kmax + 1];
        int[] lastSmall = new int[kmax + 1];

        for (int k = 1; k <= kmax; ++k) {
            rem[k] = 4L * k * k + 1L;
            lastSmall[k] = 1;
        }

        for (int p : primes) {
            if (p == 2 || (p & 3) != 1)
                continue;

            long root = tonelliShanks(p - 1, p);
            if (root == 0)
                continue;

            long inv2 = (p + 1) / 2L;
            int k1 = (int) ((root * inv2) % p);
            int k2 = (k1 == 0) ? 0 : (p - k1);

            int start1 = k1;
            if (start1 <= 0)
                start1 += p;
            for (int k = start1; k <= kmax; k += p) {
                long v = rem[k];
                if (v % p != 0)
                    continue;
                do {
                    v /= p;
                } while (v % p == 0);
                rem[k] = v;
                lastSmall[k] = p;
            }

            if (k2 != k1) {
                int start2 = k2;
                if (start2 <= 0)
                    start2 += p;
                for (int k = start2; k <= kmax; k += p) {
                    long v = rem[k];
                    if (v % p != 0)
                        continue;
                    do {
                        v /= p;
                    } while (v % p == 0);
                    rem[k] = v;
                    lastSmall[k] = p;
                }
            }
        }

        long ans = 0;
        for (int k = 1; k <= kmax; ++k) {
            long tail = rem[k];
            long lpf = lastSmall[k];
            if (tail > lpf)
                lpf = tail;

            ans += lpf % kMod18;
            if (ans >= kMod18)
                ans %= kMod18;
        }

        return ans;
    }

    public static String solve() {
        long ans = solveCase(10000000);
        return Long.toString(ans);
    }

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