I just finished playing around with Problem 3 from Project Euler. My rough idea for the problem was to build a list of the primes up to the number in question and start testing them, since I already knew the basics of the Sieve of Eratosthenes and I didn’t know any good algorithms for factoring. After I finished, I looked up some of those algorithms (see Wikipedia’s (incomplete) list), and I probably wouldn’t have bothered with anything that complex anyway. Here’s my original “solution” (I use scare quotes since I know it has bugs, and I knew it at the time).

"""Solves Problem 3 from Project Euler."""

import math

def e_sieve(upper_bound):
    """Uses the Sieve of Eratosthenes to get a list of the primes up to max."""
    primes = []
    candidates = range(2, upper_bound)
    while candidates:
        head = candidates[0]
        primes.append(head)
        candidates = [n for n in candidates[1:] if n % head]

    return primes

def find_highest_prime_factor(to_factor):
    """Find the highest prime factor of to_factor."""
    highest_prime = None
    primes = e_sieve(int(math.sqrt(to_factor)))
    for prime in primes:
        if not to_factor % prime:
            quotient = to_factor / prime
            if quotient in primes:
                return quotient
            highest_prime = prime

    if highest_prime:
        return highest_prime

    return to_factor

if __name__ == "__main__":
    print find_highest_prime_factor(600851475143)

After I coded up my solution and verified I had the right answer (after what seemed like forever but was really ~15 minutes), I just knew there had to be a faster way. A little work with cProfile quickly revealed I was spending all my time in the Sieve. It turns out that even a “naïve” algorithm for factoring is better than generating the primes for a number this size. My new, “naïve” solution ran in 44 ms.

"""Solves Problem 3 from Project Euler."""

def factorize(to_factor):
    """Use trial division to factorize to_factor and return all the resulting \
    factors."""
    factors = []
    divisor = 2
    while (divisor < to_factor):
        if not to_factor % divisor:
            to_factor /= divisor
            factors.append(divisor)
            # Note we don't bump the divisor here; if we did, we'd have
            # non-prime factors.
        elif divisor == 2:
            divisor += 1
        else:
            # Trivial optimization: skip even numbers that aren't 2.
            divisor += 2
    if not to_factor % divisor:
        # Don't forget the last factor
        factors.append(to_factor)
    return factors

if __name__ == "__main__":
    print max(factorize(600851475143))

As an added bonus, this code has fewer bugs. One of these days, I’ll learn to follow KISS. One of these days….

Back to flipping out…