Math of LLR(written by Jean Penne)

Jean writes,
The algorithm used by LLR is a generalization of the Lucas-Lehmer algorithm, used to test the primality of the Mersenne numbers. I name this algorithm the Lucas-Lehmer-Riesel algorithm because it is explained and demonstrated in an article by Hans Riesel : "Lucasian Criteria for the Primality of N=h*2^n-1" Mathematics of Computation, Vol. 23 #108, pp. 869-875, Oct. 1969

The main theoretical fact is contained in the Theorem 5
(Lucas'Criteria for h*2^n-1) :

Suppose that n=2, h is odd
A =( (a+b*sqrt(D))^2)/r, Jacobi(D,N) = -1, and Jacobi(r,N)*sign(a^2-b^2*D) = -1. Then, a necessary and sufficient condition that N shall be prime is that

u(n-2) == 0 (mod. N)

if u(n) = u^2(n-1) - 2 with u(0) = A^h + A^-h.

How to use that ?



The number u(0) can be computed using a well known recursion formula :

v(0) = 2, v(1) = A+A^-1, v(k) = v(1)*v(k-1) - v(k-2).
so, we obtain u(0) = v(h) .

The remaining problem is to found a value for v(1) .

The numbers A and A^-1 are units of the quadratic field K(sqrt(D))
(that is to say units of the ring of the integers of this field...).
So, they are powers of the fundamental unit of the field.
Instead of choosing a square free integer D and searching for units
satisfying the conditions of theorem 5, Riesel takes increasing values
for v(1), and, remarking that A and A^-1 are the roots of the
equation :
A^2 - v(1)*A + 1 = 0

computes D as the square free part of v^2(1) - 4.

It remains to verify that the resulting D, a, b and r values
satisfy the conditions of theorem 5. The value of v(1)
so found is the smallest possible.

Regards,
Jean.


Chris Nash on the subject

MATH

[T.Rex]

Hi,
In the previous explaination of the LLR's Math, it is unclear what a, b, and D are.
Also, I was unable to find a version of the article of Riesel on the Web.
Can someone provide me with an electronic or paper version ?
I'm also interested in any paper dealing with "generalization of the Lucas Lehmer test", like the paper "John H. Jaroma, Austin College : A Generalization of the Lucas-Lehmer Test", presented at the "INTEGERS CONFERENCE 2005".
Can you help ?
Thanks,
Tony

The EFF's resource on prime numbers lists more detail about the test.
Here is a snippet:

A Mersenne number 2n-1 which is prime is called a Mersenne prime. If m divides n, then 2m-1 divides 2n-1, so a Mersenne prime has a prime exponent. However, very few of the numbers of the form 2p-1 (p prime) are prime. Mersenne Numbers are the easiest type of number to prove prime (because of the Lucas-Lehmer test), so are usually the largest primes on the list of largest known primes).
Primes of this form were first studied by Euclid who explored their relationship with the even perfect numbers. They were named after Mersenne because he wrote to so many mathematicians encouraging their study and because he sparked the interest of generations of mathematicians by claiming in 1644 that Mp was prime for 2, 3, 5, 7, 13, 17, 19, 31, 67, 127, 257 and for no other primes p less than 257. It took three centuries to completely test his bold claim, and when done, it was discovered that he was wrong about M67 and M257 being prime, and he omitted M61, M89, and M107. See the entry on Mersenne's conjecture for more information.


/*
* lucas - perform a Lucas primality test on h*2^n-1
*
* Copyright (C) 1999 Landon Curt Noll
*
* Calc is open software; you can redistribute it and/or modify it under
* the terms of the version 2.1 of the GNU Lesser General Public License
* as published by the Free Software Foundation.
*
* Calc is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General
* Public License for more details.
*
* A copy of version 2.1 of the GNU Lesser General Public License is
* distributed with calc under the filename COPYING-LGPL. You should have
* received a copy with calc; if not, write to Free Software Foundation, Inc.
* 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA.
*
* @(#) $Revision: 29.6 $
* @(#) $Id: lucas.cal,v 29.6 2002/07/10 09:43:46 chongo Exp $
* @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/lucas.cal,v $
*
* Under source code control: 1990/05/03 16:49:51
* File existed as early as: 1990
*
* chongo <was here> /\oo/\ http://www.isthe.com/chongo/
* Share and enjoy! :-) http://www.isthe.com/chongo/tech/comp/calc/
*/

/*
* NOTE: This is a standard calc resource file. For information on calc see:
*
* http://www.isthe.com/chongo/tech/comp/calc/index.html
*
* To obtain your own copy of calc, see:
*
* http://www.isthe.com/chongo/tech/com...-download.html
*/

/*
* HISTORICAL NOTE:
*
* On 6 August 1989 at 00:53 PDT, the 'Amdahl 6', a team consisting of
* John Brown, Landon Curt Noll, Bodo Parady, Gene Smith, Joel Smith and
* Sergio Zarantonello proved the following 65087 digit number to be prime:
*
* 216193
* 391581 * 2 -1
*
* At the time of discovery, this number was the largest known prime.
* The primality was demonstrated by a program implementing the test
* found in these routines. An Amdahl 1200 takes 1987 seconds to test
* the primality of this number. A Cray 2 took several hours to
* confirm this prime. As of 31 Dec 1995, this prime was the 3rd
* largest known prime and the largest known non-Mersenne prime.
*
* The same team also discovered the following twin prime pair:
*
* 11235 11235
* 1706595 * 2 -1 1706595 * 2 +1
*
* At the time of discovery, this was the largest known twin prime pair.
*
* See:
*
* http://www.isthe.com/chongo/tech/mat...e/amdahl6.html
*
* for more information on the Amdahl 6 group.
*
* NOTE: Both largest known and largest known twin prime records have been
* broken. Rather than update this file each time, I'll just
* congratulate the finders and encourage others to try for
* larger finds. Records were made to be broken afterall!
*/

/* ON GAINING A WORLD RECORD:
*
* The routines in calc were designed to be portable, and to work on
* numbers of 'sane' size. The Amdahl 6 team used a 'ultra-high speed
* multi-precision' package that a machine dependent collection of routines
* tuned for a long trace vector processor to work with very large numbers.
* The heart of the package was a multiplication and square routine that
* was based on the PFA Fast Fourier Transform and on Winograd's radix FFTs.
*
* Having a fast computer, and a good multi-precision package are
* critical, but one also needs to know where to look in order to have
* a good chance at a record. Knowing what to test is beyond the scope
* of this routine. However the following observations are noted:
*
* test numbers of the form h*2^n-1
* fix a value of n and vary the value h
* n mod 2^x == 0 for some value of x, say > 7 or more
* h*2^n-1 is not divisible by any small prime < 2^40
* 0 < h < 2^39
* h*2^n+1 is not divisible by any small prime < 2^40
*
* The Mersenne test for '2^n-1' is the fastest known primality test
* for a given large numbers. However, it is faster to search for
* primes of the form 'h*2^n-1'. When n is around 200000, one can find
* a prime of the form 'h*2^n-1' in about 1/2 the time.
*
* Critical to understanding why 'h*2^n-1' is to observe that primes of
* the form '2^n-1' seem to bunch around "islands". Such "islands"
* seem to be getting fewer and farther in-between, forcing the time
* for each test to grow longer and longer (worse then O(n^2 log n)).
* On the other hand, when one tests 'h*2^n-1', fixes 'n' and varies
* 'h', the time to test each number remains relatively constant.
*
* It is clearly a win to eliminate potential test candidates by
* rejecting numbers that that are divisible by 'small' primes. We
* (the "Amdahl 6") rejected all numbers that were divisible by primes
* less than '2^40'. We stopped looking for small factors at '2^40'
* when the rate of candidates being eliminated was slowed down to
* just a trickle.
*
* The 'n mod 128 == 0' restriction allows one to test for divisibility
* of small primes more quickly. To test of 'q' is a factor of 'k*2^n-1',
* one check to see if 'k*2^n mod q' == 1, which is the same a checking
* if 'h*(2^n mod q) mod q' == 1. One can compute '2^n mod q' by making
* use of the following:
*
* if
* y = 2^x mod q
* then
* 2^(2x) mod q == y^2 mod q 0 bit
* 2^(2x+1) mod q == 2*y^2 mod q 1 bit
*
* The choice of which expression depends on the binary pattern of 'n'.
* Since '1' bits require an extra step (multiply by 2), one should
* select value of 'n' that contain mostly '0' bits. The restriction
* of 'n mod 128 == 0' ensures that the bottom 7 bits of 'n' are 0.
*
* By limiting 'h' to '2^39' and eliminating all values divisible by
* small primes < twice the 'h' limit (2^40), one knows that all
* remaining candidates are relatively prime. Thus, when a candidate
* is proven to be composite (not prime) by the big test, one knows
* that the factors for that number (whatever they may be) will not
* be the factors of another candidate.
*
* Finally, one should eliminate all values of 'h*2^n-1' where
* 'h*2^n+1' is divisible by a small primes. The ideas behind this
* point is beyond the scope of this program.
*/


global pprod256; /* product of "primes up to 256" / "primes up to 46" */

/*
* lucas - lucas primality test on h*2^n-1
*
* ABOUT THE TEST:
*
* This routine will perform a primality test on h*2^n-1 based on
* the mathematics of Lucas, Lehmer and Riesel. One should read
* the following article:
*
* Ref1:
* "Lucasian Criteria for the Primality of N=h*2^n-1", by Hans Riesel,
* Mathematics of Computation, Vol 23 #108, pp. 869-875, Oct 1969
*
* The following book is also useful:
*
* Ref2:
* "Prime numbers and Computer Methods for Factorization", by Hans Riesel,
* Birkhauser, 1985, pp 131-134, 278-285, 438-444
*
* A few useful Legendre identities may be found in:
*
* Ref3:
* "Introduction to Analytic Number Theory", by Tom A. Apostol,
* Springer-Verlag, 1984, p 188.
*
* This test is performed as follows: (see Ref1, Theorem 5)
*
* a) generate u(0) (see the function gen_u0() below)
*
* b) generate u(n-2) according to the rule:
*
* u(i+1) = u(i)^2-2 mod h*2^n-1
*
* c) h*2^n-1 is prime if and only if u(n-2) == 0 Q.E.D. :-)
*
* Now the following conditions must be true for the test to work:
*
* n >= 2
* h >= 1
* h < 2^n
* h mod 2 == 1
*
* A few misc notes:
*
* In order to reduce the number of tests, as attempt to eliminate
* any number that is divisible by a prime less than 257. Valid prime
* candidates less than 257 are declared prime as a special case.
*
* In real life, you would eliminate candidates by checking for
* divisibility by a prime much larger than 257 (perhaps as high
* as 2^39).
*
* The condition 'h mod 2 == 1' is not a problem. Say one is testing
* 'j*2^m-1', where j is even. If we note that:
*
* j mod 2^x == 0 for x>0 implies j*2^m-1 == ((j/2^x)*2^(m+x))-1,
*
* then we can let h=j/2^x and n=m+x and test 'h*2^n-1' which is the value.
* We need only consider odd values of h because we can rewrite our numbers
* do make this so.

Am I reading this correctly?

* By limiting 'h' to '2^39' and eliminating all values divisible by
* small primes < twice the 'h' limit (2^40), one knows that all
* remaining candidates are relatively prime. Thus, when a candidate
* is proven to be composite (not prime) by the big test, one knows
* that the factors for that number (whatever they may be) will not
* be the factors of another candidate.
*
* Finally, one should eliminate all values of 'h*2^n-1' where
* 'h*2^n+1' is divisible by a small primes. The ideas behind this
* point is beyond the scope of this program.
*/

I may be able to make use of this, if it is not already used by Newpgen and LLR.