In my previous article, I discovered the traditional brute force method of gnawing through a list of prime factors wasn’t going to make it much farther than the 9th Mersenne prime. Going any farther beyond the region of 10^20 would seem impossible with this simple algorithm; I needed a better one.

Version 4 of this mersenne.py takes advantage of a unique test that’s only available for this specific type of prime number — the Lucas-Lehmer test. The results were nothing short of incredible. Here’s the code, with a brand new method. Since it’s sort of long, I’ve hidden it behind ye olde “more” tag:

mersenne.py (version 4):

#!/usr/bin/env python

__metaclass__ = type
import time
from math import sqrt, log

def is_prime(n):
    """Determines if a positive number is prime or not."""
    if n == 2:
        return True
    if n % 2 == 0:
        return False
    max = int(sqrt(n))
    factor = 1
    while factor < max:
        factor += 2
        if n % factor == 0:
            return False
    return True

def lucas_lehmer_test(prime, mersenne_prime):
    """
    The actual Lucas-Lehmer test only works for odd exponents.

    http://en.wikipedia.org/wiki/Lucas-Lehmer_test_for_Mersenne_numbers
    """
    s = 4
    for x in range(prime-2):
        s = ((s*s) - 2) % mersenne_prime
    if s == 0:
        return True
    else:
        return False

def mersenne(limit):
    """
    Returns a list of Mersenne Primes up to a given limit.
    This new version uses the Lucas-Lehmer test:

    http://en.wikipedia.org/wiki/Lucas-Lehmer_test_for_Mersenne_numbers

    Because it only works for odd exponents, there is one Mersenne prime
    which would habitually be left off -- M(2) = 3.  Therefore, I'm starting
    with 3 as a given.  An alternative would be to add ten lines of code
    to catch a math problem that could probably be worked out by a third
    grader using pencil and paper.
    """
    cache = [3]
    base_limit = int(log(limit)/log(2))
    for x in xrange(2,base_limit):
        if is_prime(x):
            candidate = (2 ** x) - 1
            if lucas_lehmer_test(x, candidate):
                cache.append(candidate)
    return cache

if __name__ == '__main__':
    print "Start:", time.ctime()
    print mersenne(1e200)
    print "Finish:", time.ctime()

That’s right. This version is capable of going up to 10^200, that’s ten with 200 zeroes behind it.

What is hard for me to believe is how hardy Python is. The real beauty of this demonstration is not finding the next few primes, but watching this program reliably and accurately execute math such as 2^607 - 1. Here’s the results:

$ time ./mersenne4.py
Start: Mon May 21 19:09:57 2007
[3, 7, 31, 127, 8191, 131071, 524287, 2147483647L,
2305843009213693951L, 618970019642690137449562111L,
162259276829213363391578010288127L,
170141183460469231731687303715884105727L,
686479766013060971498190079908139321726943530014330540939446
345918554318339765605212255964066145455497729631139148085803
7121987999716643812574028291115057151L,
531137992816767098689588206552468627329593117727031923199444
138200403559860852242739162502265229285668889329486246501015
346579337652707239409519978766587351943831270835393219031728
127L]
Finish: Mon May 21 19:09:58 2007

real 0m0.68
user 0m0.62
sys 0m0.009s

That’s right — it made it up to the 14th Mersenne prime, a 183 digit monster.

And look at that speedup. The phrase “orders of magnitude” is too weak to describe it. Realize that Python just conducted reliable operations on enormous numbers to get this result — in less than a second! Version 3 took around 19 minutes to reach the 20 digit region — and this one has gone 180 orders of magnitude further along the number line, in an incredibly short time!

Ah, but there was trouble ahead. Look what happens when I move the upper search limit to 1e400, trying to hit that elusive 15th prime:

$ time ./mersenne4.py
Start: Mon May 21 19:30:42 2007
Traceback (most recent call last):
  File "./mersenne4.py", line 59, in ?
    print mersenne(1e400)
  File "./mersenne4.py", line 49, in mersenne
    base_limit = int(log(limit)/log(2))
OverflowError: math range error

Oops. I’ve hit a new limit — the maximum floating point value! There is more precision available in a Python long int than even a “double” precision float!

Version 5 of this program will remove this restriction.

Software examples on this page are licensed under the CC-GNU GPL.