Math: Lucas-Lehmer Details

# LUCAS-LEHMER DETAILS

The Lucas-Lehmer primality test is remarkably simple. It states that for P > 2, 2P-1 is prime if and only if Sp-2 is zero in this sequence: S0 = 4, SN = (SN-12 - 2) mod (2P-1). For example, to prove 27 - 1 is prime:

> S0 = 4
> S1 = (4 * 4 - 2) mod 127 = 14
> S2 = (14 * 14 - 2) mod 127 = 67
> S3 = (67 * 67 - 2) mod 127 = 42
> S4 = (42 * 42 - 2) mod 127 = 111
> S5 = (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 210000139-1 for 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. This simple approach isn't quite right. It would give a formula of how_far_factored divided by (exponent divided by 2). However, more rigorous work has shown the formula to be (how_far_factored-1) / (exponent times Euler's constant (0.577...)). In this case, 1 in 91623. Even these more rigorous formulas are unproven.