Infinite prime number generator

Thomas Jollans thomas at jollans.com
Tue Jun 29 12:51:24 EDT 2010


I've been toying with Haskell a bit, and after implementing
(essentially) the Sieve of Eratosthenes as an infinite list, thus:

primes = 1 : foldr elim_mult [] [2..]
    where elim_mult n l = n : filter ((/=0) . (`mod` n)) l

I wondered how easy it would be to do the same thing in Python.
Obviously, Python doesn't have Haskell's infinite lists as such, and
normally evaluates expressions eagerly instead of lazily, but it's
still, with generators, relatively simple do create infinite sequences
(which aren't, strictly speaking, sequences)

So, here's the same thing as a Python generator, in a simple imperative
style:

def primes():
    yield 1
    i = 2
    old = set()
    while True:
        for p in old:
            if (i % p) == 0:
                break
        else:
            old.add(i)
            yield i
        i += 1

Haskell:
*Main> take 10 primes
[1,2,3,5,7,11,13,17,19,23]
*Main>

Python:
>>> from itertools import islice
>>> list(islice(primes(), 0, 10))
[1, 2, 3, 5, 7, 11, 13, 17, 19, 23]
>>>

So far, so good. But then I thought, how close to the haskell version
can a Python generator get?
foldr is like functools.reduce, except that it's right-associative. In
Haskell, it works for infinite lists because Haskell is a lazy bastard,
but it's not that difficult to express the exactly same thing with
recursive generators:

if isinstance(filter(None, []), list): # Python 2
    from itertools import ifilter as filter

def rprimes():
    def elim_mult(n):
        yield n
        for p in filter((lambda x:x%n != 0), elim_mult(n+1)): yield p
    yield 1
    for p in elim_mult(2): yield p

This sort of works:
>>> list(islice(rprimes(), 0, 10))
[1, 2, 3, 5, 7, 11, 13, 17, 19, 23]
>>>

But it has a bit of a problem:
>>> sl = islice(rprimes(), None, None, 150) #step=150
>>> next(sl)
1
>>> next(sl)
863
>>> next(sl)

 [......]

  File "primes.py", line 21, in elim_mult
    for p in filter((lambda x:x%n != 0), elim_mult(n+1)): yield p
  File "primes.py", line 21, in elim_mult
    for p in filter((lambda x:x%n != 0), elim_mult(n+1)): yield p
  File "primes.py", line 21, in elim_mult
    for p in filter((lambda x:x%n != 0), elim_mult(n+1)): yield p
RuntimeError: maximum recursion depth exceeded
>>>

If you happen to have a copy of stackless lying around, I'd expect it to
work!

-- Thomas



More information about the Python-list mailing list