Introduction

A while ago in this post, I proved that the sum of divisors of an integer can be computed as a multiplicative function on primes. Thus, knowing the prime factorization of an integer results in a direct algorithmic method for computing the sum of divisors of the integer without any combinatorial hassle.

This little trick, as well as many other factorization techniques (e.g., greatest common divisor of two integers), pops up quite a bit in the Project Euler problems. They pop up so often that I eventually found myself wanting to collect the routines in a Python module instead of copying and pasting the code every time I needed to do factorizations in Python.

So, today I pushed this little module. (There’s no setup script yet. Just place it in your working directory and ‘import primes’ for now.)

There are other ways we could go about this (using Sage to coordinate factorization methods, for instance). And, of course, for the more complicated problems (e.g., higher numbered Project Euler problems), such an approach in Python won’t do. But, for a surprising number of Project Euler problems below 150 or so, this module should help a ton.

Factoring a Range of Integers

First, you should understand that there’s no elliptic curve method and no quadratic sieve. This module is designed to provide factorization tools (sum of divisors, gcd, lcm, etc.) for a large range or list of consecutive integers. While certain fancy methods are fast at factoring large random integers, they aren’t the best approach when factoring a big list of consecutive integers. A typical sieve and some careful planning works better in such a situation, as strange as that may feel.

Examples: Prime Sieve Class

Using IPython, we’ll explore this module a bit and see what we can do with it. Then, we’ll show how it can be used to solve a Project Euler problem. For now, let’s start by creating a prime sieve. We’ll ask Python to compute all primes below one million.

In [1]: import primes   In [2]: time S = primes.Sieve(1000000) CPU times: user 0.12 s, sys: 0.02 s, total: 0.13 s Wall time: 0.12 s

How many primes are now included in this sieve? That is, how many primes are there below one million?

In [3]: S.numPrimes Out[3]: 78498

We can use the sieve to test any integer in the given range to see if it is a prime.

In [4]: S.isPrime(873) Out[4]: False

We can make a list of primes from the given sieve, in case we wish to use an iterable over primes.

In [5]: time L = S.listPrimes() CPU times: user 0.06 s, sys: 0.00 s, total: 0.06 s Wall time: 0.05 s

Now, L is a list containing all prime values less than one million. We can also make a set. The set will be slightly more expensive to create, but will provide faster membership access once created.

In [6]: time M = S.setPrimes() CPU times: user 0.09 s, sys: 0.01 s, total: 0.10 s Wall time: 0.07 s

If you forget what the limit point of the sieve is, it is contained as a variable inside the Sieve class.

In [7]: S.limit Out[7]: 1000000

Examples: RangeFactor Class

The RangeFactor class uses a prime sieve to factor a list of integers. It then makes various methods available that use this collection of factorization records. Let’s start by factoring all natural numbers less than one million.

In [8]: time F = primes.RangeFactor(1000000) CPU times: user 2.13 s, sys: 0.06 s, total: 2.19 s Wall time: 2.18 s

We can get a more detailed time analysis by asking for verbose output.

In [9]: time F = primes.RangeFactor(1000000, verbose=True) primes : creating sieve primes : sieve created in 0.140506982803 seconds primes : creating factor records primes : factor records created in 0.232191085815 seconds primes : populating factor records primes : factor records filled in 1.8158788681 seconds CPU times: user 2.22 s, sys: 0.08 s, total: 2.29 s Wall time: 2.27 s

One of the most simple and obvious things to do now is to ask for the factorization of an integer in the given range. Currently, the factors are in a list with each integer having a dictionary of prime keys and exponent value pairs. This example should make that relationship a bit more clear.

In [10]: F.factors[5] Out[10]: {5: 1}   In [11]: F.factors[50] Out[11]: {2: 1, 5: 2}

The next most obvious thing we could do is to examine the greatest common divisor and least common multiple of two integers. We could also simply ask if two integers are coprime (which is slightly more efficient than computing the gcd or lcm).

In [12]: F.gcd(100,29754) Out[12]: 2   In [13]: F.gcd(100,6634) Out[13]: 2   In [14]: F.gcd(100,500) Out[14]: 100   In [15]: F.lcm(774,384) Out[15]: 49536   In [16]: F.areCoprime(774,384) Out[16]: False

Nothing here so far is too special. Really, we’re just making standard prime tools available more easily. But, now we can use the multiplicative function on primes to compute the sum of divisors of any integers in the given range. (Recall that a perfect number is a number that is the sum of all of its proper divisors.)

In [17]: F.sumDivisors(496) Out[17]: 992   In [18]: F.sumProperDivisors(496) Out[18]: 496   In [19]: F.isPerfect(496) Out[19]: True

Now, we can do some slightly more complicated computations. Let’s sum all the perfect numbers under one million. First, let’s see what they are. Then, we’ll take the sum.

In [20]: for i in range(2,1000000): ....: if F.isPerfect(i): print i ....: 6 28 496 8128

What’s a bit crazy here is that we’re computing the sum of divisors for ALL integers less than one million to check if each integer is perfect or not. All things considered, that’s going on pretty quickly due to our use of that nice multiplicative function.

In [28]: time sum([i for i in range(2,F.sieve.limit) if F.isPerfect(i)]) CPU times: user 1.60 s, sys: 0.03 s, total: 1.62 s Wall time: 1.58 s Out[28]: 8658

A Sample Project Euler Problem

Problem 21 is stated as follows:

Let d(n) be defined as the sum of proper divisors of n (numbers less than n which divide evenly into n).
If d(a) = b and d(b) = a, where a ≠ b, then a and b are an amicable pair and each of a and b are called amicable numbers.

For example, the proper divisors of 220 are 1, 2, 4, 5, 10, 11, 20, 22, 44, 55 and 110; therefore d(220) = 284. The proper divisors of 284 are 1, 2, 4, 71 and 142; so d(284) = 220.

Evaluate the sum of all the amicable numbers under 10000.

Python Solution

This is solved using the just discussed ‘primes’ module with the following code.

#!/usr/bin/python2   import time import primes   start = time.time()   # factor all numbers from 1 to 10000 F = primes.RangeFactor(10000)   # compute the sum of all proper divisors for integers from 1 to 10000 S = [F.sumProperDivisors(i) for i in range(10000)]   # find amicable pairs and add to a set pairs = set() for i in range(2,10000): if i in pairs: pass else: if S[i] < 10000 and S[i]!=i and S[S[i]]==i: pairs.add(i) pairs.add(S[i])   elapsed = time.time() - start   print "Found result %s in %s seconds" % (sum(pairs), elapsed)

When executed, we get the following output.

Found result 31626 in 0.0326828956604 seconds