Finding as many perfect numbers as possible

Question

I have tried to make a perfect number finder as efficient as possible in python. If anyone has any things I can add to make it more efficient or any problems with my code please say.

def is_prime(num):
 if num <= 1: return False
 elif num <= 3: return True
 elif num % 2 == 0 or num % 3 == 0: return False
 i = 5
 while i * i <= num:
 if num % i == 0 or num % (i + 2) == 0:
 return False
 i += 6
 return True
def is_mersenne_prime(prime, mersenne_prime):
 s = 4
 for i in range(prime - 2):
 s = (s*s - 2) % mersenne_prime
 if s == 0: return True
 else: return False
def calculate_perfects():
 print(6)
 prime = 3
 while True:
 if is_mersenne_prime(prime, 2**prime-1):
 print(2**(2*prime-1)-2**(prime-1))
 while True:
 prime += 2
 if is_prime(prime):
 break
calculate_perfects()

• Ludisposed: You could speed up the is_prime by building a sieve.
• 13ros27: I looked at sieves but didn't really understand them - would you mind linking to a good example or explaining a bit
• Oscar Smith: The basic idea behind a prime sieve is that instead of asking "Is this prime?". You rule out numbers that can't be. This is more efficient because crossing out multiples uses addition, while testing factors uses mod and mod is much slower.
• Sᴀᴍ Onᴇᴌᴀ: I've rolled back to the prior version, per site guidelines: "Do not add an improved version of the code after receiving an answer. Including revised versions of the code makes the question confusing, especially if someone later reviews the newer code." - Refer to the section What should I not do? on What should I do when someone answers my question?
• Quuxplusone: There are only 49 perfect numbers known (as of 2017), so you could just use a hard-coded list. It would only take about 250 megabytes of source text to represent those 49 decimal integers. ;)

Answer 1

The biggest improvement to this code is to use a prime generator that uses an efficient sieve. I will use the one here in my answer, as it is really fast.

from itertools import count
 # ideone.com/aVndFM
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 = prime_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

This function gives two advantages: first, it is much more efficient, but it also makes the rest of the logic easier. Instead of going through odd numbers, testing if they are prime, and if they are using them, the logic becomes simply

def calculate_perfects():
 yield 6
 primes = prime_sieve()
 for prime in primes:
 if is_mersenne_prime(prime):
 yield 2**(2*prime-1)-2**(prime-1)

This is not a lot faster because almost all of the time is spent in is_mersenne_prime, but it is cleaner and about 1% faster.

If we want actually faster performance, however, we need to look at is_mersenne_prime. Profiling reveals that 25% of the time is spent in the for i in range(prime-2) line. This is unfortunate, as there is very little to do to speed this up. However, the other 75% is in s = (s*s - 2) % mersenne_prime. While this initially appears to be a dead end, it isn't quite. It is probably obvious that the expensive operation here is the mod call, and thanks to some number theorists much smarter than me, it turns out that k % 2^n-1 is the same as k & 2^n + k>>n mod n. Since this only uses bitwise opps, it is much faster. below is an implementation.

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

if we call this in is_mersenne_prime, it is over 3x as fast as before. Here is the updated is_mersenne_prime code

def is_mersenne_prime(prime):
 mersenne_prime = 2**prime - 1
 s = 4
 for _ in range(prime - 2):
 s = mod_mersenne((s*s - 2), prime, mersenne_prime)
 return s == 0

On my computer this takes 5.2 instead of 16.4 seconds to generate the first 16 perfect numbers.

The next improvement we can get comes from using multiple processes. Each time is_mersenne_prime is run, it is run with information that doesn't depend on any other run. As such, we can test several numbers at a time. Here is the code that does this.

from itertools import count, compress
from multiprocessing import Pool
def calculate_perfects():
 yield 6
 primes = prime_sieve()
 pool = Pool(processes=8)
 while True:
 next_primes = [next(primes) for _ in range(8)]
 is_mersenne = pool.map(is_mersenne_prime, next_primes)
 for prime in compress(next_primes, is_mersenne):
 yield 2**(2*prime-1)-2**(prime-1)

This code is a little uglier, two lines longer, but can calculate the first 20 mersenne primes in 2.3 seconds (1.5 on pypy3)

More speedups can be found by not running the test if 2**prime-1 has an easily findable small factor. Such factors must take the form 2*k*prime+1, and factor in (1, 7) mod(8).factor will be in (1,7) at specific points depending on whether prime=4n+1 or 4n-1. The following code checks for these factors, and is a good first check before Lucas-Lehrer

def has_small_factor(prime, limit):
 """ Does 2**prime-1 have a factor less than 2*prime*limit? """
 step = 2 * prime
 if prime % 4 == 1:
 wheel = cycle((0,0,1,1))
 else:
 wheel = cycle((1,0,0,1))
 for factor in compress(range(1 + step, step*int(limit), step), wheel):
 if factor%15 in (3, 5, 9):
 continue
 if pow(2, prime, factor)-1 in (0, factor):
 return True
 return False

At this point the slowest thing about our code is that we have to multiply large numbers together. The good news is that gmpy2 has a library that makes this faster. Importing it and modifying the code to be mersenne_prime = 2**mpz(prime) -1, yields a 3x speedup (although it doesn't work well with pypy). At this point my laptop can find the first 24 perfect numbers in 33 seconds.

• Ludisposed: if is_mersenne_prime(prime, 2**prime-1): This is unnecessary since, if the n in 2**n-1 is prime, the number will be a mersenne prime
• Oscar Smith: I'm pretty sure that it's a necessary but not sufficient condition.
• 13ros27: Oscar Smith is correct, it has been proved that if 2^n-1 is prime then n will be prime but not the opposite. 11 is prime but 2^11-1 is not.
• 13ros27: It gets stuck on after the fourth perfect number saying that postponed_sieve() doesn't exist
• Oscar Smith: It needs python 3

Answer 2

Algorithm

• Using a sieve to speed up the is_prime

Currently you calculate the same thing over and over again for each prime. You could use a sieve to calculate all primes until a limit (for instance) 2**(n_limit)-1 that way you have to calculate the primes only once.

def simple_sieve(limit):
 sieve = [True] * limit
 sieve[0] = sieve[1] = False
 for (i, isprime) in enumerate(sieve):
 if isprime:
 yield i
 for n in range(i*i, limit, i):
 sieve[n] = False

Style

• Don't put everything on the same line

if num <= 1: return False

is worse to look at then

if num <= 1: 
 return False

• Use a if __name__ == '__main__' guard
• Don't print() variables but return or yield them

• Oscar Smith: Your code has a problem. If you run calculate_perfects(15), you only get 5 perfect numbers. This is because not all 2^p-1 are prime. It is also significantly slower, as it sieves way more prime numbers than it uses. It sieves up to 2^n_limit, but it only it would be much simpler to sieve up to n_limit and test the other primes.
• Oscar Smith: As a result, it is exponentially slower than op, and memory errors before producing results that are almost instant in the op.
• 13ros27: The reason I print them is that I want my code to just keep running until you kill it
• Ludisposed: @13ros27 With yield you could achieve the same results

Answer 3

To make this easy to keep updated I am creating a community wiki post of the latest updated code but please don't update this unless you are sure the update works and makes it faster or easier to read without making it slower. With thanks to @OscarSmith and @Ludisposed

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==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 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(6)
 primes = prime_sieve()
 next(primes) #2 is barely even a prime
 for prime in primes:
 if is_mersenne_prime(prime, 2**prime-1):
 yield(2**(2*prime-1)-2**(prime-1))
if __name__ == '__main__':
 for perfect in calculate_perfects():
 print(perfect)