Introduction

A Mersenne prime is a prime number of the form M_p = 2^p – 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, M₂ = 3, M₃ = 7,
M₅ = 31 and M₇ = 127, are prime.

Mersenne claimed that 2^p – 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 M₆₇ is prime was a typographical error of M₆₁. 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 M₂₅₇, 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 M₃₁ is prime sometime before 1772. He did this by showing that any prime divisor of M₃₁ 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 M_p, 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 M_p. The functions below use this theorem to provide a list of possible factors of M_p that are less than .

We use these functions to quickly find a factor of 2⁴¹ – 1. Note that q is a factor of 2^p – 1 if and only if 2^p ≡ 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 M₁₂₇, 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 < 2^p – 1.

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

The function below encodes this method to compute k (mod 2^p – 1) using bit operations only (no division). Notice that 2ⁿ – 1 has the binary form 111 … 111₂, 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 s₀ = 4, s²_i – 1 – 2; i > 0. Then M_p = 2^p – 1 is prime iff s_p – 1 ≡ 0(modM_p).

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

To efficiently test if 2^p – 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 2^p – 1.

M₁₃ to M₂₀

The first Mersenne prime discovered by a computer running the Lucas–Lehmer test was M₅₂₁, 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.

The 20^th 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, M₄₄₂₃, 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 = 2^p – 1(2^p – 1) = 2^p – 1M_p. The Wolfram Language function PerfectNumberQ checks if a number is perfect. We demonstrate this property directly for M₃₁.

The 21^st, 22^nd and 23^rd Mersenne Primes

We proceed to rediscover #21 = M₉₆₈₉, #22 = M₉₉₄₁ and #23 = M₁₁₂₁₃. 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 2^p – 1 for primes 7,927 ≤ p ≤ 17,389.

The 24^th, 25^th and 26^th Mersenne Primes

We next extend the search to find #24 = M₁₉₉₃₇,
#25 = M₂₁₇₀₁ and #26 = M₂₃₂₀₉. The last of these was discovered in February of 1979 by Landon Curt Noll and Laura Nickel. They searched the range M₂₁₀₀₁ to M₂₄₄₉₉ 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 27^th and 28^th Mersenne Primes

We next test the range 27457 ≤ p ≤ 48611 to locate
#27 = M₄₄₄₉₇. 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 = M₈₆₂₄₃. 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.

The number M₈₆₂₄₃ 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 29^th 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 = M₁₁₀₅₀₃, 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 30^th and 31^st Mersenne Primes

The next two Mersenne primes, M₁₃₂₀₄₉, the 30^th, and M₂₁₆₀₉₁, the 31^st, 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 M_p in this range.

The Great Internet Mersenne Prime Search

The discovery of the 34^th Mersenne prime, M_1257787, 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.

We verify the 34^th 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 M_p 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
2^p – 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 2^p – 1 only if 2^p ≡ 1 (mod q), which can be quickly checked with PowerMod. At the time of this writing, all Mersenne numbers up to 2¹²⁰¹ – 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 2¹¹⁹⁹ – 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 2¹²⁰¹ – 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.

• GIMPS
• The Prime Pages
• Prime Numbers

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