Project Euler Problem 3: Largest Prime Factor

At first I decided to look up what algorithms are there for this purposes and I found out that Fermat’s Factorization Method is one of the simplest there are so I decided to go and use it. I found the algorithm description and tried to convert it into code like this:

int FermatFactor(int N)
{
    double a, b2;
    if (N % 2)
    { std::cout << "Input is divisible by 2.";
        return -1;
    }
    else
    {
        a = std::ceil(std::sqrt(N));
        b2 = (a*a) - N;
        return a - std::sqrt(b2);
    }
}

This had several problems. First, the even/odd check is inverted, N % 2 is true when N is odd, not even. I didn’t notice it at first :P. Second, there’s no loop - at first I thought about using while loop to do this but I first wanted to write the logic before the loop and then put the logic inside the loop but I, once again forgot to do so - so Fermat’s method needs to keep searching until it finds a perfect square, but I was just checking one value and returning immediately. Third, I was using int which overflows for 600851475143. And fourth — perhaps the most conceptual mistake — I was returning one factor and calling it the answer. That’s not the largest prime factor, that’s just a factor.

Before fixing the code, let me explain the actual math behind Fermat’s method because it’s genuinely elegant.

The core idea is that any odd composite number can be written as a difference of two squares:

\[N = a^2 - b^2\]

And a difference of two squares always factors as:

\[a^2 - b^2 = (a+b)(a-b)\]

So if we can find a and b, we immediately get two factors of N: (a+b) and (a-b).

Why does this always work for odd composites? If N has two factors p and q where N = p*q, we can always construct:

\[a = \frac{p+q}{2}, \quad b = \frac{p-q}{2}\]

Since N is odd, both p and q must be odd, so p+q and p-q are both even which means a and b are always whole numbers. A solution always exists.

Since we can’t just conjure a and b out of thin air, we search for them. Rearranging:

\[b^2 = a^2 - N\]

We start a at \(\lceil\sqrt{N}\rceil\) (the smallest value where \(a^2 - N\) is non-negative) and keep incrementing until \(a^2 - N\) is a perfect square.

Quick concrete example with N = 5959:

So a = 80, b = 21, factors are (80+21) = 101 and (80-21) = 59. Check: 101 * 59 = 5959.

The key insight I was missing: Fermat gives you one split. The factors it returns might not be prime themselves. You need to keep splitting recursively until everything is prime, then pick the largest.

Think of it like a tree:

600851475143
    /       \
  71      8462696823
               /       \
             839     10086647
                      /      \
                    1471      6857  ← answer

Here’s the full solution:

#include <iostream>
#include <cmath>

bool isPrime(long long n) {
    if (n < 2) return false;
    for (long long i = 2; i <= std::sqrt(n); i++)
        if (n % i == 0) return false;
    return true;
}

long long fermat(long long N) {
    if (N % 2 == 0) return 2;
    if (isPrime(N)) return N;
 
    long long a = (long long)std::ceil(std::sqrt((double)N));
    long long b2 = a * a - N;
    long long b = (long long)std::sqrt((double)b2);
 
    while (b * b != b2) {
        a++;
        b2 = a * a - N;
        b = (long long)std::sqrt((double)b2);
    }
    return a - b;
}

long long largestPrimeFactor(long long N) {
    if (isPrime(N)) return N;
 
    long long factor = fermat(N);
    long long other = N / factor;
 
    return std::max(largestPrimeFactor(factor), largestPrimeFactor(other));
}

int main() {
    long long N;
    std::cout << "N: ";
    std::cin >> N;
 
    long long largest = 1;
    while (N % 2 == 0) {
        largest = 2;
        N /= 2;
    }
 
    if (N > 1)
        largest = std::max(largest, largestPrimeFactor(N));
     
    std::cout << largest << std::endl;
    return 0;
}

A few important things to note:

Integer arithmetic over floats. In my original code I used double and compared floating point square roots. For large numbers like 600851475143, a double only has ~15 significant digits of precision and starts making rounding errors. The fix is to compute b = (long long)sqrt(b2) and then check b*b ! b2= — integer comparison, exact every time.

Strip factors of 2 first. Fermat’s method only works on odd numbers. Before doing anything we divide out all factors of 2 in a loop. This also avoids the worst case for Fermat where one factor is very small (like 2 * huge_prime) and a would need to travel very far from \(\sqrt{N}\).

Why trial division for isPrime is fine here. Trial division up to \(\sqrt{n}\) sounds slow, but remember we’re calling isPrime on numbers that have already been split by Fermat repeatedly. They’re much smaller than the original input, so it’s fast enough.

The answer for 600851475143 is 6857.

While I was working through this I had to convince myself why trial division up to \(\sqrt{n}\) is sufficient for isPrime. The reasoning is clean: if d divides n, then n/d also divides n — factors come in pairs. Can both d and n/d be larger than \(\sqrt{n}\) simultaneously? No, because then \(d \cdot (n/d) > \sqrt{n} \cdot \sqrt{n} = n\), but \(d \cdot (n/d) = n\) by definition. Contradiction. So every factor pair has at least one member \(\leq \sqrt{n}\), and checking up to there is sufficient.

Simple once you see it, but worth writing down.