Prime number generator

[Ian Kelly]

On Wed, Jul 10, 2013 at 10:16 AM, Ian Kelly <ian.g.kelly at gmail.com> wrote:
> The other interesting thing about my sieve is that it's a recursive
> generator.  I'll dig it up later and share it.

As promised.  Apologies for the excessive commenting.  As noted, this
implementation is a recursive generator, which is done so that the
primes in the sieve can go only up to the square root of the current
prime, rather than tossing in every prime seen so far.

from collections import defaultdict
from itertools import count, islice

def primes():
    yield 2
    # Set up the data structures.  multiples implements the sieve as a dict
    # that marks upcoming composite numbers and records the prime factors
    # that disqualified them.
    multiples = defaultdict(list)
    # subgen is a recursive sub-generator that supplies prime factors to be fed
    # into the sieve.  islice is used to skip 2 (since we only consider odd
    # numbers) and 3 (since we can't get it from subgen before we've yielded it
    # below).
    subgen = islice(primes(), 2, None)
    # Initialize the sieve with the first prime factor q.
    q = 3
    # When we reach threshold (the square of the last prime added to the sieve),
    # it's time to add another prime to the sieve.  By construction, threshold
    # is always an odd composite.
    threshold = q * q
    multiples[threshold].append(q + q)
    # The main loop tries every odd integer, starting with 3.
    for p in count(3, 2):
        if p in multiples:
            # p is composite.  Discard the entry from multiples and mark the
            # next odd multiple of each prime factor that disqualified p.
            # Because the multiples dict actually stores the doubles of the
            # prime factors, they need only be added once to get to the next
            # odd multiple.
            for q in multiples[p]:
                multiples[p + q].append(q)
            del multiples[p]
            if p == threshold:
                # All primes up to the previous q ** 2 have been generated.
                # Get the next prime factor q and add it to the sieve.
                q = subgen.next()
                threshold = q * q
                multiples[threshold].append(q + q)
        else:
            # p is prime.  Yay!
            yield p