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 &lt;= 1: return False
    elif num &lt;= 3: return True
    elif num % 2 == 0 or num % 3 == 0: return False
    i = 5
    while i * i &lt;= 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()
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 &lt; 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 &amp; 2^n + k&gt;&gt;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 &gt; mersenne_prime:
        n = (n &amp; mersenne_prime) + (n &gt;&gt; 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.