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&gt;=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)&lt;=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&gt;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)

&gt;&gt;&gt;
78618449 = (7^1) x (13^1) x (29^1) x (31^3)
31
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:

&gt;&gt;&gt; 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 &gt; 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 &gt; 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 &gt; 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 &gt; 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.