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

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.

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

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)