Mersenne Primes.

Introduction

Update: I recorded myself reading 419 digits from the 52nd Mersenne prime for Say the prime.

The 52nd Mersenne prime was discovered on October 21, 2024 (Source: GIMPS). Despite intense effort, it took six years to surpass the 51st! It is the largest prime number ever found so far! A prime number is a number that is only divisible by 1 and itself. A Mersenne prime is a prime number of the form \(m = 2^p - 1\), where \(p\) itself is a prime number.

The first few Mersenne primes (source: Wikipedia) are:

3, 7, 31, 127, 8191, 131071, 524287, 2147483647

The specific form of Mersenne primes enables the use of a fast prime-testing algorithm called the Lucas-Lehmer test. Mersenne primes are closely connected to perfect numbers — numbers that are equal to the sum of their proper divisors. Perfect numbers take the form \(N = 2^{p-1}(2^p-1)\), where \(2^p - 1\) is a Mersenne prime.

Mersenne primes are rare, and the numbers grow exponentially as \(p\) increases, making them increasingly difficult to discover. The computational effort needed to verify whether a number of the form \(2^p - 1\) is prime also grows significantly, requiring sophisticated algorithms and powerful hardware, especially for large candidates. In this blog post, I will generate some Mersenne primes and explore how fast they grow.

Python setup

For this project, I use the standard Python library, SymPy for fast prime number calculations, and tabulate for formatting tables in a readable way. I use the following imports:

import time 
import sympy
import tabulate

To easily print tables, I created tabuliterate, which runs a function on each value of an iterable and optionally applies a condition.

def tabuliterate(fun, iterable, cond):
    table = []
    t0 = time.perf_counter() 
    for i, j in enumerate(iterable):
        n = fun(j)
        if cond is not None and not cond(n):
            continue
        t = time.perf_counter()
        table.append([i, j, t - t0, n])
        t0 = t
    return table

Checking for prime numbers

One can check if a candidate number \(n\) is prime by verifying if it is divisible by any smaller prime up to \(\sqrt{n}\). Here’s a simple python implementation of this:

def is_prime(n):
    if n < 2:
        return False
    for i in range(2, int(n**0.5) + 1):
        if n % i == 0:
            return False
    return True

Let’s use this to print all primes in the range (1, 50).

%%time
prime_table = tabuliterate(is_prime, range(1, 50), cond=lambda x: x)
print(
    tabulate.tabulate(
        prime_table,
        headers=["i", "prime(i)", "t (s)", "is_prime"],
        tablefmt="github"))

The tabuliterate function runs is_prime for each value in the range, and saves the results if is_prime returns True. I then use tabulate to print a well-formatted, blog-ready Markdown table. To better understand performance, I use time.perf_counter to measure the time between iterations (in column t (s)). Additionally, I use the Jupyter magic %%time to print timing information for the entire cell. While these are not substitutes for proper profiling, they provide sufficient insight for this blog, optimization is not the main focus.

CPU times: user 91 μs, sys: 7 μs, total: 98 μs
Wall time: 101 μs

Now let’s try a few bigger numbers, including powers of ten and some known Mersenne primes:

%%time 
candidates = 1000, 10000, 100000, 1000000, 131071, 524287, 2147483647
is_prime_table = tabuliterate(is_prime, candidates)
print(tabulate.tabulate(
    is_prime_table,
    headers=["i", "n", "t (s)", "is_prime"],
    tablefmt="github"))

CPU times: user 8.42 ms, sys: 367 μs, total: 8.79 ms
Wall time: 8.92 ms

Using this very basic and unoptimized prime checker, we can still test prime numbers relatively quickly. However, by the time we reach the 9th Mersenne prime, the computational effort required is no longer feasible for this method without significant optimization.

Generating Prime Numbers

Before attempting to find Mersenne primes, we first need to generate prime numbers for the Mersenne exponents. Specifically, we need a function to generate the nth prime: p = prime(n). We could use the is_prime function and loop over all possible values of n, but as Mersenne primes grow very quickly, this approach will soon face performance limitations. The is_prime function doesn’t store previously checked values, leading to redundant calculations.

To overcome this, we can use a cache dictionary {n: prime(n)} to store previously computed primes, allowing us to generate the nth prime more efficiently by avoiding unnecessary recomputation. Here is an implementation to compute the n-th prime:

def nth_prime(n, cache):
    if n in cache: return cache[n]
    c, num = max(cache.keys()), cache[max(cache.keys())] + 2
    while c < n:
        is_prime = True
        limit = int(num**0.5) + 1
        for p in cache.values():
            if p > limit: break
            if num % p == 0:
                is_prime = False
                break
        if is_prime: c, cache[c + 1] = c + 1, num
        num += 2
    return cache[n]

Now let’s generate some prime numbers:

%%time
cache = {1: 2, 2: 3, 3: 5}
nth_primes_table = tabuliterate(
    lambda x: nth_prime(x, cache=cache),
    (10**i for i in range(3, 7))
)
print(
    tabulate.tabulate(
        nth_primes_table,
        headers=["i", "n", "t (s)", "prime(n)"],
        tablefmt="github"))

CPU times: user 20.6 s, sys: 124 ms, total: 20.8 s
Wall time: 20.8 s

We are able to generate the millionth prime, but the approach becomes impractical for larger primes. Notably, the 8th Mersenne prime is already 2,147,483,647, making this method unsuitable without significant optimization for larger primes.

Finding Mersenne primes

Instead of trying to reinvent a faster and smarter prime generator and tester, here is the same test implemented with SymPy:

%%time
sympy_nth_primes_table = tabuliterate(
    lambda x: sympy.prime,
    (10**i for i in range(3, 7))
)
print(tabulate.tabulate(
    sympy_nth_primes_table, 
    headers=["i", "n", "t (s)", "prime(n)"],
    tablefmt="github"))

CPU times: user 47 μs, sys: 1 μs, total: 48 μs
Wall time: 57 μs

This approach is significantly faster. Now, let’s try finding some Mersenne primes! To do this, we use sympy.isprime as a filter to determine if numbers of the form are prime:

%%time
N=600
mersenne_primes_table = tabuliterate(
    lambda x: 2**x-1,
    (sympy.prime(n) for n in range(1, N)),
    sympy.isprime
)

CPU times: user 1min 14s, sys: 574 ms, total: 1min 14s
Wall time: 1min 14s

It takes a little over one minute to test the first 600 primes. However, we can optimize this further by leveraging SymPy’s `is_mersenne_prime, which is specifically designed for testing Mersenne primes:

%%time
N=600
mersenne_primes_table = tabuliterate(
    lambda x: 2**x-1,
    (sympy.prime(n) for n in range(1, N)),
    sympy.ntheory.primetest.is_mersenne_prime
)

CPU times: user 1min 7s, sys: 94.7 ms, total: 1min 7s
Wall time: 1min 7s

We managed to test up to and found the first 28 Mersenne primes in about a minute!

Conclusion

Mersenne primes are fascinating mathematical objects with deep connections to number theory. Discovering them requires both computational ingenuity and efficient algorithms, as they grow exponentially large and exceedingly rare. I was amazed I could generate a 25962 digit Mersenne prime on my laptop in a minute. However, noticing that the last number took more than twice as much as the previous one, and the 52 is another 24 steps away… the resent discovery really is something.