11 Jun 2011 maragato   » (Master)

Euler 10 in Python

I decided to take on Project Euler’s problem #10. Its statement goes like this:

The sum of the primes below 10 is 2 + 3 + 5 + 7 = 17.

Find the sum of all the primes below two million.

My first attempt used brute force and testing primality using only previously found primes:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
primes = []
 
def is_prime(n):
        if not (n < 2 or any(n % x == 0 
           for x in filter(lambda x: x < math.ceil(n ** 2), primes))):
                primes.append(n)
                return True
        return False
 
sum_primes = 0
n = 1   
while n < 2000000:
    if (is_prime(n)):
        sum_primes += n
    n += 1  
print sum_primes

It worked if you consider an execution time of over 7 minutes reasonable. Clearly a bad solution. So I decided to rethink it a bit and tried again. This time, I decided to use a little creative thinking. For every number I tested for primality, I took the time to mark its multiples as not primes (when I first test, say, 3, I already know that 6, 9, 12, ... n * 3 < 2000000 will not be a prime numbers, so I don’t have to test their primality later.)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
max_primes = 2000000
numbers = [0] * max_primes
 
def is_prime(n):
        if n <= 2:
                return True
        if numbers[n] == 1:
                return False
        c = 2
        m = 0
        while True:
                m = n * c
                if m < max_primes:
                        numbers[m] = 1
                else:
                        break
                c+= 1
 
        #if any(n % x == 0 for x in xrange(2, int(n ** 0.5) + 1)):
        i = 2
        while i < int(n** 0.5 + 1):
                if n % i == 0:
                        return False
                i += 1
 
        numbers[n] = 1
        return True
 
def sum_primes():
        soma = 0
        for i in range (2, max_primes):
                if is_prime(i):
                        soma += i
        return soma
 
print sum_primes()

This is a much better algorithm and it gave the correct result in 54 seconds. Project Euler has a “one-minute rule”, which state that all its problems should be solvable “on a modestly powered computer in less than one minute,” which means the solution applies. Still, it’s a brute force solution and I wanted a better one.

So today I picked my copy of Concrete Mathematics and read about primes. Honestly, I’ve been reading more about primes lately than I even thought I would.

Anyways, Knuth—yes, I know the book has three authors, but I can’t help thinking of it as a Knuth book—talks about several strategies to test primality, but it also mentions sieve algorithms that are used to generate lists of prime numbers. Exactly what I wanted!

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
def prime_list(limit):
        if limit < 2:
                return []
        sieve_size = limit / 2
        sieve = [True] * sieve_size
 
        for i in range(0, int(limit ** 0.5) / 2):
                if not sieve[i]:
                        continue
                for j in range((i * (i + 3) * 2) + 3, sieve_size, (i * 2) + 3):
                        sieve[j] = False
 
        primes = [2]
        primes.extend([(i * 2) + 3 for i in range(0, sieve_size) if sieve[i]])
 
        return primes
 
print reduce(lambda x,y: x+y, prime_list(2000000))

Et voilà! Correct result in 0.456s! The code is based on the description of the Sieve of Eratosthenes found in Concrete Mathematics. Also interesting to note, my idea in the second algorithm above is actually the basis of this sieve algorithm, I was just thinking “backwards.”

Syndicated 2010-08-06 01:07:00 from Roberto Teixeira's blog

Latest blog entries     Older blog entries

New Advogato Features

New HTML Parser: The long-awaited libxml2 based HTML parser code is live. It needs further work but already handles most markup better than the original parser.

Keep up with the latest Advogato features by reading the Advogato status blog.

If you're a C programmer with some spare time, take a look at the mod_virgule project page and help us with one of the tasks on the ToDo list!