A compact prime generator that runs forever
Prime sieves usually come in two flavors: fast but memory-hungry, or flexible and lazy. This article explores the latter—a compact Python generator that produces an endless sequence of primes, originally shared in a Stack Overflow answer and popularized by David Eppstein's 2002 ActiveState recipe. The core code is just a few lines:
def gen_primes():
"""Generate an infinite sequence of prime numbers."""
D = {}
q = 2
while True:
if q not in D:
D[q * q] = [q]
yield q
else:
for p in D[q]:
D.setdefault(p + q, []).append(p)
del D[q]
q += 1
The generator's useful qualities are its terseness and clarity. It does not require knowing a limit in advance, but it will consume memory indefinitely as it runs. Unlike the classic sieve, its memory footprint grows with the number of primes found—not with the range of numbers considered.
How the generator works
The algorithm relies on the dictionary D. It stores previously discovered primes as values, but rather than keying them by the prime itself, each prime is keyed by the next composite number it divides. This enables the generator to decide whether a candidate x is prime by checking whether x appears as a key in D—no division by all known primes is necessary.
When x is not in D, it is a new prime, and its square is added to D. When x is found in D, it is composite; each prime value associated with that key is moved to the next multiple, ensuring every composite eventually surfaces. The mechanism parallels the classic sieve: multiplications mark composites, and the values in D step forward in the same way the sieve's marking moves along the table.
The approach is best understood by adding print statements and tracing a few iterations. It should be noted that each prime appears in D exactly once, so its size is the value of the prime-counting function at the current number, roughly n / ln(n).
Speeding it up
The generator is compact, but the convenience costs performance—it runs over 5x slower than a classic sieve with a fixed boundary. A faster version builds on optimizations suggested by Alex Martelli, Tim Hochberg, and Wolfgang Beneicke in the original ActiveState discussion:
def gen_primes_opt():
yield 2
D = {}
for q in itertools.count(3, step=2):
p = D.pop(q, None)
if not p:
D[q * q] = q
yield q
else:
x = q + p + p # get odd multiples
while x in D:
x += p + p
D[x] = p
Three key changes make it faster:
- Replace the list stored with each composite key with a single number. When multiple primes divide the same composite, the collision is handled by searching for the next multiple of one of the witnesses, similar to linear probing in a hash table.
- Skip even numbers entirely by starting with 2 and then iterating only odd candidates.
- When advancing a witness to its next multiple, jump directly to
q + p + pto stay on odd numbers.
These adjustments bring the generator's speed to within 40% of the bounded sieve while keeping the code short. Even more elaborate approaches—such as the wheel-based optimization by Will Ness and Tim Peters, building on Sorenson's paper—offer additional speed and lower memory usage but come at the cost of readability.
For many tasks, faster execution might be more important than avoiding upfront boundary figures or exploring algorithmic refinements. For example, implementing the same generator in Go showing its experimental iterator support, even before optimizations, usually runs several times faster than the improved Python variant—the iterator and yield syntax being what makes the approach convenient in the first place.
A segmented alternative worth recalling
The Wikipedia article on the sieve of Eratosthenes also covers its segmented variant, described in section 5 of Sorenson's paper. The main insight: to sieve numbers up to N, you only need the primes up to its square root, and you can process the interval in consecutive chunks. Memory usage drops to O(√n) per segment, making the segmented sieve a practical middle ground when dealing with large ranges. A commented Python implementation shows a possible approach:
def gen_primes_upto_segmented(n):
"""Generates a sequence of primes < n.
Uses the segmented sieve or Eratosthenes algorithm with O(√n) memory.
"""
# Simplify boundary cases by hard-coding some small primes.
if n < 11:
for p in [2, 3, 5, 7]:
if p < n:
yield p
return
# We break the range [0..n) into segments of size √n
segsize = int(math.ceil(math.sqrt(n)))
# Find the primes in the first segment by calling the basic sieve on that
# segment (its memory usage will be O(√n)). We'll use these primes to
# sieve all subsequent segments.
baseprimes = list(gen_primes_upto(segsize))
for bp in baseprimes:
yield bp
for segstart in range(segsize, n, segsize):
# Create a new table of size √n for each segment; the old table
# is thrown away, so the total memory use here is √n
# seg[i] represents the number segstart+i
seg = [True] * segsize
for bp in baseprimes:
# The first multiple of bp in this segment can be calculated using
# modulo.
first_multiple = (
segstart if segstart % bp == 0 else segstart + bp - segstart % bp
)
# Mark all multiples of bp in the segment as composite.
for q in range(first_multiple, segstart + segsize, bp):
seg[q % len(seg)] = False
# Sieving is done; yield all composites in the segment (iterating only
# over the odd ones).
start = 1 if segstart % 2 == 0 else 0
for i in range(start, len(seg), 2):
if seg[i]:
if segstart + i >= n:
break
yield segstart + i
The complete source code for this post—along with tests and benchmarks—is hosted on GitHub. For discussions, email is the preferred channel.



