# Sieve of Eratosthenes - Python

I've been doing a lot of Project Euler lately and just wanted to make sure my implementation was as good as it could be. Does anyone have any suggestions to speed this up?

def sieve(upperlimit):
# mark off all multiples of 2 so we can use 2*p as the step for the inner loop
l = [2] + [x if x % 2 != 0 else 0 for x in range(3, upperlimit + 1)]

for p in l:
if p ** 2 > upperlimit:
break
elif p:
for i in range(p * p, upperlimit + 1, 2 * p):
l[i - 2] = 0
# filter out non primes from the list, not really that important i could work with a list full of zeros as well
return [x for x in l if x]

-
For most project euler problems this version of the sieve together with an implementation of the Miller-Rabin primality test should be all you need. I'm not sure saving 2 places in the array is worth the number of times your going to adjust the index with l[i - 2] though. – Jackson Feb 21 '14 at 16:11

1. Here's my starting point for computing the performance improvement due to the various revisions below: how long does it take to sieve for the prime numbers below $10^8$?

>>> from timeit import timeit
>>> test = lambda f: timeit(lambda:f(10**8), number=1)
>>> t1 = test(sieve)


The exact number is going to depend on how fast your computer is, so I'm going to compute performance ratios, but for the record, here it is:

>>> t1
78.9875438772142

2. Your initialization of the list l takes more than half the time, so let's try a cheaper approach. Let's also give this array a better name, and make it a Boolean array while we're about it.

def sieve2(n):
"""Return a list of the primes below n."""
prime = [True] * n
for p in range(3, n, 2):
if p ** 2 > n:
break
if prime[p]:
for i in range(p * p, n, 2 * p):
prime[i] = False
return [2] + [p for p in range(3, n, 2) if prime[p]]


When optimizing a function like this, it's always worth keeping the un-optimized version around to check the correctness of the optimized version:

>>> sieve(10**6) == sieve2(10**6)
True


This already runs in less than a third of the time:

>>> test(sieve2) / t1
0.30390444573149544

3. We could avoid the test for p ** 2 > n by computing a tighter limit for the loop. Note that I've used n ** .5 here as this is slightly faster than math.sqrt(n).

def sieve3(n):
"""Return a list of the primes below n."""
prime = [False, False, True] + [True, False] * (n // 2)
for p in range(3, int(n ** .5) + 1, 2):
if prime[p]:
for i in range(p * p, n, 2 * p):
prime[i] = False
return [p for p in range(2, n) if prime[p]]


This makes little difference to the overall runtime:

>>> test(sieve3) / t1
0.2971086436068156

4. We can accumulate the result as we go, instead of in a separate iteration at the end. Note that I've cached result.append in a local variable to avoid looking it up each time round the loop.

def sieve4(n):
"""Return a list of the primes below n."""
prime = [False, False, True] + [True, False] * (n // 2)
result = [2]
append = result.append
sqrt_n = (int(n ** .5) + 1) | 1    # ensure it's odd
for p in range(3, sqrt_n, 2):
if prime[p]:
append(p)
for i in range(p * p, n, 2 * p):
prime[i] = False
for p in range(sqrt_n, n, 2):
if prime[p]:
append(p)
return result


Again, this makes very little difference:

>>> test(sieve4) / t1
0.286016401170129

5. We can use Python's slice assignment instead of a loop when setting the sieve entries to False. This looks wasteful since we create a large list and then throw it away, but this avoids an expensive for loop and the associated Python interpreter overhead.

def sieve5(n):
"""Return a list of the primes below n."""
prime = [True] * n
result = [2]
append = result.append
sqrt_n = (int(n ** .5) + 1) | 1    # ensure it's odd
for p in range(3, sqrt_n, 2):
if prime[p]:
append(p)
prime[p*p::2*p] = [False] * ((n - p*p - 1) // (2*p) + 1)
for p in range(sqrt_n, n, 2):
if prime[p]:
append(p)
return result


This gives a small but noticeable improvement:

>>> test(sieve5) / t1
0.2617646381557855

6. For big improvements to the performance of numerical code, we can use NumPy.

import numpy

def sieve6(n):
"""Return an array of the primes below n."""
prime = numpy.ones(n, dtype=numpy.bool)
prime[:2] = False
prime[4::2] = False
sqrt_n = int(n ** .5) + 1
for p in range(3, sqrt_n, 2):
if prime[p]:
prime[p*p::2*p] = False
return prime.nonzero()[0]


This is more than 6 times as fast as sieve5, and more than 25 times as fast as your original code:

>>> test(sieve6) / t1
0.03726392181902129

7. We could avoid allocating space for the even numbers, improving memory locality:

def sieve7(n):
"""Return an array of the primes below n."""
prime = numpy.ones(n // 2, dtype=numpy.bool)
sqrt_n = int(n ** .5) + 1
for p in range(3, sqrt_n, 2):
if prime[p // 2]:
prime[p*p // 2::p] = False
result = 2 * prime.nonzero()[0] + 1
result[0] = 2
return result

>>> test(sieve7) / t1
0.029220096670965198

8. And finally, an implementation that sieves separately for primes of the form $6i − 1$ and $6i + 1$, due to Robert William Hanks:

def sieve8(n):
"""Return an array of the primes below n."""
prime = numpy.ones(n//3 + (n%6==2), dtype=numpy.bool)
for i in range(3, int(n**.5) + 1, 3):
if prime[i // 3]:
p = (i + 1) | 1
prime[       p*p//3     ::2*p] = False
prime[p*(p-2*(i&1)+4)//3::2*p] = False
result = (3 * prime.nonzero()[0] + 1) | 1
result[0] = 3
return numpy.r_[2,result]


This is about 40 times as fast as the original sieve:

>>> test(sieve8) / t1
0.023447068662434022

-
I'm curious, if n ** 0.5 is faster than math.sqrt(n), why isn't sqrt() implemented using n ** 0.5? Does math.sqrt(n) have better accuracy or something? – 200_success Feb 21 '14 at 20:59
@200_success: sqrt(n) involves a function call and n ** .5 doesn't. – Gareth Rees Feb 21 '14 at 21:00
One function call, not in a loop, isn't really worth fussing about. – 200_success Feb 21 '14 at 21:05
7 can be generalized using wheel factorization, which gets you another large improvement. Filling the sieve in chunks that fit into L1 or L2 cache gets another factor 5 or more according to my experience – Niklas B. Feb 22 '14 at 0:14
@NiklasB.: Sounds interesting! Can you post your code? – Gareth Rees Feb 22 '14 at 17:09

Your first setting of all the even numbers to 0 is not very efficient, the whole point of sieving is to avoid those costly modulo operations. Try the following:

l = range(2, upperlimit+1) # use list(range(...)) in Python 3
l[2::2] = [0] * ((len(l) - 3) // 2 + 1)


You can do a similar thing for the setting of zeros of the sieve for other prime numbers, but it gets complicated to figure out how many zeros to add on the right.

-
Oh nice catch, yea that saved me a decent chunk of time when generating the list – Igglyboo Feb 21 '14 at 16:27
Note that this fails in python3 (and the question is tagged python3). – Bakuriu Feb 21 '14 at 17:22
Yep, Python 3 requires wrapping the range in a list for it to work, added it as a comment above. – Jaime Feb 21 '14 at 17:59

You have correctly implemented the optimization p ** 2 > upperlimit, but, there is a more efficient way to do it...

p ** 2 > upperlimit calculates the square of p for every iteration through the loop. The product is not hugely expensive, but it is also totally unnecessary....

You can calculate the square-root of the upperlimit just once, and then reuse that calculated value as a direct comparison for p. Consider:

rootlimit = math.sqrt(upperlimit)
for p in l:
if p > rootlimit:
break;
....


Additionally, it is a small thing, but if you have break, continue, or return, inside a conditional inside a loop, then there is no need to use elif... It is just a small thing, but your code could be:

rootlimit = math.sqrt(upperlimit)
for p in l:
if p > rootlimit:
break;
if p:
....

-
for p in l:
if p ** 2 > upperlimit:
break
elif p:
...


In this loop p is often zero, while p ** 2 > upperlimit is always false until the break. That means you evaluate both conditions always. If you put if p: first, only one condition is evaluated whenever p is not prime.

-