Problem 642: Sum of Largest Prime Factors

Project Euler Problem 642 Solution

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

Problem Summary For $n \ge 2$, let $P(n)$ denote the largest prime factor of $n$. We must compute $$S(N)=\sum_{n=2}^{N} P(n)$$ for $N=201820182018$, modulo $10^9$. A direct sieve up to $N$ is completely impractical, so the solution reorganizes the sum by largest prime factor and evaluates the resulting prime sums with a Min_25-style prime-sum transform plus a sparse recursion over prime powers. Mathematical Approach Every integer contributes exactly one prime, namely its largest prime factor. The key is to separate that terminal prime from the rest of the factorization and then count the remaining cofactor in a structured way. Step 1: Separate the terminal prime If $P(n)=p$, then $n$ can be written as $$n=p\,m,$$ where every prime factor of $m$ is at most $p$. Conversely, any product $p\,m$ with $p$ prime and all prime factors of $m$ at most $p$ has largest prime factor $p$. Therefore $$S(N)=\sum_{p\le N} p\,\Psi\left(\left\lfloor\frac{N}{p}\right\rfloor,p\right),$$ where $\Psi(x,y)$ is the number of positive integers $\le x$ whose prime factors are all $\le y$. So the whole problem becomes a weighted count of smooth cofactors....

Detailed mathematical approach

Problem Summary

For $n \ge 2$, let $P(n)$ denote the largest prime factor of $n$. We must compute

$$S(N)=\sum_{n=2}^{N} P(n)$$

for $N=201820182018$, modulo $10^9$. A direct sieve up to $N$ is completely impractical, so the solution reorganizes the sum by largest prime factor and evaluates the resulting prime sums with a Min_25-style prime-sum transform plus a sparse recursion over prime powers.

Mathematical Approach

Every integer contributes exactly one prime, namely its largest prime factor. The key is to separate that terminal prime from the rest of the factorization and then count the remaining cofactor in a structured way.

Step 1: Separate the terminal prime

If $P(n)=p$, then $n$ can be written as

$$n=p\,m,$$

where every prime factor of $m$ is at most $p$. Conversely, any product $p\,m$ with $p$ prime and all prime factors of $m$ at most $p$ has largest prime factor $p$.

Therefore

$$S(N)=\sum_{p\le N} p\,\Psi\left(\left\lfloor\frac{N}{p}\right\rfloor,p\right),$$

where $\Psi(x,y)$ is the number of positive integers $\le x$ whose prime factors are all $\le y$. So the whole problem becomes a weighted count of smooth cofactors.

Step 2: Expand each smooth cofactor uniquely

Write the cofactor in increasing prime order:

$$m=r_1^{a_1}r_2^{a_2}\cdots r_t^{a_t},\qquad r_1<r_2<\cdots<r_t\le p.$$

Then every $n\le N$ has a unique representation

$$n=r_1^{a_1}r_2^{a_2}\cdots r_t^{a_t}p,$$

where $p=P(n)$. This uniqueness is the reason the recursion can enumerate factorizations without double counting: primes are introduced strictly in increasing order.

If the largest prime appears with exponent $b\ge 1$, then the cofactor contains $p^{b-1}$, while the final factor $p$ is the one that contributes to the sum.

Step 3: Derive the recursion

Let $p_1<p_2<\cdots$ be the primes. For a fixed index $j$, define

$$\Sigma_j(x)=\sum_{p_j\le q\le x,\ q\text{ prime}} q.$$

Now suppose we have already fixed a prefix made of smaller prime powers, and the remaining quotient bound is $R$. Let $F(R,j)$ be the total contribution from this state when the next admissible prime is $p_j$. Then

$$F(R,j)=\Sigma_j(R)+\sum_{i\ge j,\ p_i^2\le R}\ \sum_{e\ge 1,\ p_i^{e+1}\le R}\left(F\left(\left\lfloor\frac{R}{p_i^e}\right\rfloor,i+1\right)+p_i\right).$$

The first term counts the case where we stop immediately with one terminal prime $q\ge p_j$. The double sum chooses the next prime $p_i$, inserts $p_i^e$ into the smooth part, and leaves one more copy of $p_i$ as the possible largest prime. The extra $+p_i$ is exactly the contribution from numbers whose largest prime factor is $p_i$ itself.

The final answer is

$$S(N)=F(N,1).$$

Step 4: Fast prime sums with a Min_25-style transform

The recursion constantly needs values of $\Sigma_j(R)$, so it must be able to query prime sums up to many different bounds very quickly. The implementation prepares two synchronized table families:

one indexed directly by $x\le \lfloor\sqrt N\rfloor$, and one indexed by quotient values $\left\lfloor N/x\right\rfloor$.

They start from the naive quantities

$$C(x)=x-1,\qquad T(x)=\sum_{k=2}^{x}k=\frac{x(x+1)}{2}-1,$$

and then remove composite contributions prime by prime in the standard Min_25 manner. After preprocessing, the tables can return both the number of primes up to $x$ and the sum of primes up to $x$:

$$\pi(x)=\#\{q\le x:q\text{ prime}\},\qquad \Sigma(x)=\sum_{q\le x,\ q\text{ prime}} q.$$

Interval prime sums are obtained by subtraction:

$$\Sigma_j(R)=\Sigma(R)-\Sigma(p_j-1).$$

Step 5: Why every integer is counted once

Take any $n\le N$ with factorization

$$n=r_1^{a_1}r_2^{a_2}\cdots r_t^{a_t}p^b,\qquad r_1<r_2<\cdots<r_t<p.$$

The recursion follows one and only one path: it picks $r_1^{a_1}$, then $r_2^{a_2}$, and so on, always moving to larger primes. After all smaller primes have been fixed, the terminal prime $p$ is counted either by an interval prime-sum term when $b=1$, or by one of the direct $+p$ terms when $b\ge 2$. No other path can create the same factorization because the prime order is strictly increasing.

Worked Example: $N=10$

The largest prime factors are

$$P(2)=2,\ P(3)=3,\ P(4)=2,\ P(5)=5,\ P(6)=3,\ P(7)=7,\ P(8)=2,\ P(9)=3,\ P(10)=5,$$

$$S(10)=2+3+2+5+3+7+2+3+5=32.$$

The recursion sees the same total as follows. The root prime-sum term gives

$$\Sigma(10)=2+3+5+7=17.$$

Then the branch for $2$ contributes $2$ for $4=2^2$, another $2$ for $8=2^3$, and $3+5$ for $6=2\cdot 3$ and $10=2\cdot 5$. The branch for $3$ contributes $3$ for $9=3^2$. Hence

$$17+2+2+(3+5)+3=32.$$

How the Code Works

The C++, Python, and Java implementations follow the same algorithmic structure. They first set $M=\lfloor\sqrt N\rfloor$ and build two families of tables, one for direct arguments $x\le M$ and one for quotient arguments $\lfloor N/x\rfloor$. These tables are initialized with counts of integers from $2$ onward and with the corresponding arithmetic-series sums.

Next, the implementation sweeps through the primes up to $M$. Each newly discovered prime updates both families of tables so that composite contributions created by that prime are removed. After this preprocessing, the implementation can answer prime-count and prime-sum queries for every bound that appears later in the recursion.

The recursive stage then explores prime powers in increasing prime order. For a current bound $R$, it first adds the sum of all admissible terminal primes in the interval $[p_j,R]$. It then tries each next prime $p_i$ with $p_i^2\le R$, repeatedly divides the bound by $p_i$, recurses on the reduced bound with a stricter lower bound on the next prime, and adds one direct contribution of $p_i$ for each extra copy of that prime. Every addition is reduced modulo $10^9$.

The Python implementation is only a thin bridge to the same compiled algorithm, so the mathematics and the execution strategy are identical in all three languages.

Complexity Analysis

Let $M=\lfloor\sqrt N\rfloor$. The preprocessing stores $O(M)$ values. Its running time follows the usual Min_25-style profile: it is sublinear in $N$ and grows roughly like $O(N^{3/4}/\log N)$ for this kind of frontier update. The recursive prime-power search is sparse because it only continues while $p^2\le R$, so it is far smaller than a full scan up to $N$. The overall memory usage is $O(\sqrt N)$.

Footnotes and References

Project Euler problem page: https://projecteuler.net/problem=642
Largest prime factor: Wikipedia - Largest prime factor
Smooth number: Wikipedia - Smooth number
Prime-counting function: Wikipedia - Prime-counting function
Min_25 sieve overview: OI Wiki - Min_25 sieve

Problem 642 source code

C++

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

namespace {

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

constexpr i64 kMod = 1'000'000'000LL;

inline i64 mod_norm(i64 x) {
    x %= kMod;
    if (x < 0) x += kMod;
    return x;
}

inline i64 mod_add(i64 a, i64 b) {
    a += b;
    if (a >= kMod) a -= kMod;
    return a;
}

inline i64 mod_sub(i64 a, i64 b) {
    a -= b;
    if (a < 0) a += kMod;
    return a;
}

inline i64 mod_mul(i64 a, i64 b) {
    return static_cast<i64>((static_cast<__int128>(a) * b) % kMod);
}

u64 isqrt_u64(u64 n) {
    u64 r = static_cast<u64>(std::sqrt(static_cast<long double>(n)));
    while ((r + 1) * (r + 1) <= n) ++r;
    while (r * r > n) --r;
    return r;
}

inline i64 sum_2_to_mod(u64 x) {
    if (x < 2ULL) return 0;
    const u128 s = static_cast<u128>(x) * static_cast<u128>(x + 1ULL) / 2U;
    return static_cast<i64>((s - 1U) % static_cast<u128>(kMod));
}

struct Solver {
    u64 n = 0;
    u64 m = 0;

    std::vector<i64> pre_cnt;
    std::vector<i64> pre_sum;
    std::vector<i64> hou_cnt;
    std::vector<i64> hou_sum;
    std::vector<int> primes;

    i64 dfs(u64 res, int last, int from_prime) {
        i64 ret = 0;
        if (from_prime > 0) {
            const i64 base = (res > m) ? hou_sum[static_cast<std::size_t>(n / res)]
                                       : pre_sum[static_cast<std::size_t>(res)];
            ret = mod_sub(base, pre_sum[static_cast<std::size_t>(primes[static_cast<std::size_t>(last)] - 1)]);
        } else {
            ret = hou_sum[1];
        }

        for (int i = last; i < static_cast<int>(primes.size()); ++i) {
            const u64 p = static_cast<u64>(primes[static_cast<std::size_t>(i)]);
            if (p * p > res) break;

            u64 nres = res;
            for (u64 q = p; q * p <= res; q *= p) {
                nres /= p;
                ret = mod_add(ret, dfs(nres, i + 1, static_cast<int>(p)));
                ret = mod_add(ret, static_cast<i64>(p % static_cast<u64>(kMod)));
            }
        }
        return ret;
    }

    i64 solve(u64 input_n) {
        n = input_n;
        if (n <= 1ULL) return 0;
        m = isqrt_u64(n);

        pre_cnt.assign(static_cast<std::size_t>(m + 1), 0);
        pre_sum.assign(static_cast<std::size_t>(m + 1), 0);
        hou_cnt.assign(static_cast<std::size_t>(m + 1), 0);
        hou_sum.assign(static_cast<std::size_t>(m + 1), 0);
        primes.clear();
        primes.reserve(static_cast<std::size_t>(m / 10 + 16));

        for (u64 i = 1; i <= m; ++i) {
            pre_cnt[static_cast<std::size_t>(i)] = static_cast<i64>(i - 1ULL);
            pre_sum[static_cast<std::size_t>(i)] = sum_2_to_mod(i);
            const u64 v = n / i;
            hou_cnt[static_cast<std::size_t>(i)] = static_cast<i64>(v - 1ULL);
            hou_sum[static_cast<std::size_t>(i)] = sum_2_to_mod(v);
        }

        for (u64 p = 2; p <= m; ++p) {
            const std::size_t ps = static_cast<std::size_t>(p);
            if (pre_cnt[ps] == pre_cnt[ps - 1]) continue;

            primes.push_back(static_cast<int>(p));

            const u64 p2 = p * p;
            const u64 q = n / p;
            const i64 p_cnt = pre_cnt[ps - 1];
            const i64 p_sum = pre_sum[ps - 1];
            const u64 mid = m / p;
            const u64 end = std::min(m, n / p2);
            const i64 p_mod = static_cast<i64>(p % static_cast<u64>(kMod));

            for (u64 i = 1; i <= mid; ++i) {
                const std::size_t is = static_cast<std::size_t>(i);
                const std::size_t ips = static_cast<std::size_t>(i * p);
                hou_cnt[is] -= (hou_cnt[ips] - p_cnt);
                hou_sum[is] = mod_sub(hou_sum[is], mod_mul(mod_sub(hou_sum[ips], p_sum), p_mod));
            }
            for (u64 i = mid + 1; i <= end; ++i) {
                const std::size_t is = static_cast<std::size_t>(i);
                const std::size_t qi = static_cast<std::size_t>(q / i);
                hou_cnt[is] -= (pre_cnt[qi] - p_cnt);
                hou_sum[is] = mod_sub(hou_sum[is], mod_mul(mod_sub(pre_sum[qi], p_sum), p_mod));
            }
            for (u64 i = m; i >= p2; --i) {
                const std::size_t is = static_cast<std::size_t>(i);
                const std::size_t ip = static_cast<std::size_t>(i / p);
                pre_cnt[is] -= (pre_cnt[ip] - p_cnt);
                pre_sum[is] = mod_sub(pre_sum[is], mod_mul(mod_sub(pre_sum[ip], p_sum), p_mod));
            }
        }

        primes.push_back(static_cast<int>(m + 1ULL));
        return mod_norm(dfs(n, 0, 0));
    }
};

i64 brute(u64 n) {
    std::vector<u64> lp(static_cast<std::size_t>(n + 1), 0ULL);
    for (u64 i = 2; i <= n; ++i) {
        if (lp[static_cast<std::size_t>(i)] != 0ULL) continue;
        for (u64 j = i; j <= n; j += i) {
            lp[static_cast<std::size_t>(j)] = i;
        }
    }
    i64 s = 0;
    for (u64 i = 2; i <= n; ++i) {
        s = mod_add(s, static_cast<i64>(lp[static_cast<std::size_t>(i)] % static_cast<u64>(kMod)));
    }
    return s;
}

}  // namespace

int main() {
    Solver solver;
    assert(solver.solve(10) == 32);
    assert(solver.solve(100) == 1915);
    assert(solver.solve(10'000) == 10'118'280);
    assert(solver.solve(200'000) == brute(200'000));

    constexpr u64 n = 201'820'182'018ULL;
    std::cout << solver.solve(n) << '\n';
    return 0;
}

Python

from __future__ import annotations

import re
import shutil
import subprocess
from pathlib import Path

ANSWER_RE = re.compile(r"answer\s*:\s*(.+)$", re.IGNORECASE)
EQUAL_RE = re.compile(r"=\s*(.+)$")


def parse_output(stdout: str) -> str:
    lines = [line.strip() for line in stdout.splitlines() if line.strip()]
    if not lines:
        return ""
    answers = []
    equals = []
    for line in lines:
        m1 = ANSWER_RE.search(line)
        if m1:
            answers.append(m1.group(1).strip())
        m2 = EQUAL_RE.search(line)
        if m2:
            equals.append(m2.group(1).strip())
    if answers:
        return answers[-1]
    if equals:
        return equals[-1]
    return lines[-1]


def should_skip_cpp_checkpoints(src: Path) -> bool:
    try:
        text = src.read_text(encoding="utf-8", errors="ignore")
    except OSError:
        return False
    return "--skip-checkpoints" in text


def run_cpp(binary: Path, src: Path, root: Path) -> str:
    cmd = [str(binary)]
    if should_skip_cpp_checkpoints(src):
        cmd.append("--skip-checkpoints")

    try:
        return subprocess.check_output(cmd, text=True, cwd=root)
    except subprocess.CalledProcessError:
        return subprocess.check_output(cmd, text=True, cwd=src.parent)


def solve() -> str:
    problem_id = __file__.split("Euler")[-1].split(".")[0]
    root = Path(__file__).resolve().parent.parent
    src = root / "solutionsCpp" / f"Euler{problem_id}.cpp"
    binary = root / "solutionsCpp" / f".euler{problem_id}_py_bridge"

    if not binary.exists() or src.stat().st_mtime > binary.stat().st_mtime:
        compiler = shutil.which("clang++") or shutil.which("g++")
        if not compiler:
            raise RuntimeError("No C++ compiler found (clang++/g++).")
        subprocess.check_call([compiler, "-std=c++17", "-O2", str(src), "-o", str(binary)])

    output = run_cpp(binary=binary, src=src, root=root)
    parsed = parse_output(output)
    if not parsed:
        raise RuntimeError(f"Euler{problem_id} bridge produced empty output.")
    return parsed


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

Java

import java.util.ArrayList;

public class Euler642 {

    static final long kMod = 1000000000L;

    static long mod_norm(long x) {
        x %= kMod;
        if (x < 0)
            x += kMod;
        return x;
    }

    static long mod_add(long a, long b) {
        a += b;
        if (a >= kMod)
            a -= kMod;
        return a;
    }

    static long mod_sub(long a, long b) {
        a -= b;
        if (a < 0)
            a += kMod;
        return a;
    }

    static long mod_mul(long a, long b) {
        return (a * b) % kMod;
    }

    static long isqrt(long n) {
        long r = (long) Math.sqrt(n);
        while ((r + 1) * (r + 1) <= n)
            ++r;
        while (r * r > n)
            --r;
        return r;
    }

    static long sum_2_to_mod(long x) {
        if (x < 2)
            return 0;
        long p1 = x;
        long p2 = x + 1;
        if (p1 % 2 == 0)
            p1 /= 2;
        else
            p2 /= 2;
        p1 %= kMod;
        p2 %= kMod;
        long s = (p1 * p2) % kMod;
        return mod_sub(s, 1);
    }

    static class Solver {
        long n;
        long m;
        long[] pre_cnt;
        long[] pre_sum;
        long[] hou_cnt;
        long[] hou_sum;
        ArrayList<Integer> primes;

        long dfs(long res, int last, int from_prime) {
            long ret = 0;
            if (from_prime > 0) {
                long base = (res > m) ? hou_sum[(int) (n / res)] : pre_sum[(int) res];
                ret = mod_sub(base, pre_sum[primes.get(last) - 1]);
            } else {
                ret = hou_sum[1];
            }

            for (int i = last; i < primes.size(); ++i) {
                long p = primes.get(i);
                if (p * p > res)
                    break;

                long nres = res;
                for (long q = p; q * p <= res; q *= p) {
                    nres /= p;
                    ret = mod_add(ret, dfs(nres, i + 1, (int) p));
                    ret = mod_add(ret, p % kMod);
                }
            }
            return ret;
        }

        long solve(long input_n) {
            n = input_n;
            if (n <= 1)
                return 0;
            m = isqrt(n);

            pre_cnt = new long[(int) (m + 1)];
            pre_sum = new long[(int) (m + 1)];
            hou_cnt = new long[(int) (m + 1)];
            hou_sum = new long[(int) (m + 1)];
            primes = new ArrayList<>((int) (m / 10 + 16));

            for (long i = 1; i <= m; ++i) {
                pre_cnt[(int) i] = i - 1;
                pre_sum[(int) i] = sum_2_to_mod(i);
                long v = n / i;
                hou_cnt[(int) i] = v - 1;
                hou_sum[(int) i] = sum_2_to_mod(v);
            }

            for (long p = 2; p <= m; ++p) {
                if (pre_cnt[(int) p] == pre_cnt[(int) (p - 1)])
                    continue;

                primes.add((int) p);

                long p2 = p * p;
                long q = n / p;
                long p_cnt = pre_cnt[(int) (p - 1)];
                long p_sum = pre_sum[(int) (p - 1)];
                long mid = m / p;
                long end = Math.min(m, n / p2);
                long p_mod = p % kMod;

                for (long i = 1; i <= mid; ++i) {
                    hou_cnt[(int) i] -= (hou_cnt[(int) (i * p)] - p_cnt);
                    hou_sum[(int) i] = mod_sub(hou_sum[(int) i],
                            mod_mul(mod_sub(hou_sum[(int) (i * p)], p_sum), p_mod));
                }
                for (long i = mid + 1; i <= end; ++i) {
                    long qi = q / i;
                    hou_cnt[(int) i] -= (pre_cnt[(int) qi] - p_cnt);
                    hou_sum[(int) i] = mod_sub(hou_sum[(int) i], mod_mul(mod_sub(pre_sum[(int) qi], p_sum), p_mod));
                }
                for (long i = m; i >= p2; --i) {
                    long ip = i / p;
                    pre_cnt[(int) i] -= (pre_cnt[(int) ip] - p_cnt);
                    pre_sum[(int) i] = mod_sub(pre_sum[(int) i], mod_mul(mod_sub(pre_sum[(int) ip], p_sum), p_mod));
                }
            }

            primes.add((int) (m + 1));
            return mod_norm(dfs(n, 0, 0));
        }
    }

    public static String solve() {
        Solver solver = new Solver();
        long ans = solver.solve(201820182018L);
        return Long.toString(ans);
    }

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