Find largest prime factor

Question

Here is my attempt at solving Project Euler problem 3. It works for test numbers up to 9 digits long, but overflows when I input the real 12-digit behemoth.

I've looked at other methods of solution, but I wanted to challenge myself to code one that doesn't search with unnecessary numbers - only primes - and never any of their factors.

I would sincerely appreciate it if someone with more coding experience could show me where in my program I am taking up a lot of memory/computational efficiency/etc.

#The prime factors of 13195 are 5, 7, 13 and 29.
#What is the largest prime factor of the number 600851475143 ?
import numpy as np
import math

def lpf(n):
    orig=n #Store the original n as a reference.
    ceiling=int(math.floor(np.sqrt(n)))+1 #Don't check for prime factors larger than sqrt(n).
    thelist=[]
    for i in xrange(2,ceiling):
        thelist.append(i) #Make a list of integers from 2 up to floor(sqrt(n))+1.
    plist=[2]  #Initialize a running list of primes.
    pfpairs=[] #Initialize a list to store (prime factor, multiplicity) pairs.
    removed=[] #Initialize a list to store removed prime multiples.
    for p in plist:
        if p>=ceiling: #Again, don't check for prime factors larger than sqrt(n).
            break
        i=2
        pexp=0      #Set the multiplicity of the prime to zero initially.
        if n%p==0:  #If the prime divides n, set the new n to n/p and see if it divides again.
            #print n, p
            pexp=1
            n=n//p
            while True:
                if n%p==0: #As long as there is no remainder after division by p, keep dividing n by p,
                    n=n//p #making sure to redefine n and increment p's multiplicity with each successful divison.
                    pexp+=1
                else:
                    pfpairs.append((p,pexp)) #Once a divison fails, store the prime and its multiplicity.
                    break
        while (p*i)<=ceiling: #As long as the prime multiples dont exceed the ceiling,
            if (p*i) in removed: #if the current multiple has already been removed, move on.
                i+=1
            elif (p*i) not in thelist:
                i+=1
            else:
                removed.append(p*i) #Else, add that multiple to the removed pile, and then remove it from the list.
                thelist.remove(p*i)
                i+=1
        for number in thelist: #After all prime multiples (not including p*1) are deleted from the list,
            if number>p:       #the next number in the list larger than p is guaranteed to be prime, so set
                plist.append(number) #p to this new number and go through the logic again.
                break
    print '%d =' % orig, #Print the prime factorization of n to the user.
    for pair in pfpairs:
        if pair != pfpairs[-1]:
            print '(%d^%d) x' % (pair[0], pair[1]),
        if pair == pfpairs[-1]:
            print '(%d^%d)' % (pair[0], pair[1])
    return max(pfpairs, key = lambda pair: pair[0])[0] #Return the largest prime factor.

print lpf(78618449)

>>>
78618449 = (7^1) x (13^1) x (29^1) x (31^3)
31

Answer 1

The problem is:

What is the largest prime factor of the number 600851475143?

Your solution has a function named lpf, which computes the complete prime factorization of the number and prints it and returns the largest prime factor. That is doing too many things! Your function should do one thing: find the largest prime factor.

Doing that, we can simplify your code a lot. Let's start with:

ceiling=int(math.floor(np.sqrt(n)))+1 #Don't check for prime factors larger than sqrt(n).
thelist=[]
for i in xrange(2,ceiling):
    thelist.append(i) #Make a list of integers from 2 up to floor(sqrt(n))+1.

You imported numpy just for sqrt, but math already has it. Is there a reason for that? int also truncates already, so you don't need floor on top of that. Furthermore, why loop through xrange when we can just use range? So a simpler one would be:

ceiling = int(math.sqrt(n)) + 1
thelist = range(2, ceiling)

Now, as saulspatz points out and I totally missed, sqrt(n) is a ceiling for determining if a number is prime. But it's not a ceiling for determining all the factors. For example, the largest prime factor of 10 is 5... but sqrt(10) ~ 3.16. A valid ceiling would be:

ceiling = n/2 + 1

We can refine the ceiling once we find the first prime factor as an optimization.

Next, what are you actually doing with thelist? It looks like you're sort of doing the Sieve of Eratosthenes. But not really. You're doing all sorts of removal from this list, and that's not entirely necessary. You could just use the sieve directly:

def lpf_me(n):
    largest = None
    ceiling = n/2 + 1 
    thelist = [True] * ceiling
    for i in xrange(2, ceiling):
        if thelist[i]:                        # i is prime
            for j in xrange(2*i, ceiling, i): # turn off multiples of i
                thelist[j] = False
            if n % i == 0:
                largest = i                   # new largest prime factor
    return largest

Which when run:

>>> lpf(78618449)
31

Now, this is valid, but pretty slow (run with timeit, 10 times):

me                        148.2s 

That's pretty bad. Can we build on this? Sure we can!

Let's go back to the drawing board with the sqrt(n). Any factor larger than that is either prime or divisible by something less than sqrt(n). But we've already put in the work to find all the primes less than sqrt(n)! So we can just reuse that work:

def lpf_me(n):
    largest = None
    ceiling = int(math.sqrt(n)) + 1 
    thelist = [True] * ceiling
    others = []

    for i in xrange(2, ceiling):
        if thelist[i]:                        # i is prime
            for j in xrange(2*i, ceiling, i): # turn off multiples of i
                thelist[j] = False
            if n % i == 0:
                largest = i                   # new largest prime factor
                other = n/i 
                if other > ceiling:           # hold onto this for later
                    others.append(other)

    for factor in others:
        for i in xrange(2, ceiling):
            if thelist[i] and factor % i == 0: # only have to check divisibility
                break                          # by primes.
        else:
            # this is prime and necessarily larger than largest
            return factor

    return largest

How fast is this? Pretty speedy. There aren't going to be very many candidates, and we can run through those pretty quickly:

up to n/2       148.22s
up to sqrt(n)
- holding on      0.02s
OP solution      21.83s ** note this is incorrect as stops at sqrt(n)
                        ** instead of n/2

I tried running on the original Euler problem value but yours never finished. But, to be fair, mine does less work. So let's actually do the prime factorization. We'll return a list of pairs like you did in the original:

def prime_factorization(n):

Now, again, we can use the Sieve approach. But instead of modifying n throughout, we can just write a helper function to do it for us:

def factor(num, p):
    exp = 0
    while num % p == 0:
        exp += 1
        num /= p
    return (p, exp) if exp > 0 else None

Simple. We don't have to touch n, keep an orig, or pollute our main function with extraneous logic. The flow of the function is the same - same sieve, except instead of keeping track of largest, we keep track of all_factors:

def factorize(n):
    def factor(num, p): 
        exp = 0 
        while num % p == 0:
            exp += 1
            num /= p
        return (p, exp) if exp > 0 else None

    ceiling = int(math.sqrt(n)) + 1 
    thelist = [True] * ceiling
    all_factors = []
    others = []

    for i in xrange(2, ceiling):
        if thelist[i]:                        # i is prime
            for j in xrange(2*i, ceiling, i): # turn off multiples of i
                thelist[j] = False

            next_factor = factor(n, i)
            if next_factor:
                all_factors.append(next_factor)
                other = n/i 
                if other > ceiling:           # hold onto this for later
                    others.append(other)

    for f in reversed(others):
        for i in xrange(2, ceiling):
            if thelist[i] and f % i == 0: # only have to check divisibility
                break                          # by primes.
        else:
            # this is prime
            all_factors.append(factor(n, f))


    return all_factors

Now it's a fair comparison, but mine still runs in 0.027s. It's all about algorithmic complexity.

Summary

• Avoid erasing from lists. It's bad.
• Break large functions up into smaller functions. It'll make your code easier to read and, in this case, even make the logic a lot simpler.
• When you need to solve a problem, just solve the problem. Worry about extraneous similar problems when, and if, they come up.

Answer 2

Your code, and the accepted solution, bath suffer from the misapprehension that the largest prime factor of a number does not exceed the square root of that number. This happens to be true of the number in the Project Euler problem, but it's not true in general. (The largest prime factor of 10 is 5, but the square root of 10 is about 3.2).

You are correct though, that you will have to find the complete factorization of the number (or come pretty close) in order to answer the question in general.

Here's one way to do it. First we guess that the largest prime factor will be less than one million, and we use the sieve of Eratosthenes to find all the primes less than one million. Here it's okay to use the square root. A composite number must have SOME prime factor less than or equal to its square root.

import math
import functools

given=600851475143
x = given

limit = 1000000
primes =  limit * [True]   # limit
primes[0] = False
primes[1] = False
for n in range(2, int(math.sqrt(limit)+1)):
    if primes[n]:
        for m in range(2*n, limit, n):
            primes[m] = False

Now that we know all the primes, we can check which of them divide the given number, and we divide out the large power of that prime that divides the given. If we have found all the primes factors, we'll eventually get down to 1.

primes = [i for i in range(limit) if primes[i]]
factors = []

for p in primes:
    while x > 1:
        quot, rem = divmod(x, p)
        if rem != 0:
            break
         else:
            x = quot
            largest = p
            factors.append(p)

Technically, it isn't necessary to keep track of the factors, or even to reduce x down to 1. Once x <= limit we know we have found all the factors, and we could simply take

max([p for p in primes if x % p == 0])

but of course, it's a good idea to check that we really have found all the factors.

if x > 1:
    print("Increase limit")
else:
    print(largest)
    product = functools.reduce(lambda x, y:x*y, factors, 1)
    assert product == given

The accepted solution gives the correct answer to the Project Euler problem, but that's fortuitous, because all the prime factors are small.

EDIT: Now that the accepted solution has been changed, it's better than this one. There can be at most one prime factor greater than the square root, so we can make the limit the square root of the given to begin with. Then if x is not 1 when when we finish dividing out the prime factors, the largest prime factor must be x. If x is 1, we proceed as above.

Answer 3

I've looked at other methods of solution, but I wanted to challenge myself to code one that doesn't search with unnecessary numbers - only primes - and never any of their factors.

I think this is a fine restriction, but I'd recommend carving up your code into at least two functions: one to find primes by searching over integers, and another to find factors by searching over primes.

def find_primes_up_to_n(n):
    # code goes here
    # ...
    return list_of_primes

def find_prime_factors_of_n(n):
    # code goes here
    # ...
    return list_of_prime_factors

Having separate functions separates the multiple algorithms you want to use more cleanly. This is valuable not only for your own thinking, but also if you want to profile your code to isolate specific time- or memory-intensive steps. A Python line-profiling tool that I use and recommend is called line_profiler.

For bonus points, you could use yield instead of return in both functions, which would lower memory requirements, although this could require modifications to the algorithms you use.

Answer 4

I think the accepted solution is not taking full advantage of the known bounds to the problem, so I'm going to give it a shot as well. One of the nice things of factorization is that, despite the initial bounds not being too optimistic, you can typically improve them as you go.

If you do trial division of a number \$n\$ by the sequence of increasingprime numbers, you will know that the number is itself as prime if it is not divisible by all primes smaller than or equal to its square root. Similarly, if you remove all smaller prime factors from a number, you know that what's left is a prime itself, once the next trial prime is larger than its square root. So despite the fact that the largest possible prime factor of a number is the number itself, you only need to do trial division by primes smaller than or equal to the square root of the remaining factor.

If you couple this with a sieve of Erathostenes like algorithm, you will find yourself only needing to sieve for prime factor smaller than the quartic root of the number.

The whole idea can be implemented very compactly:

def largest_prime_factor(n):
    max_prime_factor = int(n**0.5 + 0.5)
    # 0 and 1 are not primes
    sieve = [False, False] + [True] * (max_prime_factor - 1)
    lpf = None
    prime = 0
    while n > 1:
        while prime < len(sieve) and not sieve[prime]:
            prime += 1
        prime2 = prime * prime
        if prime2 > n:
            # This means that n, or what's left of it after removing
            # smaller prime factors, is a prime itself.
            lpf = n
            break
        # Remove all instances of this prime factor from n
        while not n % prime:
            n //= prime
            lpf = prime
        if prime2 * prime2 <= n:
            # We need to sieve.  We really do not need to go up to
            # max_prime_factor, and could instead do it to the sqrt of
            # n only, not sure it's worth the effort.
            for idx in range(prime2, max_prime_factor + 1, prime):
                sieve[idx] = False
        prime += 1
    return lpf

There are a number of obvious optimizations that would probably make the code more obscure, like keeping only odd primes in the sieve, or dynamically updating the limit up to which sieving is done. But they would not affect the overall complexity of the algorithm.

Answer 5

If you are going for running time, but you do not need the algorithm to be deterministic (you only need it to be "good enough"), you may go with the Miller-Rabin primality test. If you need more about its reliability, you can read more about it at SO: How many iterations of Rabin-Miller should I use for cryptographic safe primes?

Answer 6

Well, you'd do better to store your lists of primes and not primes as sets, rather than lists. Lists have an \$O(n)\$ access time, whereas a set should have an \$O(1)\$ access time

Also it seems to me that once you have done all the tests for 2, you do them for 2 and 3, then 2 and 3 and 5. But you've already determined that n is not divisible by 2 (I haven't had time to copy this and fiddle with it so I may be way off here). A while loop might be better than a for p in plist.

Answer 7

The usual way to implement the algorithm you're using (the sieve of Eratosthenes to find all primes up to a chosen bound) is to make thelist an array of boolean values, so that thelist[i] == True means i is in the list, and thelist[i] == False means that i is not in the list.

The point is that removing an element from the list only costs O(1) time, so it can be done much faster than erasing an element in the middle of the list.

(and also, you don't even have to bother checking if the element is in the list before you remove it)