Prime factorization for integers with up to 9-digit long prime factors

Question

I've become interested prime factorization since solving Project Euler problem 3 (finding the largest prime factor of 600851475143). Learning here that initializing lists with many elements to later prune them is the computational equivalent of swimming with bricks, I challenged myself to code a function that returns all prime factors as efficiently as possible, without the need for initializing a big list beforehand.

The routine takes the input n and divides it by 2 as many times as possible before leaving a remainder, redefining n to be the result of the division each time. When 2 no longer cleanly divides n, it moves to 3, but it is after 3 that the crux of the inefficiency arises.

After 3, if n does not equal 1 and there are still factors to be found, the next factor tried is the next consecutive odd integer after 3. For example, take the number \$129373200 = (2^4) × (3^5) × (5^2) × (11^3)\$. After the algorithm divides all 2's, 3's and 5's away, the next test-number is 7, which I view as efficient because 7 is prime. However, as the routine iterates, 9 is tested before 11, and I see this as inefficient because 9 is composite.

If I edit the code such that it stores all tested numbers in a list, and checks if the next test-number is a multiple of a previously tested one or not before iterating, the runtime slows down considerably. Is there an efficient way to do this sort of check via elementary functions without storing tested numbers in a list (or set, dictionary, etc.)?

TLDR: I want to understand why factoring numbers with 8-digit-long prime factors takes the code 7 seconds versus numbers with 9-digit-long prime factors that take nearly 2.5 minutes to solve. I want to reduce this large jump in runtime.

def pf(n):
    startTime=datetime.now()
    factors = [] #Initialize a list to store prime factors.
    while n % 2 == 0:  #While n/2 continues to yield no remainder:
        factors += [2] #Append 2 to the factor list.
        n /= 2         #Redefine n as n/2.
    if n == 1:         #If n is 1, all of its prime factors have been found,
        print datetime.now() - startTime
        return factors #so return the factor list to the user.
    p = 3 #Initialize a count at 3, the next prime after 2
    while p*p <= n:        #While n is greater than or equal to p*p:
        if n%p == 0:       #If p divides n:
            factors += [p] #Append p to the factor list.
            n /= p         #Redefine n as n/p.
        else:       #If p doesn't divide n:
            p += 2  #See if the next consecutive odd number up divides n.
                    #Once all smaller factors are found, and n is smaller than p*p,
    factors += [n]  #append n to the factor list,
    print datetime.now() - startTime,
    return factors  #and then return it to the user.

print pf(121)
print pf(42768)
print pf(19440)
print pf(97200)
print pf(129373200) #inefficiency example
print pf(600851475143)
print pf(31610054640417607788145206291543662493274686990) #consecutive primes
print pf(4383898882371133212190175441147530134182228613257) #5-6 digit primes
print pf(815145012617325671714771027149) #8-digit primes
print pf(9657874875862260078751562987967607300225789) #9-digit primes

Output:

0:00:00 [11, 11]
0:00:00 [2, 2, 2, 2, 3, 3, 3, 3, 3, 11]
0:00:00 [2, 2, 2, 2, 3, 3, 3, 3, 3, 5]
0:00:00 [2, 2, 2, 2, 3, 3, 3, 3, 3, 5, 5]
0:00:00 [2, 2, 2, 2, 3, 3, 3, 3, 3, 5, 5, 11, 11, 11]
0:00:00 [71, 839, 1471, 6857L]
0:00:00 [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113L]
0:00:00.023000 [37489, 59617, 63577, 63841, 74771, 75521, 82217, 97283, 102181, 104717L]
0:00:06.935000 [15485867, 32452843, 32452867, 49979687L]
0:02:24.362000 [122949829, 314606891, 393342743, 674506111, 941083987L]

Answer 1

Code Comments

First, that is a thoroughly excessive amount of commenting. I appreciate your zeal in documentation, but you do not need to comment every line. Just a couple comments here and there is thoroughly sufficient.

Second, your function should do its job. If you want to time it, that is an external job. Use the timeit library:

def pf(n):
    # stuff that doesn't involve datetime
    return factors

print timeit.timeit('pf(val)',
                    setup='from __main__ import pf; val = 145721*(145721+2)',
                    number=10)

That will print how long it took your function to run 10 times. Great library, highly recommended.

Thirdly, factors += [p] should be spelled factors.append(p).

Algorithm Comments

The reason your code performs so poorly is because of the its algorithmic complexity. You are trying every odd number from 2 to sqrt(N), so your algorithm is O(sqrt(N)). But here N is the actual number... if we express N in terms of digits, it's O(sqrt(10^N)) = O(10^(N/2)). So when we go from 30-digit number to a 43-digit number, we would expect to see an increase in time on the order of a million in the worst case (if we were looking at twin primes). The fact that we see an increase of "only" 20x is an artifact of case selection more than anything else. I tried picking multiples of (arbitrary) twin primes, and the runtimes I got after 10 runs each were:

145721 * 145723        0.536s
1117811 * 1117813      2.716s
18363797 * 18363799   45.436s

Roughly 10x each time. Yay predictions.

So what can we do to speed this up? When our times are this bad, it's either because we have a really stupid bug or we're using a bad algorithm. There isn't anything I see that will save this particular algorithm. At best, we can improve by changing the iteration to avoid multiples of 3 and 5:

jump = 4
while p*p <= n:
    ...
    else:
        p += jump
        jump ^= 6 #so it alternates between 2 and 4

That gets us to:

145721 * 145723        0.263s
1117811 * 1117813      2.172s
18363797 * 18363799   34.740s

Still exponential, but better constants at least. Not sure how much better we're going to do than that. So let's just google some prime factorization algorithms. First one I came across was Pollard's rho. Simply copying the algorithm gives produces:

def pollard(n):
    from fractions import gcd

    def get_factor(n):
        x_fixed = 2
        x = 2
        cycle_size = 2
        factor = 1

        while factor == 1:
            for count in xrange(cycle_size):
                if factor > 1: break
                x = (x * x + 1) % n
                factor = gcd(x - x_fixed, n)

            cycle_size *= 2
            x_fixed = x

        return factor

    factors = []
    while n > 1:
        next = get_factor(n)
        factors.append(next)
        n //= next
    return factors

That gives me the correct response on all three twin primes, with this performance:

145721 * 145723          0.021s (12.5x improvement)
1117811 * 1117813        0.295s ( 7.4x)
18363797 * 18363799      1.131s (30.7x)
1500450269 * 1500450271  5.733s (est. ~1000x)

Even jumping up to 10-digits, we still see good performance (we'd expect over an hour with the previous algorithm, so we're talking 1000x improvement), and we see better growth to.

Sometimes, you just need a better algorithm.

Answer 2

timing a function

Generally it is not a good idea to combine timing code with your function. There are better ways to determine how much time a function call takes:

1. Use %timeit in IPython
2. Write a wrapper function.

An example of a wrapper function:

def time_pf(n):
  startTime=datetime.now()
  result = pf(n)   # call the original function
  delta = datetime.now() - startTime
  print delta, "pf for n =", n, "returned:", result
  return result

time_pf(4383898882371133212190175441147530134182228613257)

Then your original function pf can still be called when you don't want the timing information.

You an even generalize this wrapper idea and encapsulate it into a decorator. See this article or this other one or this SO question.

better factoring

The algorithm is slow because, well, in some cases you are testing a lot of numbers which are not divisors (or even prime.)

A nice explanation of how some advanced factoring techniques work is given in this "Ask Dr. Math" article. Some of the techniques include:

• Express \$N\$ as \$(x-y)(x+y)\$
• Express \$N\$ as sums of squares in two ways: \$N = a^2 + b^2 = c^2 + d^2\$, where \$0 < a < c < d < b\$
• Pollard P-1 method
• Elliptic Curve Method
• Number Field Sieve
• Fermat's Method
• Continued Fraction Method
• SQUFOF
• … and others

Implementing one of those algorithms is a good exercise.