Non distinct prime factors of larger numbers

Question

I wrote and use this function to produce prime factors of a number:

import numpy as np
from math import sqrt

def primesfrom3to(n):
    """ Returns a array of primes, p < n """
    assert n>=2
    sieve = np.ones(n/2, dtype=np.bool)
    for i in xrange(3,int(n**0.5)+1,2):
        if sieve[i/2]:
            sieve[i*i/2::i] = False
    return np.r_[2, 2*np.nonzero(sieve)[0][1::]+1]    

def primefactors(tgt,verbose=True):
    if verbose: 
        print '\n\nFinding prime factors of: {:,}'.format(tgt)

    primes=primesfrom3to(sqrt(tgt)+1)

    if verbose:
        print ('{:,} primes between 2 and square root of tgt ({:.4})'.
                      format(len(primes),sqrt(tgt))) 

    return [prime for prime in primes if not tgt%prime]

If I call this with the value from Project Euler #3, it successfully produces the list of distinct primes:

>>> print primefactors(600851475143)
Finding prime factors of: 600,851,475,143
62,113 primes between 2 and square root of tgt (7.751e+05)
[71, 839, 1471, 6857]

This agrees with what Wolfram Alpha produces for the prime factors. (And the largest there is the correct answer to Project Euler #3)

Now let's say that I want to factors of that number x 1e6:

>>> print primefactors(600851475143*1000000)
Finding prime factors of: 600,851,475,143,000,000
39,932,602 primes between 2 and square root of tgt (7.751e+08)
[2, 5, 71, 839, 1471, 6857]

For this larger number, Wolfram Alpha produces:

2**6 * 5**6 * 71 * 839 * 1471 * 6857

Is there an easy way to modify my code that I can calculate the magnitude of the 2 and the 5 as prime factors of the bigger number?

(I am interested in the raw code or algorithm of this -- not a pointer to a library that will do it for me, thanks!)

Answer 1

The traditional way to do this is to divide out each prime factor in turn and then recurse on your factorisation method. This will in general be faster than sieving for all of the primes, because you only care about the (few) primes that actually divide your number.

Of course, there are many, many better prime factorisation algorithms than trial division; people usually use something like the quadratic sieve for a wide range of numbers, with Pollard's rho method on the small end and the number field sieve on the large. These are all significantly more complicated.

---

Since you're sieving for all the primes beforehand, you don't care about the efficiency of your algorithm. Given that, it's easiest just to work out the multiplicities ex post facto, which is what @tobias_k has written. You could also break it out into a separate function as

def multiplicity(n, p):
    i = 0
    while not n % p:
        i, n = i+1, n/p
    return i

and then

>>> n = 600,851,475,143,000,000
>>> n = 600851475143000000
>>> factors = [2, 5, 71, 839, 1471, 6857]
>>> [(f, multiplicity(n,f)) for f in factors]
[(2, 6), (5, 6), (71, 1), (839, 1), (1471, 1), (6857, 1)]

Answer 2

Once you have the distinct prime factors, you can do something like this:

factors = []
for f in distinct_prime_factors:
    while n % f == 0:
        factors.append(f)
        n /= f

Now factors will hold the list of all the prime factors.