3

I need to generate prime numbers using generator in Python. Here is my code:

def genPrimes():
    yield 2
    x=2
    while True:
        x+=1
        for p in genPrimes():
            if (x%p)==0:
                break
        else:
            yield x

I have a RuntimeError: maximum recursion depth exceeded after the 2nd prime.next() when I run it.

MattDMo
  • 90,104
  • 20
  • 213
  • 210
user1347096
  • 123
  • 2
  • 4
  • 10

7 Answers7

8

genPrimes() unconditionally calls itself with exactly the same arguments. This leads to infinite recursion.

Here is one way to do it using a (non-recursive) generator:

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

Note that this is optimized for simplicity and clarity rather than performance.

NPE
  • 438,426
  • 93
  • 887
  • 970
  • thank you. what would do primes=gen_primes() compare to primes=set()? – user1347096 Mar 29 '13 at 17:28
  • is `set` guaranteed to produce its contents to the `for` enumeration in the order of their addition into the set? – Will Ness Apr 01 '13 at 12:18
  • @WillNess: I don't think it is, but I don't see how that matters. – NPE Apr 01 '13 at 13:13
  • thank you. It only matters for efficiency, but as you said, you're not concerned with it here, so here it doesn't matter. – Will Ness Apr 01 '13 at 13:15
  • @WillNess: That's a good point. I've tested this function with both `list` and `set` using CPython 2.7.3, and there was no measurable difference (probably due to the properties of the hash function CPython uses for integers). – NPE Apr 01 '13 at 13:20
  • 1
    on this version it won't matter of course as it has a ***much*** bigger problem, not stopping at the square root of the number being tested. And apparently, it *can't* stop, as the primes come in *unordered*. :) – Will Ness Apr 01 '13 at 13:21
8

The fastest way to generate primes is with a sieve. Here we use a segmented Sieve of Eratosthenes to generate the primes, one by one with no maximum, in order; ps is the list of sieving primes less than the current maximum and qs is the offset of the smallest multiple of the corresponding ps in the current segment.

def genPrimes():
    def isPrime(n):
        if n % 2 == 0: return n == 2
        d = 3
        while d * d <= n:
            if n % d == 0: return False
            d += 2
        return True
    def init(): # change to Sieve of Eratosthenes
        ps, qs, sieve = [], [], [True] * 50000
        p, m = 3, 0
        while p * p <= 100000:
            if isPrime(p):
                ps.insert(0, p)
                qs.insert(0, p + (p-1) / 2)
                m += 1
            p += 2
        for i in xrange(m):
            for j in xrange(qs[i], 50000, ps[i]):
                sieve[j] = False
        return m, ps, qs, sieve
    def advance(m, ps, qs, sieve, bottom):
        for i in xrange(50000): sieve[i] = True
        for i in xrange(m):
            qs[i] = (qs[i] - 50000) % ps[i]
        p = ps[0] + 2
        while p * p <= bottom + 100000:
            if isPrime(p):
                ps.insert(0, p)
                qs.insert(0, (p*p - bottom - 1)/2)
                m += 1
            p += 2
        for i in xrange(m):
            for j in xrange(qs[i], 50000, ps[i]):
                sieve[j] = False
        return m, ps, qs, sieve
    m, ps, qs, sieve = init()
    bottom, i = 0, 1
    yield 2
    while True:
        if i == 50000:
            bottom = bottom + 100000
            m, ps, qs, sieve = advance(m, ps, qs, sieve, bottom)
            i = 0
        elif sieve[i]:
            yield bottom + i + i + 1
            i += 1
        else: i += 1

A simple isPrime using trial division is sufficient, since it will be limited to the fourth root of n. The segment size 2 * delta is arbitrarily set to 100000. This method requires O(sqrt n) space for the sieving primes plus constant space for the sieve.

It is slower but saves space to generate candidate primes with a wheel and test the candidates for primality with an isPrime based on strong pseudoprime tests to bases 2, 7, and 61; this is valid to 2^32.

def genPrimes(): # valid to 2^32
    def isPrime(n):
        def isSpsp(n, a):
            d, s = n-1, 0
            while d % 2 == 0:
                d /= 2; s += 1
            t = pow(a,d,n)
            if t == 1: return True
            while s > 0:
                if t == n-1: return True
                t = (t*t) % n; s -= 1
            return False
        for p in [2, 7, 61]:
            if n % p == 0: return n == p
            if not isSpsp(n, p): return False
        return True
    w, wheel = 0, [1,2,2,4,2,4,2,4,6,2,6,4,2,4,\
        6,6,2,6,4,2,6,4,6,8,4,2,4,2,4,8,6,4,6,\
        2,4,6,2,6,6,4,2,4,6,2,6,4,2,4,2,10,2,10]
    p = 2; yield p
    while True:
        p = p + wheel[w]
        w = 4 if w == 51 else w + 1
        if isPrime(p): yield p

If you're interested in programming with prime numbers, I modestly recommend this essay at my blog.

user448810
  • 16,364
  • 2
  • 31
  • 53
  • what are other short lists of bases for `isSpsp` and their corresponding ranges of validity? Where can one find those? Thanks. – Will Ness May 12 '13 at 14:14
  • @WillNess: The "best solution" is the smallest number that fools the test. For instance, given the three-prime set 2, 7, 61, the smallest composite number that the test reports as prime is 4759123141. Or else it's the largest number that doesn't fool the test. I forget which, but it would be easy for you to check. Charles Greathouse and Jeff Gilchrist have also done work in this area, if you want to google for their websites, but if you just want the answer, the page I linked is the simplest place to look. – user448810 May 12 '13 at 15:39
  • thanks! also found the simple list for "small" numbers in the WP page linked to from that page you mentioned. :) – Will Ness May 12 '13 at 15:46
2

A good, fast way to find primes. n is the upper limit to stop searching.

def prime(i, primes):
    for prime in primes:
        if not (i == prime or i % prime):
            return False
    primes.add(i)
    return i

def find_primes(n):
    primes = set([2])
    i, p = 2, 0
    while True:
        if prime(i, primes):
            p += 1
            if p == n:
                return primes
        i += 1
MattDMo
  • 90,104
  • 20
  • 213
  • 210
2

Here's a fast and clear imperative prime generator not using sieves:

def gen_primes():
    n = 2
    primes = []
    while True:
        is_prime = True
        for p in primes:
            if p*p > n:
                break
            if n % p == 0:
                is_prime = False
                break
        if is_prime:
            primes.append(n)
            yield n
        n += 1

It is almost identical to NPE's answer but includes a root test which significantly speeds up the generation of large primes. The upshot is the O(n) space usage for the primes list.

Rabih Kodeih
  • 8,540
  • 10
  • 44
  • 54
1

Very nice! You just forgot to stop taking primes from the inner genPrimes() when the square root of x is reached. That's all.

def genPrimes():
    yield 2 
    x=2
    while True:
        x+=1
        for p in genPrimes():
            if p*p > x:        # 
                yield x        #
                break          # 
            if (x%p)==0:
                break
        # else:
        #    yield x

Without it, you slide off the deep end, or what's the expression. When a snake eats its own tail, it must stop in the middle. If it reaches its head, there's no more snake (well, the direction of eating here is opposite, but it still applies...).

Will Ness
  • 62,652
  • 8
  • 86
  • 167
1

Just a bit more concise:

import itertools

def check_prime(n, primes):
    for p in primes:
        if not n % p:
            return False
    return True

def prime_sieve():
    primes = set()
    for n in itertools.count(2):
        if check_prime(n, primes):
            primes.add(n)
            yield n

And you can use it like this:

g = prime_sieve()
   g.next()
=> 2
   g.next()
=> 3
   g.next()
=> 5
   g.next()
=> 7
   g.next()
=> 11
bogs
  • 2,156
  • 15
  • 22
1

Here is a script that generates list of prime number from 2 to given number

from math import sqrt
next = 2
n = int(raw_input())
y = [i for i in range(2, n+1)]
while y.index(next)+1 != len(y) and next < sqrt(n):
    y = filter(lambda x: x%next !=0 or x==next, y)
    next = y[y.index(next)+1]
print y

This is just another implementation of Sieve of Eratosthenes.

Santosh Linkha
  • 13,692
  • 17
  • 72
  • 113