# Finding as many perfect numbers as possible

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()

• You could speed up the is_prime by building a sieve. – Ludisposed Dec 5 '17 at 21:26
• I looked at sieves but didn't really understand them - would you mind linking to a good example or explaining a bit – 13ros27 Dec 5 '17 at 21:27
• 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. – Oscar Smith Dec 5 '17 at 21:40
• 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? – Sᴀᴍ Onᴇᴌᴀ Dec 6 '17 at 17:48
• 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. ;) – Quuxplusone Dec 25 '17 at 6:10

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.

• 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 – Ludisposed Dec 5 '17 at 22:00
• I'm pretty sure that it's a necessary but not sufficient condition. – Oscar Smith Dec 5 '17 at 22:03
• 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 Dec 6 '17 at 7:54
• It gets stuck on after the fourth perfect number saying that postponed_sieve() doesn't exist – 13ros27 Dec 6 '17 at 8:45
• It needs python 3 – Oscar Smith Dec 6 '17 at 12:52

# 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 = sieve = 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
• 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 Dec 5 '17 at 22:33
• As a result, it is exponentially slower than op, and memory errors before producing results that are almost instant in the op. – Oscar Smith Dec 6 '17 at 1:26
• The reason I print them is that I want my code to just keep running until you kill it – 13ros27 Dec 6 '17 at 9:00
• @13ros27 With yield you could achieve the same results – Ludisposed Dec 6 '17 at 9:10

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)