LUCAS-LEHMER DETAILS

 

This program uses the Lucas-Lehmer primality test to see if 2P-1 is prime.

 

The Lucas sequence is defined as:

                S (1) = 4

                S (n+1) = (S (n)2 - 2) mod (2P-1)

 

2P-1 is prime if and only if S (p-1) = 0.

 

This program uses a discrete weighted transform (see Mathematics of Computation, January 1994) to square numbers in the Lucas-Lehmer sequence.

 

The Lucas-Lehmer primality test is remarkably simple. For example, to prove 27 - 1 is prime:

 
                                         S (1) = 4
                                         S (2) = (4 * 4 – 2) mod 127 = 14
                                         S (3) = (14 * 14 – 2) mod 127 = 67
                                         S (4) = (67 * 67 – 2) mod 127 = 42
                                         S (5) = (42 * 42 – 2) mod 127 = 111
                                         S (6) = (111 * 111 – 2) mod 127 = 0

 

To implement the Lucas-Lehmer test efficiently, one must find the fastest way to square huge numbers modulo 2P-1. Since the late 1960's the fastest algorithm for squaring large numbers is to split the large number into pieces forming a large array, then perform a Fast Fourier Transform (FFT), a squaring, and an Inverse Fast Fourier Transform (IFFT). See the "How Fast Can We Multiply?" section in Knuth's Art of Computer Programming vol. 2. In a January 1994 Mathematics of Computation article by Richard Crandall and Barry Fagin titled "Discrete Weighted Transforms and Large-Integer Arithmetic", the concept of using an irrational base FFT was introduced. This improvement more than doubled the speed of the squaring by allowing us to use a smaller FFT and it performs the mod 2P-1 step for free. Although GIMPS uses a floating point FFT for reasons specific to the Intel Pentium architecture, Peter Montgomery showed that an all-integer weighted transform can also be used.

 

As mentioned in the last paragraph, GIMPS uses floating point FFTs written in highly pipelined, cache friendly assembly language. Since floating-point computations are inexact, after every iteration the floating point values are rounded back to integers. The discrepancy between the proper integer result and the actual floating-point result is called the convolution error. If the convolution error ever exceeds 0.5 then the rounding step will produce incorrect results - meaning a larger FFT should have been used. One of GIMPS' error checks is to make sure the maximum convolution error does not exceed 0.4. Unfortunately, this error check is fairly expensive and is not done on every squaring. There is another error check that is fairly cheap. One property of FFT squaring is that:

 
                    (Sum of the input FFT values) 2 = (sum of the output IFFT values)

 

Since we are using floating-point numbers we must change the "equals sign" above to "approximately equals". If the two values differ by a substantial amount, then you get a SUMINP! = SUMOUT error as described in the readme.txt file. If the sum of the input FFT values is an illegal floating point value such as infinity, then you get an ILLEGAL SUMOUT error. Unfortunately, these error checks cannot catch all errors, which brings us to our next section.

 

What are the chances that the Lucas-Lehmer test will find a new Mersenne prime number? A simple approach is to repeatedly apply the observation that the chance of finding a factor between 2X and 2X+1 is about 1/x. For example, you are testing 2 10000139-1 which trial factoring has proved there are no factors less than 264. The chance that it is prime is the chance of no 65-bit factor * chance of no 66-bit factor * ... * chance of no 5000070-bit factor. That is:

 
                    64   65         5000069
                    -- * -- * ... * -------
                    65   66         5000070

 

This simplifies to 64 / 5000070 or 1 in 78126. The general formula is “how far factored” divided by (exponent divided by 2). There are more rigorous and accurate analyses, but even these are unproven.