Sieve of Eratosthenes - Python

Question

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]

Answer 1

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
   

Answer 2

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.

Answer 3

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

Answer 4

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.