# Find largest prime factor

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

• Have you profiled your code? – Curt F. Sep 10 '15 at 14:44
• @yaegermenjensen see stackoverflow.com/q/582336/1310566 and docs.python.org/2/library/profile.html . Profiling is a useful tool to detect what part of your program slows it down. – Simon Forsberg Sep 10 '15 at 15:03
• What makes you think that a number can have no prime factor larger than its square root? The SMALLEST prime factor is less than or equal to the square root. For example, the largest prime factor of 10 is 5, and the square root of 10 is about 3.2. – saulspatz Sep 10 '15 at 15:45
• @yaegermenjensen I was going to mention this in your recently deleted question, but your code has a lot of comments in it. I assume that this is because you're coming from another language where this is more commonplace, and that's fine, but in Python your code should be written such that its purpose is self-evident. You should only really have to comment when the code does something obscure or hard to understand. Here's a link to a Python style guide. – Zenohm Sep 10 '15 at 20:02
• @saulspatz but for efficiency, it's only worth checking those less than or equal to the square root, as once you find 2, you can get 5 by doing (10/2) – user2813274 Sep 10 '15 at 23:52

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.
• You, sir, are a baller of the highest caliber. This type of answer was above and beyond my expectation for help. Walking me through my inefficiencies like this has been an invaluable learning experience. I'm a college student, and I feel as if I went to a CS professor's office hours with the same question, I wouldn't even receive a fraction of a smidgen of the real guidance exhibited in this response. I never thought a random stranger on the internet would be so kind. Thank you sir! I took all of your advice to heart. The only thing I'm still a bit confused on is your timeit reference. – yaegermenjensen Sep 10 '15 at 15:20
• I think you're mistaken about the size of the largest prime factor. The SMALLEST prime factor must be less than or equal to the square root. The largest prime factor of 10, for example, is 5, and the square root of 10 is about 3.2. I admit I haven't tried running your code, so maybe I'm overlooking something. – saulspatz Sep 10 '15 at 15:54
• You only need to run to root of the number. Findthe other multiple from the one under the square root – exussum Sep 10 '15 at 18:36
• You can stop updating your sieve once you get to sqrt(ceiling). For example if n was 10000 so ceiling was 100, then you need to find all primes below 100, so you only need to update your sieve for primes below ten. Multiples of eleven won't give you any new composite numbers below 100, the first one is 11*11=121. – Brian Moths Sep 10 '15 at 18:51
• Note that in: xrange(2*i, ceiling, i) you can replace 2*i with i*i or i**2. (All numbers below i**2 are either prime or have divisors smaller than i that you have already checked). You could also special-case 2 since it's the only even prime and then you could use 2*i as last argument (and avoid iterating over evens in the first place). – Bakuriu Sep 10 '15 at 19:15

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.

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.

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.

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?

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.

• "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." I initially thought of this as well, and to check to ensure this was not happening, I added a debugging 'print p' statement just after the 'for p in plist' line. The console prints consecutive primes up to sqrt(n). I will take your advice and learn a bit about sets though and get back to you. Thank you so much. – yaegermenjensen Sep 10 '15 at 13:29
• aah. my apologies, it's the adding values to the list you are scanning through. other languages won't guarantee that works. it's probably a good idea not to alter the list you're looping over and a while would be easier to read (and avoid some of the list manipulation) – Tom Tanner Sep 10 '15 at 13:37
• There is one additional complication: you need to be able to find the smallest element of thelist greater than p. – user14393 Sep 11 '15 at 7:34

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)