Problem 526: Largest Prime Factors of Consecutive Numbers

Project Euler Problem 526 Solution

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

Problem Summary Let $P(m)$ denote the largest prime factor of $m$. For nine consecutive integers define $$S(n)=\sum_{j=0}^{8} P(n+j).$$ The task is to determine the maximum value of $S(n)$ for starting points up to $10^{16}$. A direct scan is hopeless: the search range is enormous, and each trial value of $n$ would require factoring nine nearby integers. The implementations therefore do something more structured: they force several offsets to contain prescribed small factors, then check whether the remaining quotients are prime. Mathematical Approach The central idea is to build arithmetic progressions in which the numbers $n,n+1,\dots,n+8$ already have a favorable factorization pattern. If the large quotients are prime, then those quotients are automatically the largest prime factors we want to sum....

Detailed mathematical approach

Problem Summary

Let $P(m)$ denote the largest prime factor of $m$. For nine consecutive integers define

$$S(n)=\sum_{j=0}^{8} P(n+j).$$

The task is to determine the maximum value of $S(n)$ for starting points up to $10^{16}$. A direct scan is hopeless: the search range is enormous, and each trial value of $n$ would require factoring nine nearby integers. The implementations therefore do something more structured: they force several offsets to contain prescribed small factors, then check whether the remaining quotients are prime.

Mathematical Approach

The central idea is to build arithmetic progressions in which the numbers $n,n+1,\dots,n+8$ already have a favorable factorization pattern. If the large quotients are prime, then those quotients are automatically the largest prime factors we want to sum.

Step 1: Force a Useful Divisibility Pattern

The search works modulo

$$M=16\cdot 9\cdot 5\cdot 7\cdot 11\cdot 13=720720.$$

First impose

$$n\equiv 5 \pmod 9,\qquad n\equiv 1 \pmod 5,\qquad n\equiv 3 \pmod 7.$$

Then $n+4$ is divisible by $9$, $5$, and $7$, so

$$315 \mid (n+4).$$

Next choose one of the two residues

$$n\equiv 1 \pmod{16}\qquad \text{or}\qquad n\equiv 7 \pmod{16}.$$

If $n\equiv 1 \pmod{16}$, then the odd offsets satisfy

$$2\mid(n+1),\qquad 4\mid(n+3),\qquad 2\mid(n+5),\qquad 8\mid(n+7).$$

If $n\equiv 7 \pmod{16}$, the powers of two are redistributed as

$$8\mid(n+1),\qquad 2\mid(n+3),\qquad 4\mid(n+5),\qquad 2\mid(n+7).$$

Because $n\equiv 5 \pmod 9$, both $n+1$ and $n+7$ are also divisible by $3$. So the odd offsets always receive the divisor multiset

$$\{6,4,2,24\},$$

although the exact placement of these four divisors depends on which of the two residues modulo $16$ is chosen.

Step 2: Exclude More Small Prime Divisors with CRT

The construction also requires

$$n\equiv 1\text{ or }2 \pmod{11},\qquad n\equiv 1,2,3,\text{ or }4 \pmod{13}.$$

These choices ensure that none of the nine numbers $n,n+1,\dots,n+8$ is divisible by $11$ or $13$: the corresponding residue intervals never hit $0$ modulo those primes. Combining all independent choices gives

$$2\cdot 2\cdot 4=16$$

admissible residue classes modulo $M$. The Chinese Remainder Theorem turns the congruences above into 16 explicit residues $r$, and every candidate examined by the program satisfies

$$n\equiv r \pmod M$$

for one of those residues.

Step 3: Convert Largest Prime Factors into Quotients

For one representative residue family the nine numbers take the form

$$n+j=d_j q_j,\qquad d=(1,6,1,4,315,2,1,24,1).$$

The mirrored family simply permutes the four odd-offset divisors from $(6,4,2,24)$ to $(24,2,4,6)$. In either case the small divisor $d_j$ contains only primes from $\{2,3,5,7\}$. If the quotient $q_j$ is prime, then $q_j$ is automatically larger than every prime dividing $d_j$, so

$$P(n+j)=q_j=\frac{n+j}{d_j}.$$

Therefore, on a fully successful candidate, the target sum becomes a linear expression in $n$. For the displayed family,

$$S(n)=n+\frac{n+1}{6}+(n+2)+\frac{n+3}{4}+\frac{n+4}{315}+\frac{n+5}{2}+(n+6)+\frac{n+7}{24}+(n+8).$$

The mirrored family has the same coefficient of $n$, because it uses the same reciprocal multiset

$$\left\{1,\frac16,1,\frac14,\frac1{315},\frac12,1,\frac1{24},1\right\}.$$

This matters because the score increases as $n$ increases, so once a residue class is fixed we only care about the nearest admissible candidates below the top of the range.

Step 4: Reparameterize the Search by a Single Index

Let

$$N_0=M\left\lfloor\frac{10^{16}}{M}\right\rfloor.$$

For each admissible residue $r$, the implementations examine candidates of the form

$$n=N_0+r-Mk,\qquad 0\le k<10^7.$$

This converts the original search over $n$ into a one-dimensional scan over $k$. Since every step decreases $n$ by exactly $M$, and the score is affine in $n$ once the divisor pattern is fixed, smaller $k$ means a better candidate inside the same residue class.

Step 5: Mark Forbidden Indices with a Prime Sieve

The sieve uses every prime $p$ with

$$17\le p<10^8.$$

These are exactly the primes not already handled by the residue construction. Because $p\nmid M$, the inverse $M^{-1}\pmod p$ exists. For any offset $j$, the condition that the quotient at that offset be divisible by $p$ is equivalent to

$$N_0+r+j-Mk\equiv 0 \pmod p,$$

so the bad values of $k$ lie in the arithmetic progression

$$k\equiv (N_0+r+j)M^{-1}\pmod p.$$

Thus each prime marks at most one residue class of $k$ for each of the nine offsets. Unmarked indices are precisely the candidates for which none of the nine reduced quotients has a prime divisor between $17$ and $10^8-1$.

Step 6: Why Sieving up to $10^8$ Is Enough

The largest reduced quotient is on the order of $10^{16}$, so its square root is about $10^8$. The smallest reduced quotient comes from dividing by $315$, and

$$\sqrt{\frac{10^{16}}{315}}<10^8.$$

Hence every quotient checked by the program has square root below the sieve limit. If a positive quotient survives all prime divisors up to $10^8$, it has no divisor up to its square root and must be prime. That turns the modular sieve into a full primality certificate for the reduced numbers.

Worked Example: The Embedded Checkpoint at $n=100$

The implementations contain a small verification example before the full search. Direct factorization gives

$$\begin{aligned} P(100)&=5,\quad P(101)=101,\quad P(102)=17,\quad P(103)=103,\\ P(104)&=13,\quad P(105)=7,\quad P(106)=53,\quad P(107)=107,\quad P(108)=3. \end{aligned}$$

Therefore

$$S(100)=5+101+17+103+13+7+53+107+3=409.$$

The same checkpoint scan also verifies that

$$\max_{2\le n\le 100} S(n)=417.$$

This example is small enough to check directly, and it confirms that the largest-prime-factor definition is being interpreted correctly before the residue-class sieve is run on the real bound.

How the Code Works

The C++, Python, and Java implementations all use the same mathematical plan. First they build the 16 admissible residues produced by the CRT constraints. Next they generate all primes from $17$ up to $10^8-1$, and for each such prime they precompute the modular inverse of $M$.

Each worker then processes one subset of the residue classes. For one residue class it allocates a byte array of length $10^7$, one slot for each index $k$. For every prime and every offset $j=0,\dots,8$, it marks the arithmetic progression of indices that would make the corresponding reduced quotient divisible by that prime. After the marking phase, the unmarked indices are the surviving candidates.

Finally the implementation evaluates the quotient-based score for each surviving candidate and keeps the maximum over all residue classes. The Python entry point delegates to the same native search logic, while the Java version implements the same sieve-and-scan strategy directly.

Complexity Analysis

Let $B=10^8$ and $K=10^7$. Building the prime table up to $B$ costs $O(B\log\log B)$ time and $O(B)$ memory. For one residue class, the marking work is

$$O\!\left(9K\sum_{p<B}\frac1p\right)=O(K\log\log B),$$

followed by one linear scan of $K$ flags. Since there are only 16 residue classes, the overall asymptotic cost is still $O(B\log\log B+K\log\log B)$ up to a fixed constant factor. In practice the dominant memory cost is the prime sieve plus one $K$-sized marker array per worker thread.

Footnotes and References

Problem page: https://projecteuler.net/problem=526
Prime factor: Wikipedia — Prime factor
Chinese Remainder Theorem: Wikipedia — Chinese remainder theorem
Modular multiplicative inverse: Wikipedia — Modular multiplicative inverse
Sieve of Eratosthenes: Wikipedia — Sieve of Eratosthenes

Problem 526 source code

C++

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

using u8 = std::uint8_t;
using u32 = std::uint32_t;
using u64 = std::uint64_t;
using i64 = std::int64_t;

constexpr u64 L = 10000000000000000ULL;
constexpr int K = 10000000;
constexpr u64 MOD = 16ULL * 9ULL * 5ULL * 7ULL * 11ULL * 13ULL;
constexpr u64 START = (L / MOD) * MOD;
constexpr std::array<u64, 9> D = {1, 6, 1, 4, 315, 2, 1, 24, 1};

u64 egcd_inv(u64 a, u64 mod) {
    i64 t = 0, nt = 1;
    i64 r = static_cast<i64>(mod), nr = static_cast<i64>(a % mod);
    while (nr != 0) {
        i64 q = r / nr;
        i64 tt = t - q * nt;
        t = nt;
        nt = tt;
        i64 rr = r - q * nr;
        r = nr;
        nr = rr;
    }
    if (r != 1) return 0;
    if (t < 0) t += static_cast<i64>(mod);
    return static_cast<u64>(t);
}

u64 crt_two(u64 m1, u64 r1, u64 m2, u64 r2) {
    const u64 inv = egcd_inv(m1, m2);
    const u64 delta = (r2 >= r1 % m2) ? (r2 - (r1 % m2)) : (r2 + m2 - (r1 % m2));
    const u64 t = (static_cast<__uint128_t>(delta) * inv) % m2;
    return r1 + m1 * t;
}

std::vector<u64> build_residues() {
    std::vector<u64> rems = {1, 7};
    u64 prod = 16;

    for (u64& r : rems) r = crt_two(prod, r, 9, 5);
    prod *= 9;

    for (u64& r : rems) r = crt_two(prod, r, 5, 1);
    prod *= 5;

    for (u64& r : rems) r = crt_two(prod, r, 7, 3);
    prod *= 7;

    {
        std::vector<u64> nxt;
        nxt.reserve(rems.size() * 2);
        for (u64 r : rems) {
            nxt.push_back(crt_two(prod, r, 11, 1));
            nxt.push_back(crt_two(prod, r, 11, 2));
        }
        rems.swap(nxt);
    }
    prod *= 11;

    {
        std::vector<u64> nxt;
        nxt.reserve(rems.size() * 4);
        for (u64 r : rems) {
            nxt.push_back(crt_two(prod, r, 13, 1));
            nxt.push_back(crt_two(prod, r, 13, 2));
            nxt.push_back(crt_two(prod, r, 13, 3));
            nxt.push_back(crt_two(prod, r, 13, 4));
        }
        rems.swap(nxt);
    }

    std::sort(rems.begin(), rems.end());
    rems.erase(std::unique(rems.begin(), rems.end()), rems.end());
    return rems;
}

struct PrimeInv {
    u32 p;
    u32 inv;
};

std::vector<PrimeInv> build_primes_and_inverses() {
    constexpr int LIM = 100000000;
    std::vector<u8> is_prime(LIM, 1);
    is_prime[0] = 0;
    is_prime[1] = 0;
    for (int i = 2; static_cast<i64>(i) * i < LIM; ++i) {
        if (!is_prime[i]) continue;
        for (int j = i * i; j < LIM; j += i) is_prime[j] = 0;
    }

    std::vector<PrimeInv> out;
    out.reserve(6000000);
    for (int p = 17; p < LIM; ++p) {
        if (!is_prime[p]) continue;
        if (p == 2 || p == 3 || p == 5 || p == 7 || p == 11 || p == 13) continue;
        const u64 inv = egcd_inv(MOD % static_cast<u64>(p), static_cast<u64>(p));
        out.push_back({static_cast<u32>(p), static_cast<u32>(inv)});
    }
    return out;
}

u64 solve_residue_block(const std::vector<u64>& rems,
                        const std::vector<PrimeInv>& pinv,
                        int thread_id,
                        int thread_count) {
    std::vector<u8> flag(K, 255);
    u64 best_sum = 0;

    for (int rid = thread_id; rid < static_cast<int>(rems.size()); rid += thread_count) {
        const u64 rem = rems[static_cast<std::size_t>(rid)];
        const u8 rid_tag = static_cast<u8>(rid);

        for (const auto [p, inv] : pinv) {
            u64 base = ((START + rem) % p);
            base = (static_cast<__uint128_t>(base) * inv) % p;
            for (int i = 0; i < 9; ++i) {
                for (u64 k = base; k < K; k += p) flag[static_cast<size_t>(k)] = rid_tag;
                base += inv;
                if (base >= p) base -= p;
            }
        }

        for (u64 i = 0; i < K; ++i) {
            if (flag[static_cast<size_t>(i)] == rid_tag) continue;
            const u64 cand = START + rem - MOD * i;
            u64 sum = 0;
            for (int j = 0; j < 9; ++j) sum += (cand + static_cast<u64>(j)) / D[j];
            if (sum > best_sum) best_sum = sum;
        }
    }

    return best_sum;
}

struct WorkerCtx {
    const std::vector<u64>* rems;
    const std::vector<PrimeInv>* pinv;
    int thread_id;
    int thread_count;
    u64 local_best;
};

void* worker_main(void* ptr) {
    auto* ctx = static_cast<WorkerCtx*>(ptr);
    ctx->local_best = solve_residue_block(*ctx->rems, *ctx->pinv, ctx->thread_id, ctx->thread_count);
    return nullptr;
}

unsigned choose_thread_count(std::size_t tasks) {
    if (tasks <= 1) {
        return 1U;
    }
    long cpu_count = sysconf(_SC_NPROCESSORS_ONLN);
    unsigned threads = (cpu_count > 0) ? static_cast<unsigned>(cpu_count) : 1U;
    if (threads > 8U) {
        threads = 8U;
    }
    if (threads > tasks) {
        threads = static_cast<unsigned>(tasks);
    }
    return (threads == 0U) ? 1U : threads;
}

u64 solve() {
    const auto rems = build_residues();
    const auto pinv = build_primes_and_inverses();
    const unsigned thread_count = choose_thread_count(rems.size());

    if (thread_count == 1U) {
        return solve_residue_block(rems, pinv, 0, 1);
    }

    std::vector<pthread_t> threads(thread_count);
    std::vector<WorkerCtx> ctx(thread_count);
    unsigned created = 0U;
    bool failed = false;

    for (unsigned t = 0U; t < thread_count; ++t) {
        ctx[t] = WorkerCtx{&rems, &pinv, static_cast<int>(t), static_cast<int>(thread_count), 0ULL};
        if (pthread_create(&threads[t], nullptr, worker_main, &ctx[t]) != 0) {
            failed = true;
            break;
        }
        ++created;
    }

    for (unsigned t = 0U; t < created; ++t) {
        pthread_join(threads[t], nullptr);
    }

    if (failed) {
        return solve_residue_block(rems, pinv, 0, 1);
    }

    u64 best_sum = 0ULL;
    for (const auto& c : ctx) {
        if (c.local_best > best_sum) {
            best_sum = c.local_best;
        }
    }
    return best_sum;
}

u64 lpf_small(u64 n) {
    u64 best = 1;
    for (u64 p = 2; p * p <= n; ++p) {
        while (n % p == 0) {
            best = p;
            n /= p;
        }
    }
    if (n > 1) best = n;
    return best;
}

bool checks() {
    if (lpf_small(100) != 5ULL) return false;
    u64 g100 = 0;
    for (u64 n = 100; n <= 108; ++n) g100 += lpf_small(n);
    if (g100 != 409ULL) return false;
    u64 h100 = 0;
    for (u64 x = 2; x <= 100; ++x) {
        u64 s = 0;
        for (u64 n = x; n <= x + 8; ++n) s += lpf_small(n);
        h100 = std::max(h100, s);
    }
    if (h100 != 417ULL) return false;
    return true;
}

int main() {
    if (!checks()) return 1;
    std::cout << solve() << '\n';
}

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;
import java.util.Arrays;
import java.util.Collections;
import java.util.List;
import java.util.stream.Collectors;

public class Euler526 {

    static final long L = 10000000000000000L;
    static final int K = 10000000;
    static final long MOD = 16L * 9L * 5L * 7L * 11L * 13L;
    static final long START = (L / MOD) * MOD;
    static final long[] D = { 1, 6, 1, 4, 315, 2, 1, 24, 1 };

    static long egcdInv(long a, long mod) {
        long t = 0, nt = 1;
        long r = mod, nr = a % mod;
        while (nr != 0) {
            long q = r / nr;
            long tt = t - q * nt;
            t = nt;
            nt = tt;
            long rr = r - q * nr;
            r = nr;
            nr = rr;
        }
        if (r != 1)
            return 0;
        if (t < 0)
            t += mod;
        return t;
    }

    static long crtTwo(long m1, long r1, long m2, long r2) {
        long inv = egcdInv(m1, m2);
        long delta = (r2 >= r1 % m2) ? (r2 - (r1 % m2)) : (r2 + m2 - (r1 % m2));
        // Use exact modulo for delta * inv
        long t = 0;
        long modDiv = m2;
        // BigInteger equivalent without objects, using modulo properties:
        long p_val = (delta % modDiv) * (inv % modDiv); // may exceed long if modDiv > 3e9, but m2 <= 13 here
        t = p_val % modDiv;
        return r1 + m1 * t;
    }

    static List<Long> buildResidues() {
        List<Long> rems = new ArrayList<>(Arrays.asList(1L, 7L));
        long prod = 16;

        for (int i = 0; i < rems.size(); i++)
            rems.set(i, crtTwo(prod, rems.get(i), 9, 5));
        prod *= 9;

        for (int i = 0; i < rems.size(); i++)
            rems.set(i, crtTwo(prod, rems.get(i), 5, 1));
        prod *= 5;

        for (int i = 0; i < rems.size(); i++)
            rems.set(i, crtTwo(prod, rems.get(i), 7, 3));
        prod *= 7;

        List<Long> nxt11 = new ArrayList<>();
        for (long r : rems) {
            nxt11.add(crtTwo(prod, r, 11, 1));
            nxt11.add(crtTwo(prod, r, 11, 2));
        }
        rems = nxt11;
        prod *= 11;

        List<Long> nxt13 = new ArrayList<>();
        for (long r : rems) {
            nxt13.add(crtTwo(prod, r, 13, 1));
            nxt13.add(crtTwo(prod, r, 13, 2));
            nxt13.add(crtTwo(prod, r, 13, 3));
            nxt13.add(crtTwo(prod, r, 13, 4));
        }
        rems = nxt13;

        Collections.sort(rems);
        return rems.stream().distinct().collect(Collectors.toList());
    }

    static class PrimeInv {
        int p;
        int inv;

        PrimeInv(int p, int inv) {
            this.p = p;
            this.inv = inv;
        }
    }

    static List<PrimeInv> buildPrimesAndInverses() {
        int LIM = 100000000;
        byte[] isPrime = new byte[LIM];
        Arrays.fill(isPrime, (byte) 1);
        isPrime[0] = isPrime[1] = 0;
        for (int i = 2; (long) i * i < LIM; i++) {
            if (isPrime[i] == 1) {
                for (int j = i * i; j < LIM; j += i)
                    isPrime[j] = 0;
            }
        }

        List<PrimeInv> pinv = new ArrayList<>(6000000);
        for (int p = 17; p < LIM; p++) {
            if (isPrime[p] == 1) {
                if (p == 2 || p == 3 || p == 5 || p == 7 || p == 11 || p == 13)
                    continue;
                long inv = egcdInv(MOD % p, p);
                pinv.add(new PrimeInv(p, (int) inv));
            }
        }
        return pinv;
    }

    public static void main(String[] args) {
        List<Long> rems = buildResidues();
        List<PrimeInv> pinv = buildPrimesAndInverses();

        long bestSum = rems.parallelStream().mapToLong(rem -> {
            byte[] flag = new byte[K];
            Arrays.fill(flag, (byte) 255);

            for (PrimeInv pi : pinv) {
                long p = pi.p;
                long inv = pi.inv;
                long base = ((START + rem) % p);
                base = (base * inv) % p;
                for (int i = 0; i < 9; i++) {
                    for (long k = base; k < K; k += p)
                        flag[(int) k] = 0;
                    base += inv;
                    if (base >= p)
                        base -= p;
                }
            }

            long localBest = 0;
            for (int i = 0; i < K; i++) {
                if (flag[i] != 0) {
                    long cand = START + rem - MOD * i;
                    long sum = 0;
                    for (int j = 0; j < 9; j++)
                        sum += (cand + j) / D[j];
                    if (sum > localBest)
                        localBest = sum;
                }
            }
            return localBest;
        }).max().orElse(0);

        System.out.println(bestSum);
    }
}