0

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?

6
  • projecteuler.net/problem=12 seems different problem
    – marcadian
    Jul 3, 2014 at 4:21
  • @marcadian That is the one I'm doing. I guess I should clarify: I'm trying to get the multiplicities of the nth triangle number in order to test if the sum of those multiplicities + 1 exceeds 500. The first triangle number to meet that qualification should then be the answer. Jul 3, 2014 at 4:24
  • I know that works, but the problem I'm having is finding an efficient method to solve for each of those exponents given an arbitrary number of prime factors. Jul 3, 2014 at 4:31
  • I don't think you need the multiplicities. The problem asks for the number of factors, not all of which are always prime.
    – Blender
    Jul 3, 2014 at 4:31
  • @Blender But you can use the multiplicities of the prime factors to find the number of divisors. Jul 3, 2014 at 4:32

4 Answers 4

1

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.

1

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)
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.

0

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)

8
  • Do you have a general way to solve for an arbitrary number of ks? Jul 3, 2014 at 4:33
  • what do you mean? that is arbitrary already, as it is denoted by k1..kn
    – marcadian
    Jul 3, 2014 at 4:35
  • I mean I don't know how to go about getting each k without testing each combination of k (e.g. nesting for loops equal to the number of prime factors and testing those combinations until a^k1 * b^k2 * ... * z^kn == the triangle number in question) <- note edit, sorry Jul 3, 2014 at 4:39
  • @ThroatOfWinter57: (k_1 + 1) * ... * (k_n + 1) is the number of factors, not the number itself. All you need to do is factor the number and keep track of the multiplicities.
    – Blender
    Jul 3, 2014 at 4:48
  • @Blender Sorry, I fixed that in the edit above. I don't know how to get the multiplicities in the first place though, and can't think of a way to factor the number such that factors are repeated according to their multiplicities. Jul 3, 2014 at 4:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.