Naive version of the Lucas-Lehmer test for Mersenne primes

Ever since 1952, the largest known primes have all had the form 2ⁿ − 1, with one exception in 1989. The reason the largest known primes have this form is that it is easier to test whether these numbers are prime than other numbers.

A number of the form 2ⁿ − 1 is called a Mersenne number, denoted M_n. A Mersenne prime is a Mersenne number that is also prime. There is a theorem that says that if M_n is prime then n must be prime.

Lehmer’s theorem

Derrick Henry Lehmer proved in 1930 that M_pis prime if p is an odd prime and M_p divides s_p−2 where the s terms are defined recursively as follows. Define s₀ = 4 and s_n = s_n−1² − 2 for n > 0. This is the Lucas-Lehmer primality test for Mersenne numbers, the test alluded to above.

Let’s try this out on M₇ = 2⁷ − 1 = 127. We need to test whether 127 divides s₅. The first six values of s_n are

4, 14, 194, 37634, 1416317954, 2005956546822746114

and indeed 127 does divide 2005956546822746114. As you may have noticed, s₅ is a big number, and s₆ will be a lot bigger. It doesn’t look like Lehmer’s theorem is going to be very useful.

The missing piece

Here’s the idea that makes the Lucas-Lehmer [1] test practical. We don’t need to know s_p−2 per se; we only need to know whether it is divisible by M_p. So we can carry out all our calculations mod M_p. In our example, we only need to calculate the values of s_n mod 127:

4, 14, 67, 42, 111, 0

Much better.

Lucas-Lehmer test takes the remainder by M_p at every step along the way when computing the recursion for the s terms. The naive version does not, but computes the s terms directly.

To make the two versions of the test explicit, here are Python implementations.

Naive code

def lucas_lehmer_naive(p):
    M = 2**p - 1
    s = 4
    for _ in range(p-2):
        s = (s*s - 2) 
    return s % M  == 0

Good code

def lucas_lehmer(p):
    M = 2**p - 1
    s = 4
    for _ in range(p-2):
        s = (s*s - 2) % M
    return s == 0

The bad version

How bad is the naive version of the Lucas-Lehmer test?

The OEIS sequence A003010 consists of the s_n. OEIS gives the following formula:

$s_n = \left\lceil (2 + \sqrt{3})^{2^n}\right\rceil$

This shows that the s_n grow extremely rapidly, which suggests the naive algorithm will hit a wall fairly quickly

When I used the two versions of the test above to verify that the first few Mersenne primes M_p really are prime, the naive version gets stuck after p = 19. The better version immediately works through p = 9689 before slowing down noticeably.

Historical example

Let’s look at M₅₂₁. This was the first Mersenne number verified by a computer to be prime. This was done on the vacuum tube-based SWAC computer in 1952 using the smart version of the Lucas-Lehmer test.

Carrying out the Lucas-Lehmer test to show that this number is prime requires doing modular arithmetic with 521-bit numbers, which was well within the 9472 bits of memory in the vacuum tube-based SWAC computer that proved M₅₂₁ was prime.

But it would not be possible now, or ever, to verify that M₅₂₁ is prime using the naive version of the Lucas-Lehmer test because we could not store, much less compute, s₅₁₉.

Since

$\log_{2} s_{519} = 2^{519} \log_{2} (2 + \sqrt{3}) \approx 2^{520}$

we’d need around 2⁵²⁰ bits to store s₅₁₉. The number of bits we’d need is itself a 157-digit number. We’d need more bits than there are particles in the universe, which has been estimated to be on the order of 10⁸⁰.

The largest Mersenne prime that the SWAC could verify using the naive Lucas-Lehmer test would have been M₁₃= 8191, which was found to be prime in the Middle Ages.

[1] Édouard Lucas conjectured in 1878 the theorem that Lehmer proved in 1930.

The bad version of a good test