eratosthenes_sieve.py

#!/usr/bin/env python2
''' Eratosthenes sieve
'''
from __future__ import print_function

import array
import time


start = time.time()
limit = 2*10**6                        # Four million
sieve = array.array('l',range(limit))  
primes = [2]                           # Start with just first, then append to it
n = 0                                  # Index into current prime being eliminated
    
while n*n < limit:
    p = primes[n]
    # Eliminate all multiples of primes[n]:
    for i in range(p+p,limit, p):
        sieve[i] = 0
    # Find remaining numbers up to new upperbound (current prime squared)
    for i in range(primes[-1]+1,min(limit, p*p)):  # min() to avoid out of bounds access.
        if sieve[i]:
            primes.append(sieve[i])
    n+=1

elapsed = start - time.time()
print(', '.join(['%s' % x for x in primes[:20]]), "...", 
      ', '.join(['%s' % x for x in primes[-10:]]))
print("%d primes found in %g seconds" % (len(primes), elapsed))

This can find the first 2 million primes in about 2 seconds on my laptop. It's possible to do better with vectorized operations over NumPy ndarrays; but this is a basis for comparison.

For this particular case the use of array.array over just using a plain list seems to give less than a 2% performance benefit.

A quick and dirty test with the limit set at 20 million gives me 10 seconds for the array.array version and 13 seconds for the list version. So that's about a 25% advantage at that scale.

https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes for those who want to know more about this algorithm.