Infinite Prime Number Generator In Python
I was looking at this article (The Genuine Sieve of Eratosthenes / PDF) about a common functional prime number generator that is mistakenly called the Sieve of Eratosthenes.
It is quite complicated, section 2 deals with the performance of the naive algorithm and 3 implements the sieve properly. I could not quite grok the Haskell but the following quote helped me know what was going on
Whereas the original algorithm crosses off all multiples of a prime at once, we perform these “crossings off” in a lazier way: crossing off just-in-time. For this purpose, we will store a table in which, for each prime p that we have discovered so far, there is an “iterator” holding the next multiple of p to cross off. Thus, instead of crossing off all the multiples of, say, 17, at once (impossible, since there are infinitely many for our limit-free algorithm), we will store the first one (at 17 × 17; i.e., 289) in our table of upcoming composite numbers. When we come to consider whether 289 is prime, we will check our composites table and discover that it is a known composite with 17 as a factor, remove 289 from the table, and insert 306 (i.e., 289+17). In essence, we are storing “iterators” in a table keyed by the current value of each iterator.
Here it is using a dictionary for each composite number with a list of increments it was reached by:
from itertools import count
def gen_primes():
yield 2
def update_composites(number,increment):
try:
composites[number] += [increment]
except:
composites[number] = [increment]
composites = {4:[2]}
c = count(3)
while 1:
#print composites
next = c.next()
smallest = min(composites.keys())
incs = composites[smallest]
del composites[smallest]
for inc in incs:
update_composites(smallest+inc,inc)
if next < smallest:
yield next
c.next()
update_composites(next**2,next)
g=gen_primes()
for i in range(100):
print g.next()
(On Github)
They mention a speedup moving to a heapqueue as the data structure (as we just take the lowest each lookup) so I implemented that too, here is the code
from itertools import count
from heapq import heappush,heappop
def gen_primes():
yield 2
def update_composites(number,increment):
heappush(composites,(number,increment))
composites = [(4,2)]
c = count(3)
while 1:
#print composites
smallest,increment = heappop(composites)
update_composites(smallest+increment,increment)
if composites[0][0] == smallest:
continue
next = c.next()
if next < smallest:
yield next
c.next()
update_composites(next**2,next)
g=gen_primes()
for i in range(100):
print g.next()
(On github)