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 (0) = 4
S (1) = 4 * 4 - 2 mod 127 = 14
S (2) = 14 * 14 - 2 mod 127 = 67
S (3) = 67 * 67 - 2 mod 127 = 42
S (4) = 42 * 42 - 2 mod 127 = 111
S (5) = 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.