ActiveState Code

Recipe 117119: Sieve of Eratosthenes

Computes an infinite sequence of primes using simple generators. A Python dictionary is used to mark multiples of the generated primes, according to the Sieve of Eratosthenes.

Python 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
# Sieve of Eratosthenes
# David Eppstein, UC Irvine, 28 Feb 2002

from __future__ import generators

def eratosthenes():
	'''Yields the sequence of prime numbers via the Sieve of Eratosthenes.'''
	D = {}  # map composite integers to primes witnessing their compositeness
	q = 2   # first integer to test for primality
	while 1:
		if q not in D:
			yield q        # not marked composite, must be prime
			D[q*q] = [q]   # first multiple of q not already marked
		else:
			for p in D[q]: # move each witness to its next multiple
				D.setdefault(p+q,[]).append(p)
			del D[q]       # no longer need D[q], free memory
		q += 1

Discussion

As Wim Stolker suggested, I have modified this to use setdefault, 1 Mar 2001.

Comments

  1. 1. At 2:07 a.m. on 1 mar 2002, Wim Stolker said:

    Shortcut. Nice short and clear recipe, however,

    Instead of:

    if D.has_key(p+q):
    D[p+q].append(p)
else:
    D[p+q] = [p]

    I would use:

    D.setdefault(p+q,[]).append(p)

    See: Add an entry to a dictionary, unless the entry is already there.

    http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/66516

  2. 2. At 7:17 a.m. on 22 mar 2003, Alex Martelli said:

    slightly faster by keeping a dict of single numbers instead of one of lists. neat indeed, but we can be marginally faster (about 40%+ on my machine for primes up to 1,000,000 for example) by keeping in the dict a single value for each composite number -- the first prime factor that revealed to us that number is a composite. Here's the Python 2.3 code for the enhanced main loop:

    while True:
    p = D.pop(q, None)
    if p:
        x = p + q
        while x in D: x += p
        D[x] = p
    else:
        D[q*q] = q
        yield q
    q += 1

    this uses the new-in-2.3 "fetch and remove" method .pop of dictionaries, but that's marginal (and I put just the same small enhancement in the list-using version -- I posted both to comp.lang.python just now) -- the small speedup in this variant comes from not having to build so many small lists.

    Alex

  3. 3. At 1:37 p.m. on 25 oct 2005, N N said:

    I hapen to like haskell's. haskell's canonical version:

    primes = sieve [ 2.. ]
         where
         sieve (p:x) = p : sieve [ n | n <- x, n `mod` p > 0 ]

    translates roughly to python as:

    import itertools as it
def sieve(num=it.count(2)):
    p=num.next()
    x=sieve(n for n in num if n%p)
    yield p
    while True: yield x.next()

    Unfortunately it exceeds the recursion limit after about 3000 on my machine.

  4. 4. At 6:58 a.m. on 13 feb 2006, Tim Hochberg said:

    Slightly faster yet by skipping even numbers.

    def sieve():
    yield 2
    D = {}
    q = 3
    while True:
        p = D.pop(q, 0)
        if p:
            x = q + p
            while x in D: x += p
            D[x] = p
        else:
            yield q
            D[q*q] = 2*q
        q += 2

    It's too bad itertools.count doesn't take a step argument! This would look cleaner, and likely be marginally faster, with

    for q in count(3,step=2)

    rather than the current while True cruft.

  5. 5. At 7:27 a.m. on 13 feb 2006, Tim Hochberg said:

    Another recursive version.

    from itertools import count, ifilter
def sieve():
    seq = count(2)
    while True:
        p = seq.next()
        seq = ifilter(p.__rmod__, seq)
        yield p

    This uses ifilter to generate the recursive sequence. It's much less vulnerable to stack overflow than the version that recurses on sieve, but otherwise it's pretty similar.

    I managed to generate all primes up to 909,691 before it bombed the Python interpreter mysteriously, but it takes a while!

Sign in to comment