How to find multiplicity of prime factors in order to get number of divisors

Question

You might have guessed that I'm doing project euler #12 by the title. My brute force solution took much too long, so I went looking for optimizations that I could understand.

I'm interested in extending the strategy outlined here

The way I've tried to tackle this is by using the Sieve of Eratosthenes to get prime factors like this:

divs = []

multiples = set()
for i in xrange(2, n + 1):
    if i not in multiples:
        if n % i == 0:
            divs.append(i)
        multiples.update(xrange(2*i, n+1, i))
return divs

This itself is a problem because line 8 will yield an overflow error long before the program gets within the range of the answer (76576500).

Now, assuming I'm able to get the prime factors, how can I find their respective multiplicities efficiently?

Answer 1

Your approach for factorization is far from optimal, even if you limit yourself to relatively simple algorithms (i.e. not Brent's algorithm or anything more advanced).

Each time you find a prime factor, divide by that factor until it is no longer divisible. The number of times you can do that is the multiplicity.

Continue with the quotient after division, not your original number.

Find factors and divide until the remaining quotient is less than the square of your divisor. In that case the quotient is 1 or a prime (the last prime factor with multiplicity 1).

To find the factor it is enough to do trial division by 2 and odd numbers starting from 3. Any non-primes will not be a problem, because its prime factors will already be removed before it is reached.

Use the correct data structure to represent the prime factors together with their multiplicity (a map or multiset).

You can also compute the number of divisors directly, without storing the factorization. Each time you find a prime factor and its multiplicity, you can accumulate the result by multiplying with the corresponding factor from the formula for the number of divisors.

If you need to do many factorizations of numbers that are not too big, you can precompute an array with the smallest divisor for each array index, and use that to quickly find divisors.

Answer 2

My standard prime-numbers script is appended below; it provides the Sieve of Eratsothenes to generate primes, a Miller-Rabin primality test, a function that factors integers using a 2,3,5-wheel and Pollard's rho method, the number-theoretic function sigma that calculates the sum of the x'th powers of the divisors of an integer, using the method that you reference in your post, and a function that computes the aliquot sequence starting from a given integer. Given that script, it is easy to solve Project Euler 12, remembering that sigma with x=0 returns the count of the divisors of an integer:

$ python
Python 2.6.8 (unknown, Jun  9 2012, 11:30:32)
[GCC 4.5.3] on cygwin
Type "help", "copyright", "credits" or "license" for more information.
>>> execfile('primes.py')
>>> factors(76576500)
[2, 2, 3, 3, 5, 5, 5, 7, 11, 13, 17]
>>> sigma(0,76576500)
576
>>> i, t = 1, 1
>>> while sigma(0, t) < 500:
...     i += 1; t += i
...
>>> print t
76576500

You can run the program at http://programmingpraxis.codepad.org/V5LiI8V9, and you'll find lots of prime-number stuff at my blog. Here's the code:

# prime numbers

def primes(n): # sieve of eratosthenes
    i, p, ps, m = 0, 3, [2], n // 2
    sieve = [True] * m
    while p <= n:
        if sieve[i]:
            ps.append(p)
            for j in range((p*p-3)/2, m, p):
                sieve[j] = False
        i, p = i+1, p+2
    return ps

# from random import randint

seed = 17500728 # RIP j s bach

def random(): # float on range [0,1)
    global seed
    seed = (69069 * seed + 1234567) % 4294967296
    return seed / 4294967296.0

def randint(lo,hi): # int on range [lo,hi)
    return int((hi - lo) * random()) + lo

def isPrime(n, k=5): # miller-rabin
    if n < 2: return False
    for p in [2,3,5,7,11,13,17,19,23,29]:
        if n % p == 0: return n == p
    s, d = 0, n-1
    while d % 2 == 0:
        s, d = s+1, d/2
    for i in range(k):
        x = pow(randint(2, n-1), d, n)
        if x == 1 or x == n-1: continue
        for r in range(1, s):
            x = (x * x) % n
            if x == 1: return False
            if x == n-1: break
        else: return False
    return True

# from fractions import gcd

def gcd(a,b): # greatest common divisor
    if b == 0: return a
    return gcd(b, a % b)

def insertSorted(x, xs): # insert x in order
    i, ln = 0, len(xs)
    while i < ln and xs[i] < x: i += 1
    xs.insert(i,x)
    return xs

def factors(n, b2=-1, b1=10000): # 2,3,5-wheel, then rho
    if -1 <= n <= 1: return [n]
    if n < -1: return [-1] + factors(-n)
    wheel = [1,2,2,4,2,4,2,4,6,2,6]
    w, f, fs = 0, 2, []
    while f*f <= n and f < b1:
        while n % f == 0:
            fs.append(f)
            n /= f
        f, w = f + wheel[w], w+1
        if w == 11: w = 3
    if n == 1: return fs
    h, t, g, c = 1, 1, 1, 1
    while not isPrime(n):
        while b2 <> 0 and g == 1:
            h = (h*h+c)%n # the hare runs
            h = (h*h+c)%n # twice as fast
            t = (t*t+c)%n # as the tortoise
            g = gcd(t-h, n); b2 -= 1
        if b2 == 0: return fs
        if isPrime(g):
            while n % g == 0:
                fs = insertSorted(g, fs)
                n /= g
        h, t, g, c = 1, 1, 1, c+1
    return insertSorted(n, fs)

def sigma(x, n, fs=[]): # sum of x'th powers of divisors of n
    def add(s, p, m):
        if x == 0: return s * (m+1)
        return s * (p**(x*(m+1))-1) / (p**x-1)
    if fs == []: fs = factors(n)
    prev, mult, sum = fs.pop(0), 1, 1
    while len(fs) > 0:
        fact = fs.pop(0)
        if fact <> prev:
            sum, prev, mult = add(sum, prev, mult), fact, 1
        else: mult += 1
    return add(sum, prev, mult)

def aliquot(n): # print aliquot sequence
    s, ss, k, fs = n, [n], 0, factors(n)
    print n, k, s, fs
    while s > 1:
        s, k = sigma(1,s,fs) - s, k + 1
        fs = factors(s)
        print n, k, s, fs
        if s in ss: return "cycle"
        ss.append(s)
    return ss.pop(-2)

Answer 3

Borrowing from the other answer:

The number a1^k1*a2^k2*...an^kn has number of factor = (k1+1)*(k2+1)...(kn+1)

You can get the prime numbers below a certain number using the following code: Courtesy of Fastest way to list all primes below N

n = number

def primesfrom2to(n):
    """ Input n>=6, Returns a array of primes, 2 <= p < n """
    sieve = numpy.ones(n/3 + (n%6==2), dtype=numpy.bool)
    for i in xrange(1,int(n**0.5)/3+1):
        if sieve[i]:
            k=3*i+1|1
            sieve[       k*k/3     ::2*k] = False
            sieve[k*(k-2*(i&1)+4)/3::2*k] = False
    return numpy.r_[2,3,((3*numpy.nonzero(sieve)[0][1:]+1)|1)]

primes = primesfrom2to(n).tolist()  # list of primes.
primes = map(int, primes)
factors = {}
for prime in primes:
    n = number
    factor = 0
    while True:
        if n%prime == 0:
            factor += 1
            n /= prime
            factors[prime] = factor
        else: break

factors will give you the multiplicity of the prime factors.

Answer 4

example: 28 = 2^2 * 7 ^ 1 number of factors = (2 + 1) * (1 + 1) = 6

in general a ^ k1 * a ^ k2 .. a ^ kn number of factors = (k1 + 1) * (k2 + 1) ... (kn + 1)