Largest prime factor of a 12-digit number, the patient way

A walk through Project Euler problem #3 — why the obvious algorithm is also the wrong one, why the textbook fix is structurally necessary and not just a speed trick, and how trial division up to √N collapses from minutes to microseconds once you stop being clever in the wrong place.

Why Python

Project Euler #3 fits in 64-bit integer arithmetic and is solved by a tight inner loop of n % p and n //= p. The bottleneck for any sensible algorithm is well under a millisecond, so language choice is dominated by clarity of explanation, not throughput. Python wins on three concrete points:

1. Arbitrary-precision integers are the default. n % p is exact whether n is 12 digits or 12,000. We don't have to think about uint64 overflow, signed/unsigned mismatches, or __int128. If a reader extends the algorithm to a 30-digit composite, the same code keeps working.
2. math.isqrt(n) ships in the standard library since 3.8. It returns the exact integer square root with no floating-point round-trip — important because int(math.sqrt(n)) is wrong starting around n ≈ 2⁵² and Project Euler problems eventually wander past that.
3. The trace I want to show is iterative. I don't need closures, generators, lifetimes, or async. The algorithm fits in twelve lines of straight-line code; a faster language would have me arguing about u64::checked_mul instead of arguing about primes.

The "but it's slow!" objection doesn't survive contact with the actual workload — see the wall-clock section below. C would shave the runtime from 75 µs to maybe 5 µs. Both are imperceptible.

The problem, restated

Project Euler problem #3 (projecteuler.net/problem=3) asks for the largest prime factor of the integer

That's twelve digits, just over 6 × 10¹¹. The interesting question is not "what is the answer" — any factor calculator will tell you — but how cheaply can we reach it, and why does the cheap way work.

The fundamental fact, and what it earns us

Every integer N ≥ 2 has a unique factorization into primes. Order the primes in that factorization:

The answer is p_k. Two facts shape every reasonable algorithm:

Fact 1 (the √N bound). If a composite number N has any nontrivial factorization, at least one factor is ≤ √N. Proof: if both factors exceeded √N, their product would exceed N. So when we're hunting for divisors, we only ever need to look up to √N — never further. For N = 600,851,475,143 that ceiling is 775,146.

Fact 2 (the reduction lemma). If we've divided out every prime ≤ p from N, and the remaining quotient n still satisfies p*p > n, then n itself is prime (or 1). Proof: any composite n must have a factor ≤ √n; if n had a prime factor q, then q would have to be both > p (since we removed everything ≤ p) and ≤ √n (since n is composite). But q > p ⟹ q² > p² > n, yet q ≤ √n ⟹ q² ≤ n — contradiction. So n is prime.

Fact 2 is doing the real work in this problem. It says the loop terminates not at √N, but at √(N / extracted factors) — a much smaller number, because every extraction shrinks the remaining quotient by a multiplicative factor.

Approaches I considered and rejected

Approach 0: build a sieve of Eratosthenes up to N

Sieve up to N, then test each prime against N. Sieve of size 6 × 10¹¹ needs roughly ~75 GB of RAM at one bit per entry. Out before we start.

Approach 1: trial divide every integer from 2 to N

for p in 2..N: if N % p == 0: largest = p. That's 6 × 10¹¹ iterations. At a generous 10⁹ modulos per second in C, this is 600 seconds — ten minutes, well past Project Euler's one-minute guideline. In Python it's an hour. Killed by Fact 1.

Approach 2: trial divide up to √N, no reduction

Bound the loop at √N. About 387 K iterations of if N % p == 0. Should be fast, and Fact 1 promises correctness. So I wrote it:

# from /labs-output/compare2.py, function lpf_naive_no_reduction_broken
def lpf_naive_no_reduction_broken(n: int) -> int:
    largest = 1
    if n % 2 == 0:
        largest = 2
    p = 3
    limit = math.isqrt(n)
    while p <= limit:
        if n % p == 0:
            largest = p
        p += 2
    return largest

When I ran it on N = 600,851,475,143, it returned 486,847 in about 47 ms. That's wrong.

# transcript: python3 compare2.py
B  naive sqrt-trial, no reduction (BROKEN):  returns=486847  @ 46749.3 us/call
   71 * 6857 = 486847   <- composite divisor below sqrt(N)

What this loop actually computes is "the largest divisor of N that lies below √N", not "the largest prime factor". Composite divisors below √N count too — and one of them (the product of the smallest prime factor and the largest one, when the largest prime factor itself fits below √N) ends up being recorded last and clobbering the prime answers.

This is the moment the algorithm gets interesting. The fix is not "add an is_prime(p) check inside the loop", although that does restore correctness. The fix that also makes the loop fast is to divide N down as we go. Fact 2 says: once you've removed every prime ≤ p from N, the smallest remaining divisor of N must itself be a prime. So if we keep N reduced, every n % p == 0 hit is automatically prime — no separate primality test needed.

Approach 2': trial divide up to √N, no reduction, but verify primality of p

For a fair speed comparison I also implemented the corrected no-reduction version (verify each p is prime before recording, and check the final cofactor):

# from /labs-output/compare2.py, function lpf_no_reduction_correct
def lpf_no_reduction_correct(n: int) -> int:
    def is_prime(k):
        if k < 2: return False
        if k % 2 == 0: return k == 2
        i = 3
        while i * i <= k:
            if k % i == 0: return False
            i += 2
        return True

    largest = 1
    p = 2
    limit = math.isqrt(n)
    while p <= limit:
        if n % p == 0 and is_prime(p):
            largest = p
        p += 1 if p == 2 else 2

    cofactor = n
    for q in range(2, limit + 1):
        while cofactor % q == 0:
            cofactor //= q
    if cofactor > 1:
        largest = max(largest, cofactor)
    return largest

This gives the right answer but takes 99.2 ms per call — the cofactor pass alone is another √N loop. Two slow passes.

Approach 3: trial divide with quotient reduction

This is the one. The post's algorithm. See next section.

The algorithm

# from /labs-output/solver.py
def largest_prime_factor(n: int) -> int:
    """Return the largest prime factor of n (n >= 2)."""
    if n < 2:
        raise ValueError("n must be >= 2")

    largest = 1

    # Pull out every factor of 2.
    while n % 2 == 0:
        largest = 2
        n //= 2

    # Pull out every factor of 3.
    while n % 3 == 0:
        largest = 3
        n //= 3

    # Now n is coprime to 6. Step through candidates 5, 7, 11, 13, 17, 19, ...
    # using a 6k +/- 1 wheel: increments alternate 2, 4, 2, 4, ...
    p = 5
    step = 2
    while p * p <= n:
        while n % p == 0:
            largest = p
            n //= p
        p += step
        step = 6 - step   # toggles 2 <-> 4

    # If a prime > sqrt(original_n_after_reductions) remains, that's it.
    if n > 1:
        largest = n

    return largest

Three things are happening here, and each one is doing real work.

The outer condition p * p <= n. This is Fact 2 in code. n shrinks inside the loop; p grows. We exit the moment p² overtakes the current n, not the original. If the answer is a large prime, n will have been chiseled down to exactly that prime by the time we exit, and the final if n > 1 catches it.

The inner while n % p == 0. When p divides n, we don't just record p once — we keep dividing while we can. For N = 600,851,475,143 each prime appears exactly once so the inner loop runs once per hit, but for N = 84 = 2² × 3 × 7 the inner loop on p = 2 runs twice (84 → 42 → 21). Without this, we'd correctly identify the largest prime factor but the next outer iteration would re-find the same factor.

The 6k ± 1 wheel. After removing factors of 2 and 3, every remaining prime is congruent to 1 or 5 mod 6. So instead of stepping p += 2 (which still visits multiples of 3), we step p += 2, 4, 2, 4, .... The pattern visits 5, 7, 11, 13, 17, 19, 23, 25, 29, 31 — skipping 9, 15, 21, 27 entirely. That's a 33% reduction in candidate primes versus odd-only. The step = 6 - step flips between 2 and 4 cleanly without an extra counter.

You could push this further — a 30-stage wheel that skips multiples of 2, 3, 5 saves another ~20%; precomputing trial divisors via a sieve (Eratosthenes or Atkin) eliminates composite candidates like 25, 35, 49 that the wheel still visits — but at N ≈ 6 × 10¹¹ the wheel is already overkill. The whole solve takes 75 µs.

A worked example

Project Euler's introductory example for problem #3 walks through N = 13,195. I'll do the same trace through this algorithm, so the mechanism is concrete before we run it on the actual N.

# transcript: python3 compare2.py
worked example: factor 13195
  divide by 5: 13195 -> 2639
  divide by 7: 2639 -> 377
  divide by 13: 377 -> 29
  divide by 29: 29 -> 1

The trace as a table, with the loop state on each line:

The final if n > 1 is doing real work — by the time the outer loop exits, the largest prime (29) has never been compared to a candidate p, because p² overshot n before we got there. The cleanup line catches it.

Two things to notice. First, the candidates only went up to 13, not up to √13195 ≈ 114. The bound shrinks as n shrinks; that's the whole point. Second, the prime 29 was never tested by trial division — it was identified as prime by Fact 2 alone, because 17² > 29. For PE3 the same thing happens with the answer: it's not found by trial division, it's found by the cofactor check.

Iteration count on the real N

Same algorithm, instrumented to count outer-loop iterations on N = 600,851,475,143:

# transcript: python3 iters.py
C: outer-loop iterations           = 490
C: final remaining n               = [redacted - that is the PE3 answer]
B': bound = sqrt(N)                = 775146
B': odd-only outer iterations      = ~387573
ratio                              = ~790x

Four hundred and ninety. That's the work to factor a twelve-digit number. The reduction-free baseline does 387,573 iterations of essentially the same n % p check, 790× more work, to reach the same answer.

Where do the 490 come from? The algorithm has to scan the 6k ± 1 candidates between consecutive prime factors. After dividing out the smallest factor p₁, the loop has to scan from p₁ + (2 or 4) up to p₂, do the divide, scan up to p₃=1471, and divide — at which point the remaining n is itself p₄=6857, and the outer condition p² > n immediately falls false because p is just past p₃ and p² > p₄. The if n > 1 clause catches p₄. The grand total of candidates scanned is roughly p₃ × (2/6) (the wheel hits 2 out of every 6 integers), which for PE3 lands at ~490.

Wall-clock and the one-minute rule

# transcript: python3 solver.py
N           = 600851475143
sqrt(N)     ~ 775146
answer      = [redacted]
reps        = 10000
per-call    = 76.84 microseconds
total wall  = 768.430 ms

76 microseconds per call, averaged over 10,000 repetitions. That's four orders of magnitude under Project Euler's one-minute guideline. The first call (including Python startup, module loading, JIT cache warm-up if any) is closer to 200 µs; once the interpreter is warm it stays under 80.

Comparison across all four approaches, on the same N:

Two observations:

• C beats Pollard rho here. That's expected: at twelve digits with small prime factors, trial division is already finding factors by candidate ~1500 and ρ has Birthday-Paradox overhead it can't amortize. Pollard ρ overtakes trial division around 20+ digit semiprimes.
• B' is slower than B, even though B is wrong. B's loop is just n % p; B' adds an is_prime(p) call on every hit and a second √N cofactor pass. Wrong but fast vs. correct but slow.

Sanity tests

# transcript: python3 compare2.py
sanity: known factorisations
  lpf(         2) = 2
  lpf(         3) = 3
  lpf(         4) = 2
  lpf(        12) = 3
  lpf(        13) = 13
  lpf(        84) = 7
  lpf(       600) = 5
  lpf(      9973) = 9973
  lpf(   1000003) = 1000003
  lpf(2147483647) = 2147483647

The last one is the Mersenne prime M₃₁ = 2³¹ − 1. The algorithm walks all the way to p ≈ √M₃₁ ≈ 46,341 without ever hitting a divisor, exits the outer loop with n = M₃₁ untouched, and the if n > 1 cleanup returns it. That's the worst case for trial division — a prime input — and it still finishes in well under a millisecond.

lpf(1000003) exercises the same path: 1,000,003 is prime, the loop scans up to ~1000 candidates without a hit, returns the input untouched.

The Pollard rho aside

For inputs PE3's size, Pollard's rho is overkill but worth knowing about — it's the algorithm that takes over when trial division stops being viable, around 20 digits.

The idea: pick a random function f(x) = x² + c mod n, iterate it starting from a random seed, and track two pointers — one stepping one iteration at a time, one stepping two. Because the iterates eventually cycle (pigeonhole, there are only n possible values), the two pointers eventually land on the same value mod some divisor of n. At that point gcd(|x − y|, n) is a nontrivial factor of n with high probability. Expected cost is O(n^(1/4)) modular operations — not O(√n) — which is why ρ scales to inputs where trial division would need a year.

# from /labs-output/compare2.py
def pollard_rho(n: int) -> int:
    if n % 2 == 0:
        return 2
    while True:
        x = random.randrange(2, n)
        y = x
        c = random.randrange(1, n)
        d = 1
        while d == 1:
            x = (x * x + c) % n
            y = (y * y + c) % n
            y = (y * y + c) % n
            d = math.gcd(abs(x - y), n)
        if d != n:
            return d

Paired with a Miller-Rabin primality test (deterministic up to ~3.5 × 10¹² with the witness set {2, 3, 5, 7, 11, 13}), this is the standard recipe for factoring composites up to maybe 25 digits in interactive time. The is_probable_prime helper in compare2.py implements the deterministic-witness variant — when n is below ~3.5 × 10¹² those six witnesses are enough to prove primality, not just probabilistically suggest it.

On PE3, ρ is correct but loses to trial division (112 µs vs 75 µs). Trial division wins anywhere the smallest prime factor is below ~5000; ρ wins anywhere both prime factors are large.

Putting it all together — the solver

The full solver.py from /labs-output/:

"""
Project Euler #3 — largest prime factor of N = 600851475143.

Strategy: trial division with running quotient. Extract each small prime
fully before incrementing the divisor. Stop when p*p exceeds the remaining
quotient — anything left is itself the (largest) prime factor.
"""

import time
import math

def largest_prime_factor(n: int) -> int:
    if n < 2:
        raise ValueError("n must be >= 2")

    largest = 1

    while n % 2 == 0:
        largest = 2
        n //= 2

    while n % 3 == 0:
        largest = 3
        n //= 3

    p = 5
    step = 2
    while p * p <= n:
        while n % p == 0:
            largest = p
            n //= p
        p += step
        step = 6 - step

    if n > 1:
        largest = n

    return largest

def trace(n: int):
    steps = []
    original = n

    while n % 2 == 0:
        steps.append((2, n, n // 2))
        n //= 2
    while n % 3 == 0:
        steps.append((3, n, n // 3))
        n //= 3

    p, step = 5, 2
    while p * p <= n:
        while n % p == 0:
            steps.append((p, n, n // p))
            n //= p
        p += step
        step = 6 - step

    if n > 1:
        steps.append(("remainder is prime", n, 1))

    return original, steps

if __name__ == "__main__":
    N = 600851475143

    REPS = 10_000
    t0 = time.perf_counter_ns()
    for _ in range(REPS):
        ans = largest_prime_factor(N)
    t1 = time.perf_counter_ns()

    per_call_us = (t1 - t0) / REPS / 1_000
    print(f"N           = {N}")
    print(f"sqrt(N)     ~ {math.isqrt(N)}")
    print(f"answer      = {ans}")
    print(f"reps        = {REPS}")
    print(f"per-call    = {per_call_us:.2f} microseconds")
    print(f"total wall  = {(t1 - t0)/1e6:.3f} ms")

    print()
    print("trace of factor extractions:")
    _, steps = trace(N)
    for p, before, after in steps:
        print(f"  divide by {p}: {before} -> {after}")

Save it, run python3 solver.py, get the answer. The post stops short of printing it; Project Euler's culture is that you produce the number yourself, even when it's three keystrokes away.

Things I would change for bigger inputs

The algorithm in this post is correct for any N up to whatever you can store. It stops being fast when the smallest prime factor of N exceeds maybe a million. The natural upgrades, in order of usefulness:

1. A bigger wheel. A 2-3-5 wheel (period 30, 8 residues coprime to 30) saves another ~20% over 6k ± 1. A 2-3-5-7 wheel (period 210, 48 residues) saves a bit more, though neither reaches 50%.
2. A small-prime sieve. Pre-sieve the primes below, say, 10⁶ with Eratosthenes (fits in 1 MB), and trial-divide only by those. The wheel approach above wastes work on candidates like 35, 49, 77 that aren't prime; pre-sieving eliminates them.
3. Pollard ρ for the remainder. Trial-divide up to some threshold T (say 10⁶), then hand the cofactor to ρ. This is what GMP's factoring routines do for medium inputs.
4. ECM and the quadratic sieve for inputs above ~30 digits. These are different algorithms, not optimizations — the math changes. Out of scope here.

For PE #3 specifically, none of this matters. The smallest prime factor is two digits; the algorithm exits in 490 iterations; we're done in 75 microseconds. The hardest problem in the post was understanding why the naive version returned 486,847.

References

• Project Euler problem #3 — https://projecteuler.net/problem=3
• math.isqrt — Python 3.8+ standard library
• Miller-Rabin deterministic witnesses for n < 3,474,749,660,383 — see the table at https://oeis.org/A014233
• Pollard, J. M., "A Monte Carlo method for factorization", BIT 15 (1975)
• The 6k ± 1 wheel observation goes back to Eratosthenes' optimized sieve; for a modern treatment see Crandall & Pomerance, Prime Numbers: A Computational Perspective, Chapter 3.

Artefacts

The files referenced inline — solver.py, compare2.py, iters.py — are in the download tarball for this post. Each one is standalone Python 3.8+; no third-party packages. Run solver.py to get the answer in 75 µs; run compare2.py to reproduce the timing table.