I need help regarding a Mersenne Prime Sieve which I edited to perform this task; this code is from a perfect number generator. The problem is the program slows down at 11213 in python. Can anyone offer some speed enhancements for this sieve? The output starts out like this. 2 3 5 7 13 17 19 31 61 89 107 127 521 607 1279 2203 2281 3217 4253 4423 9689 9941 11213
from itertools import count
def postponed_sieve(): # postponed sieve, by Will Ness
yield 2; yield 3; yield 5; yield 7; # original code David Eppstein,
sieve = {} # Alex Martelli, ActiveState Recipe
#2002
ps = postponed_sieve() # a separate base Primes Supply:
p = next(ps) and next(ps) # (3) a Prime to add to dict
q = p*p # (9) its sQuare
for c in count(9,2): # the Candidate
if c in sieve: # c's a multiple of some base prime
s = sieve.pop(c) # i.e. a composite ; or
elif c < q:
yield c # a prime
continue
else: #(c==q) # or the next base prime's square:
s=count(q+2*p,2*p) # (9+6, by 6 : 15,21,27,33,...)
p=next(ps) # (5)
q=p*p # (25)
for m in s: # the next multiple
if m not in sieve: # no duplicates
break
sieve[m] = s # original test entry: ideone.com/WFv4f
def prime_sieve(): # postponed sieve, by Will Ness
yield 2; yield 3; yield 5; yield 7; # original code David Eppstein,
sieve = {} # Alex Martelli, ActiveState Recipe
#2002
ps = postponed_sieve() # a separate base Primes Supply:
p = next(ps) and next(ps) # (3) a Prime to add to dict
q = p*p # (9) its sQuare
for c in count(9,2): # the Candidate
if c in sieve: # c’s a multiple of some base prime
s = sieve.pop(c) # i.e. a composite ; or
elif c < q:
yield c # a prime
continue
else: #(c==s) # or the next base prime’s square:
s=count(q+2*p,2*p) # (9+6, by 6 : 15,21,27,33,...)
p=next(ps) # (5)
q=p*p # (25)
for m in s: # the next multiple
if m not in sieve: # no duplicates
break
sieve[m] = s # original test entry: ideone.com/WFv4f
def mod_mersenne(n, prime, mersenne_prime):
while n > mersenne_prime:
n = (n & mersenne_prime) + (n >> prime)
if n == mersenne_prime:
return 0
return n
def is_mersenne_prime(prime, mersenne_prime):
s = 4
for i in range(prime - 2):
s = mod_mersenne((s*s - 2), prime, mersenne_prime)
return s == 0
def calculate_perfects():
yield(2)
primes = prime_sieve()
next(primes) #2 is barely even a prime
for prime in primes:
if is_mersenne_prime(prime, 2**prime-1):
yield(prime)
if __name__ == '__main__':
for perfect in calculate_perfects():
print(perfect)
#edited by Tom E. O'Neil to find Mprimes
