Project Euler 3: Finding the Largest Prime Factor by Dividing It Out

This is the third problem in the Project Euler series. Prime factorization is a tool this series will keep coming back to, and this problem is the place to get its most basic form, trial division, exactly right. The code is about a dozen lines, but asking properly “why is the value left at the end prime?” and “how many times does the loop run in the worst case?” teaches more than the length suggests.

Spoiler warning: this post contains a full solution and the answer. If you want to try it yourself first, start with the original problem. Posts in this series never carry ads, because the problem text is licensed under CC BY-NC-SA 4.0, which is non-commercial.

Problem

A prime factor of an integer is a prime that divides it exactly. The problem’s example: the prime factors of 13195 are 5, 7, 13 and 29, and the largest of them is 29. The task is to find, in the same sense, the largest prime factor of 600851475143.

The input is a single 12-digit integer, and the output is a single prime.

Background

Primes and composites. An integer of at least 2 whose only divisors are 1 and itself is prime; every other integer of at least 2 is composite. The number 1 is neither, and it has no prime factors at all.

Factorization and Euclid’s lemma. Every integer of at least 2 can be written as a product of primes, and apart from the order of the factors there is only one way to do it (the fundamental theorem of arithmetic). What this post actually uses is the lemma underneath it, Euclid’s lemma: if a prime p divides a product ab, then p divides a or b (or both). So once n is split as n = a × b, every prime factor of n is a prime factor of a or of b.

A composite has a prime factor no larger than √m. If m is composite, write m = a × b with 1 < a ≤ b < m. Since a ≤ b, a × a ≤ a × b = m, so a ≤ √m. Any prime factor p of a satisfies p ≤ a ≤ √m, and p divides m too. Turned around: if an integer m ≥ 2 is not divisible by any prime up to √m, then m is prime. This fact is the heart of the solution.

Integer division and rounding down. For positive integers f and r, with r / f computed as integer division (rounding down, ⌊r / f⌋), the conditions f × f ≤ r and f ≤ ⌊r / f⌋ are equivalent. If f × f ≤ r then f ≤ r / f, and since f is an integer it is also at most r / f rounded down. Conversely, f ≤ ⌊r / f⌋ ≤ r / f, and multiplying by f gives f × f ≤ r. The C++ version uses this equivalence to avoid an overflowing multiplication.

Brute force

Translating the definition directly: look at every integer d from 2 to n, check whether d divides n and whether d is prime, and keep the largest d that passes both.

def is_prime(m: int) -> bool:
    return m >= 2 and all(m % d != 0 for d in range(2, m))

n = 13195
best = 0
for d in range(2, n + 1):
    if n % d == 0 and is_prime(d):
        best = d

On the small example this happily gives 29, but with n = 600851475143 the outer loop alone would run nearly six hundred billion times. The outer loop is already O(n), and a separate primality test for each divisor adds more. Every extra digit in n multiplies the work by ten, so this approach stops being usable as n grows.

Counting down from n and stopping at the first prime divisor does not fix it either: when the largest prime factor is much smaller than n, almost every number in between still has to be checked.

Derivation

Step 1: dividing out a small prime factor shrinks the problem. Divide 13195 by small numbers in turn. It is not divisible by 2 (it is odd) or by 3 (its digit sum is 19), and so not by 4 either. Dividing by 5 gives

13195 = 5 × 2639

Every remaining prime factor of 13195 is now a prime factor of 2639 (Euclid’s lemma), so we can forget the original number and work with the smaller 2639.

Step 2: keep going. 2639 is not divisible by 5 again, so move on to 6, which fails, and then 7:

2639 = 7 × 377

377 is not divisible by 7 again. 8, 9, 10, 11 and 12 all fail, and 13 gives

377 = 13 × 29

The value left is 29.

Step 3: see when we can stop. The next candidate is 14, and 14 × 14 = 196 > 29. Everything from 2 to 13 has been divided out as far as possible, so 29 is not divisible by any prime up to 13. Since √29 ≈ 5.4, every prime up to √29 has already been tried, so by the fact from the background section 29 is prime. There is no need to test anything against 29.

So 13195 = 5 × 7 × 13 × 29, and the 29 left at the end is larger than the 5, 7 and 13 divided out before it, so the answer is 29 — matching the problem’s example.

Step 4: composite candidates take care of themselves. Step 1 tried composites such as 4 and 6 as candidates, and they never divided. That is not luck. If a candidate f is composite, all of its prime factors are smaller than f, and those primes were completely divided out of the remaining value before f’s turn came. So the remaining value cannot be divisible by f. As a result, any candidate that actually divides is always prime, and there is no need to build a separate list of primes to test.

Final algorithm

1. If n < 2, there is no prime factor; return 0 (see “Boundary cases” below).
2. Start with remaining = n, largest = 1 and factor = 2.
3. While factor × factor ≤ remaining: as long as remaining is divisible by factor, divide it and set largest = factor. Then increase factor by 1.
4. When the loop ends, if remaining > 1, set largest = remaining.
5. Return largest.
6. Before printing the answer, check small inputs: 1 → 0, 2 → 2, 9 → 3, 64 → 2, 97 → 97, 13195 → 29. If any of them differs, exit with a nonzero status.

The C++ version writes the condition in step 3 as factor ≤ remaining / factor. As shown in the background section it is the same condition, but without a multiplication it cannot overflow; the integer-range notes below explain why that matters.

Correctness proof

Each time the loop tests its condition with candidate f (f = 2, 3, …), the following invariant holds:

• (a) n = (product of the primes divided out so far) × remaining, and every prime divided out so far is smaller than f.
• (b) remaining is not divisible by any prime smaller than f.
• (c) If anything has been divided out, largest is the largest such prime; otherwise it is 1.

At the start (f = 2) nothing has been divided out, remaining = n and largest = 1, so the invariant holds. When the condition is true and f is processed, the inner loop divides out f until remaining is no longer divisible by it. If f divides at all, it is prime by Step 4, and each division moves one f into the product, so (a) is kept. After the inner loop remaining is not divisible by f, so (b) holds for f + 1. f is larger than every prime divided out earlier, so setting largest = f keeps (c). Each division at least halves remaining, so the inner loop always ends.

In the outer loop f increases by 1 each time and remaining never grows, so eventually f × f > remaining and the loop ends. Call the final candidate F. There are two cases.

• remaining = 1: by (a), n is the product of the primes divided out. By Euclid’s lemma every prime factor of n is one of them, so largest from (c) is the largest prime factor.
• remaining > 1: by (b), remaining is not divisible by any prime smaller than F, and since F × F > remaining, every prime up to √remaining is smaller than F. By the background fact, remaining is prime. Its smallest prime factor — itself — is at least F, so it is larger than every prime divided out, all of which are smaller than F. Every prime factor of n is either a divided-out prime or remaining, so the largest prime factor is remaining.

In both cases the return value is the largest prime factor of n.

Complexity

Worst-case iteration count. The body of the outer loop runs only for f with f × f ≤ remaining ≤ n, so at most ⌊√n⌋ − 1 times. The number of divisions that actually happen in the inner loop equals the number of prime factors of n counted with multiplicity, which is at most log₂ n because each is at least 2. Apart from that, each outer iteration costs one remainder test and one condition test. On fixed-width integers these are constant-time operations, so the time is O(√n), and only three values are stored, so the auxiliary space is O(1).

That bound is actually reached. If n is prime, remaining never shrinks and f has to climb to about √n. Calling this “logarithmic because remaining keeps shrinking” is wrong: the number of divisions is bounded by a logarithm, but the number of candidates skipped over is not.

When it finishes early. The loop stops sooner as remaining gets smaller. List the prime factors with multiplicity in order of size and call the largest p₁ and the second largest p₂. The loop runs to roughly max(p₂, √p₁): it has to reach p₂ to divide it out, and after that only p₁ is left, so it stops once f passes √p₁. The more small prime factors a number has, the earlier it finishes compared with √n. The answer section counts how many iterations this problem’s input actually takes.

Integer range.

• 600851475143 is larger than 2^39 ≈ 5.5 × 10^11 and smaller than 2^40 ≈ 1.1 × 10^12, a 40-bit number. It does not fit in a 32-bit integer (at most about 4.3 × 10^9), so a 64-bit integer is required.
• The C++ code uses unsigned long long. The C++ standard guarantees it is at least 64 bits wide (it follows the C standard’s minimum for ULLONG_MAX), and unsigned arithmetic is defined modulo 2^N, so there is no undefined behaviour. Overflow would still silently give a wrong value. The literal carries the ULL suffix so it has that type from the start.
• Division and remainder inside the loop never make values larger, so they cannot overflow. The risk is the multiplication in factor * factor <= remaining. If remaining were a prime close to 2^64, factor would climb to 2^32, where factor × factor = 2^64 wraps around to 0 and the condition becomes true again, so the loop would fail to stop when it should. The C++ version therefore writes the same condition as factor <= remaining / factor: it uses only division, so it cannot overflow, and by the equivalence in the background section it gives the same result.
• While the condition is true, factor is at most √remaining < 2^32, so ++factor cannot overflow either. The C++ function is therefore exact for every unsigned long long value. A prime close to 2^64, though, makes the loop run about 2^32 ≈ 4.3 billion times, so it is correct but slow.
• Python’s int has unlimited precision and cannot overflow, so the Python version keeps the condition as factor * factor <= remaining. The algorithm is the same.

Boundary cases

• n < 2: 0 and 1 have no prime factors, so “the largest prime factor” is undefined. Both implementations return 0 in that case. 0 is not prime, so it cannot be confused with a real answer. In Python negative numbers get the same treatment; the C++ type is unsigned and cannot hold a negative value. 1 → 0 is among the small-example checks.
• n prime: the loop never divides anything, and n itself is what remains and is returned (2 → 2, 97 → 97). For 2, 2 × 2 > 2 from the start, so the loop body never runs.
• n a prime power: 64 = 2^6 becomes 1 after dividing by 2 six times, and the answer is the last factor divided out, 2. The final remaining > 1 check keeps 1 from being mistaken for the answer.
• n the square of a prime: 9 can only be divided out at f = 3 if the condition 3 × 3 ≤ 9 counts as true. Writing the condition with < by mistake would declare 9 prime and return 9, which is why this case is among the small-example checks.

C++17

#include <iostream>

// Largest prime factor of n. Values below 2 have no prime factor; return 0 for them.
// Divide out every factor from 2 upward while factor * factor <= remaining.
// The loop tests factor <= remaining / factor, the same condition without a
// multiplication that could overflow. Whatever is left above 1 at the end is
// prime and larger than every factor divided out before it.
unsigned long long largest_prime_factor(unsigned long long n) {
    if (n < 2) {
        return 0;
    }
    unsigned long long remaining = n;
    unsigned long long largest = 1;
    for (unsigned long long factor = 2; factor <= remaining / factor; ++factor) {
        while (remaining % factor == 0) {
            remaining /= factor;
            largest = factor;
        }
    }
    if (remaining > 1) {
        largest = remaining;
    }
    return largest;
}

int main() {
    if (largest_prime_factor(1) != 0 || largest_prime_factor(2) != 2 ||
        largest_prime_factor(9) != 3 || largest_prime_factor(64) != 2 ||
        largest_prime_factor(97) != 97 || largest_prime_factor(13195) != 29) {
        return 1;
    }
    std::cout << largest_prime_factor(600851475143ULL) << '\n';
    return 0;
}

The small examples are checked with if, not assert, because assert disappears in builds that define NDEBUG. If any check fails, return 1 gives a nonzero exit status.

Python 3

def largest_prime_factor(n: int) -> int:
    """Largest prime factor of n. Values below 2 have no prime factor; return 0 for them.

    Divide out every factor from 2 upward while factor * factor <= remaining.
    Whatever is left above 1 at the end is prime and larger than every factor
    divided out before it.
    """
    if n < 2:
        return 0
    remaining = n
    largest = 1
    factor = 2
    while factor * factor <= remaining:
        while remaining % factor == 0:
            remaining //= factor
            largest = factor
        factor += 1
    if remaining > 1:
        largest = remaining
    return largest


def main() -> None:
    if (
        largest_prime_factor(1) != 0
        or largest_prime_factor(2) != 2
        or largest_prime_factor(9) != 3
        or largest_prime_factor(64) != 2
        or largest_prime_factor(97) != 97
        or largest_prime_factor(13195) != 29
    ):
        raise SystemExit(1)
    print(largest_prime_factor(600851475143))


if __name__ == "__main__":
    main()

Same algorithm and same small examples as the C++ version. Division uses floor division //= so the value never turns into a float; /= would produce a floating-point result that can lose precision on large numbers. A failed small example raises SystemExit(1) for exit status 1; assert is avoided because running Python with -O removes it.

Answer and checks

The largest prime factor of 600851475143 is 6857.

Here is what the algorithm goes through:

factor = 2 … 70      no division
factor = 71          600851475143 = 71 × 8462696833
factor = 72 … 838    no division
factor = 839         8462696833   = 839 × 10086647
factor = 840 … 1470  no division
factor = 1471        10086647     = 1471 × 6857
factor = 1472        1472 × 1472 > 6857 → stop; the remaining 6857 is prime

So 600851475143 is the product of the four primes in the table. The outer loop body runs 1470 times, for factor = 2 through 1471, and only three divisions actually happen. It stops far earlier than the initial bound ⌊√600851475143⌋ = 775146 because, as the Complexity section explained, the second-largest prime factor, 1471, is small.

The result is checked three ways.

1. Small examples. Before printing the answer, both programs check that 1, 2, 9, 64, 97 and 13195 give 0, 2, 3, 2, 97 and 29, and exit with a nonzero status otherwise. 13195 → 29 was also traced by hand in the derivation.
2. Both implementations agree. In this series’ repository a check script compiles the C++ with g++ -std=c++17 -O2 -Wall -Wextra -Werror (warnings count as failures), runs it, and runs the Python with python3 -I -B. Both outputs must match the stored answer.txt byte for byte; a nonzero exit status or a timeout (1 s for C++, 10 s for Python) is a failure. By the series rules, only code that passes this check is published.
3. By hand. Multiply the four prime factors in the table back together: 71 × 839 = 59569, the last two give 10086647, and 59569 × 10086647 = 600851475143, the original number. Each of the four can be confirmed prime with the background fact; for example, the largest one has a square root of about 82.8 and is not divisible by any of the primes 2, 3, 5, …, 79.

Going further

• Half the candidates. After dividing out all the 2s, only odd candidates need testing. Which part of the correctness proof has to change, and how much does the iteration count drop? What about handling 2 and 3 first and then testing only numbers of the form 6k ± 1?
• All prime factors and exponents. Recording (prime, count) at each division in the same loop gives the whole factorization. That leads to formulas for the number and the sum of divisors, which come up within problems 1–100 of this series as well (for example problems 12 and 21).
• Many numbers at once. If you sieve the smallest prime factor of every number from 1 to N in advance, each number can be factored by following that table with only a few divisions. How would you adapt the sieve of Eratosthenes? A sieve is also useful in problems 7 and 10 of this series, which deal with many primes.
• The worst input. Which 64-bit input makes this code run longest? Compare when the loop stops for a large prime and for a product of two primes of similar size, using max(p₂, √p₁) from the Complexity section.
• An overflow experiment. Change the C++ condition to factor * factor <= remaining and think about what happens for a value near the largest prime below 2^64 (actually running it would take a long time).

Sources and license

• Project Euler, Problem 3 — the original problem. The Problem section of this post is my own summary of it and is shared under the same CC BY-NC-SA 4.0 license. Attribution for the problem: Project Euler.
• Project Euler, Copyright — the license of the problem text and how to link to problems. Because of the non-commercial condition, this post carries no ads.
• The solution, proofs and code are original. No Project Euler overview or editorial, forum post, member-only material or solution site was consulted. Motitory is not affiliated with or sponsored by Project Euler, and Project Euler has not reviewed this post.
• C++17 final working draft N4659 — §6.9.1 fundamental types (unsigned arithmetic modulo 2^N), §21.3.6 <climits> (ranges follow the C standard), §22.3 assert (NDEBUG)
• Python Standard Library: Numeric Types — unlimited precision of int, // (floor division) and / (true division)
• Python Language Reference: The assert statement — assert is removed in optimised mode