Wolfram Computation Meets Knowledge

Mersenne Primes and the Lucas–Lehmer Test

Introduction

A Mersenne prime is a prime number of the form Mp = 2p – 1, where the exponent p must also be prime. These primes take their name from the French mathematician and religious scholar Marin Mersenne, who produced a list of primes of this form in the first half of the seventeenth century. It has been known since antiquity that the first four of these, M2 = 3, M3 = 7,
M5 = 31 and M7 = 127, are prime.

Marin Mersenne

Mersenne claimed that 2p – 1 was prime for primes p ≤ 257 only for
p ∈ {2,3,5,7,13,17,19,31,67,127,257}. It is easy to verify where he was correct and where he was not, using the Wolfram Language function PrimeQ. PrimeQ uses modern prime testing methods that do not require finding a factor to prove a number to be composite.


It is possible that his claim that M67 is prime was a typographical error of M61. However, it is not hard to understand why primality testing was difficult in Mersenne’s time, since trial division was one of the few tools available. For example, for M257, the smallest factor is a 15-digit number and, even with modern factoring methods, it is not easy to find. The Wolfram Language function FactorInteger uses advanced methods that enable it to factor large integers.


Euler and the Mersenne Factor Theorem

Some of the first advances in primality testing were accomplished by the great mathematician Leonhard Euler, who verified that M31 is prime sometime before 1772. He did this by showing that any prime divisor of M31 must be congruent to 1 or 62 (mod 248).


Such a relatively short list could be checked by trial division (by hand) in a reasonable amount of time in Euler’s day. His was an application of the Mersenne factor theorem, which states that if q is a divisor of Mp, then q ≡ 1 or –1 (mod 8), q ≡ 1 (mod p) and q = 2kp + 1 for some positive integer k. These facts greatly limit the possible divisors of Mp. The functions below use this theorem to provide a list of possible factors of Mp that are less than .

We use these functions to quickly find a factor of 241 – 1. Note that q is a factor of 2p – 1 if and only if 2p ≡ 1 (mod q). This enables the use of PowerMod, which provides very efficient modular exponentiation.





The following is a Mersenne number with 161,649 digits





Lucas and Lehmer

The next major advance was the discovery by Édouard Lucas of a clever method to test the primality of numbers of this form. He used his method in 1876 to verify that M127, the largest Mersenne prime discovered before the age of computers, is prime. In the early twentieth century, after the understanding of binary arithmetic and algebra became widely known, Derek Henry Lehmer refined Lucas’ method. The resulting Lucas–Lehmer primality test provides an efficient method of testing if a number of this form is prime. It does this by using the modular equivalence

This means that k is congruent to the number represented by its lowest-order p bits plus the number represented by the remaining bits. This relation can be applied recursively until k < 2p – 1.

Consider the example that follows. Here we show that for
k = 1234567891. Note that kmod-bitand-k-2-23-1 , the lowest-order 23 bits, and , the remaining bits shifted to the lowest position.






The function below encodes this method to compute k (mod 2p – 1) using bit operations only (no division). Notice that 2n – 1 has the binary form 111 … 1112, all 1s and no 0s, so it also serves as a mask for the lower-order p bits of k.

The following function encodes the Lucas–Lehmer primality test (LLT). We define the sequence s0 = 4, s2i – 1 – 2; i > 0. Then Mp = 2p – 1 is prime iff sp – 1 ≡ 0(modMp).

Note: Experiments have shown that the runtime of these functions is dominated by the large integer arithmetic.

To efficiently test if 2p – 1 is prime, it is better to first check for small prime divisors and to perform other basic primality testing. We first use the Mersenne prime divisor theorem encoded in checkMPDivisors and then the Wolfram Language function PrimeQ. If neither of these complete in a short time, then we apply the Lucas–Lehmer test.

Here we present an extended version of PrimeQ that applies the Lucas–Lehmer test for large integers of the form 2p – 1.

M13 to M20

The first Mersenne prime discovered by a computer running the Lucas–Lehmer test was M521, found by Raphael M. Robinson on January 30, 1952, using the early vacuum tube-based computer SWAC (Standards Western Automatic Computer). The Williams tube memory unit of this computer, holding 256 words of 37 bits each, is shown below.

Williams tube

The 20th Mersenne prime was discovered by Alexander Hurwitz in November of 1961 by running the Lucas–Lehmer test for about 50 minutes on an IBM 7090. We reproduce these early results below, using about 151 seconds of single-core computing time on a modern laptop.



One feature of the Wolfram Language that makes it suitable for this kind of work is its fast, large-integer arithmetic. This was a real challenge in the early days of computerized Mersenne prime searching. Researchers quickly adopted fast Fourier transform methods to convert the problem of multiplying two large integers, essentially a convolution of two lists of digits, into a simple element-by-element product of transformed digits. Fast integer multiplication is needed for the squaring step in the Lucas–Lehmer test. The Wolfram Language uses the latest platform-optimized algorithms to work with exact integer numbers with up to billions of digits. By way of example, we verify that the last of these, M4423, is indeed a Mersenne prime and show all of its digits.


Perfect Numbers

There is an interesting connection between Mersenne primes and perfect numbers. A perfect number is a number that is equal to the sum of all of its divisors (other than the number itself). Euclid suspected, and Euler finally proved, that all even, perfect numbers have the form P = 2p – 1(2p – 1) = 2p – 1Mp. The Wolfram Language function PerfectNumberQ checks if a number is perfect. We demonstrate this property directly for M31.





The 21st, 22nd and 23rd Mersenne Primes

We proceed to rediscover #21 = M9689, #22 = M9941 and #23 = M11213. These were all discovered by Donald B. Gillies running the LLT on an ILLIAC II during the spring 1963 (the article can be found here). We use nearly 6 minutes of elapsed time to test all of the numbers of the form 2p – 1 for primes 7,927 ≤ p ≤ 17,389.






The 24th, 25th and 26th Mersenne Primes

We next extend the search to find #24 = M19937,
#25 = M21701 and #26 = M23209. The last of these was discovered in February of 1979 by Landon Curt Noll and Laura Nickel. They searched the range M21001 to M24499 using 6,000 CPU hours on a CDC Cyber 174 (that article can be found here). Our computations are becoming sufficiently intense to warrant the use of parallel processing. Since the tests of the candidate prime factors are independent, we can use ParallelMap to speed up the work. We check the range
17,393 ≤ p ≤ 27,449 in about three and a half minutes using 4 cores.





Notice how the specialized Lucas–Lehmer test is significantly faster than the more general function PrimeQ for these Mersenne primes.



The 27th and 28th Mersenne Primes

We next test the range 27457 ≤ p ≤ 48611 to locate
#27 = M44497. This was discovered in April 1979 on a Cray-1 by Harry Nelson and his team. Our search of this range runs in about 15 minutes.




The next Mersenne prime is #28 = M86243. It was discovered in September of 1982 by David Slowinski, also on a Cray-1. The Cray-1 supercomputer weighed about 5 tons, consumed about 115 kilowatts of power and delivered 160 MFLOPS of computing performance. It was supplied with 1 million 64-bit words of memory (8 megabytes), and cost about $16 million in today’s dollars. A detail of its significant cooling system is shown below. By comparison, a Raspberry Pi weighs a few ounces, runs on 4 watts, delivers about 410 MFLOPS and is provided with 1 gigabyte of RAM, all for about $40, and it comes with Mathematica.

Cray-1 supercomuter

The number M86243 has 25,962 digits. In 1 hour and 14 minutes we were able to find this value (on my laptop, not on a Raspberry Pi) by testing over the range 48,619 ≤ p ≤ 87,533.





The 29th Mersenne Prime

Since we are now using serious computer time, we also produce a timestamp for each run. We now check the range 87,557 ≤ p ≤ 110,597. In 1 hour and 44 minutes, this reveals #29 = M110503, first discovered on January 29, 1988, by Walker Colquitt and Luke Welsh running the LLT on an NEC DX-2 supercomputer (the article can be found here).







The 30th and 31st Mersenne Primes

The next two Mersenne primes, M132049, the 30th, and M216091, the 31st, were actually discovered before #29, by the same team that discovered #28. They used a Cray X-MP to find #30 in September of 1983 and #31 in September 1985. We verify #30 by searching the range 110,603 ≤ p ≤ 139,901. It took nearly 4 hours and 8 minutes to check each Mp in this range.








The Great Internet Mersenne Prime Search

The discovery of the 34th Mersenne prime, M1257787, in September 1996 ended the reign of the supercomputer in the search for Mersenne primes. The next 15 were found by volunteers of the Great Internet Mersenne Prime Search (GIMPS), which runs a variant of the Lucas–Lehmer test as a background process on personal computers. This large-scale distributed computing project currently achieves a performance equivalent to approximately 300 TFLOPS per second, harnessing the otherwise idle time of more than 1.3 million computers.

GIMPS Logo

We verify the 34th Mersenne prime by directly using the Lucas–Lehmer test. We are reaching the limits of personal computer capability. Testing thousands of Mersenne numbers in this range would take many days. It is interesting to note that the Lucas–Lehmer test is often used as a stress test for the reliability of computer hardware and software, as even one arithmetic error among the billions of computations needed for testing one large prime will produce an incorrect conclusion, miss a true Mersenne prime or falsely report that a composite is prime. The fact that we have tested every Mp for primes between 2 and 139,901 is strong evidence for the reliability of large integer arithmetic and binary operations in Mathematica.


Factoring Mersenne Numbers

As we have seen, the possible factors of numbers of the form
2p – 1 are limited by the Mersenne factor theorem. This has enabled an efficient computerized search for the factors of large integers of this form. The integer q divides 2p – 1 only if 2p ≡ 1 (mod q), which can be quickly checked with PowerMod. At the time of this writing, all Mersenne numbers up to 21201 – 1 have been fully factored. A Wolfram Demonstration that illustrates these can be found here. The GIMPS project has also led to the discovery of large factors of many Mersenne numbers as a byproduct of its primality testing. A recent paper that presents a modern approach to this problem, along with 17 new factorizations, can be found here. The largest number factorized was 21199 – 1; its smallest newly found factor has 76 decimal digits. Its smallest factor is actually 23. The total computer time for this effort was nearly 7,500 core years.

We can quickly find the first few factors of 21201 – 1 using the Wolfram Language function FactorInteger.


The Wolfram Language has cataloged of all of the Mersenne primes discovered to date, with ordering up to #44. Access to this information is provided by the functions MersennePrimeExponent and MersennePrimeExponentQ.



If you find this subject interesting, you can find more details at the following websites.

Download this post as a Computable Document Format (CDF) file. New to CDF? Get your copy for free with this one-time download.

Comments

Join the discussion

!Please enter your comment (at least 5 characters).

!Please enter your name.

!Please enter a valid email address.

2 comments

  1. Thanks John

    This is a very interesting blog. But besides trying to numerically check for primality, what is the status of attempts at proving that all Mersenne numbers are prime? How does it relate to the distribution of primes and the Riemann hypothesis about the zeros of the Zeta function?

    Regards
    Michael

    Reply