Permalink
Fetching contributors…
Cannot retrieve contributors at this time
1847 lines (1781 sloc) 48.6 KB
/*
* lucas - perform a Lucas primality test on h*2^n-1
*
* Copyright (C) 1999,2017,2018 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.
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*
* 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/
*/
/*
* For a general tutorial on how to find a new largest known prime, see:
*
* http://www.isthe.com/chongo/tech/math/prime/prime-tutorial.pdf
*
* Also see the reference code, available both C and go:
*
* https://github.com/arcetri/goprime
*/
/*
* 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/comp/calc/calc-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/math/prime/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 after all!
*/
/*
* ON GAINING A WORLD RECORD:
*
* For a general tutorial on how to find a new largest known prime, see:
*
* http://www.isthe.com/chongo/tech/math/prime/prime-tutorial.pdf
*
* 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.
*
* NOTE: While the PFA Fast Fourier Transform and Winograd's radix FFTs
* might have been optimal for the Amdahl 6 team at the time,
* they might not be optimal for your CPU architecture. See
* the above mentioned tutorial for information on better
* methods of performing multiplications and squares of very
* large numbers.
*
* 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.
*
* NOTE: Today, for world record sized h*2^n-1 primes, one might
* search for factors < 2^46 or more. By excluding h*2^n-1
* with prime factors < 2^46, where h*2^n-1 is a bit larger
* than the largest known prime, one may exclude about 96.5%
* of candidates that have "small" prime factors.
*/
pprod256 = 0; /* 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
*
* http://www.ams.org/journals/mcom/1969-23-108/S0025-5718-1969-0262163-1/
* S0025-5718-1969-0262163-1.pdf
*
* NOTE: Join the above two lines for the complete URL of the paper.
*
* 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.
*
* An excellent 5-page paper by Oystein J. Rodseth (we apologize that the
* ASCII character set does not allow us to spell his name with the
* umlaut marks on the O's):
*
* NOTE: The original Amdahl 6 method predates the publication of Ref4.
* The gen_v1() function used by lucas() uses the Ref4 method.
* See the 'Amdahl 6 legacy code' section below for the original
* method of generating v(1).
*
* Ref4:
*
* "A note on primality tests for N = h*2^n-1", by Oystein J. Rodseth,
* Department of Mathematics, University of Bergen, BIT Numerical
* Mathematics. 34 (3): pp 451-454.
*
* http://folk.uib.no/nmaoy/papers/luc.pdf
*
* This test is performed as follows: (see Ref1, Theorem 5)
*
* a) generate u(2) (see the function gen_u2() below)
* (NOTE: some call this u(0))
*
* b) generate u(n) 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) == 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 miscellaneous 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.
*
* input:
* h h as in h*2^n-1 (must be >= 1)
* n n as in h*2^n-1 (must be >= 1)
*
* returns:
* 1 => h*2^n-1 is prime
* 0 => h*2^n-1 is not prime
* -1 => a test could not be formed, or h >= 2^n, h <= 0, n <= 0
*/
define
lucas(h, n)
{
local testval; /* h*2^n-1 */
local shiftdown; /* the power of 2 that divides h */
local u; /* the u(i) sequence value */
local v1; /* the v(1) generator of u(2) */
local i; /* u sequence cycle number */
local oldh; /* pre-reduced h */
local oldn; /* pre-reduced n */
local bits; /* highbit of h*2^n-1 */
/*
* check arg types
*/
if (!isint(h)) {
quit "FATAL: bad args: h must be an integer";
}
if (h < 1) {
quit "FATAL: bad args: h must be an integer >= 1";
}
if (!isint(n)) {
quit "FATAL: bad args: n must be an integer";
}
if (n < 1) {
quit "FATAL: bad args: n must be an integer >= 1";
}
/*
* reduce h if even
*
* we will force h to be odd by moving powers of two over to 2^n
*/
oldh = h;
oldn = n;
shiftdown = fcnt(h,2); /* h % 2^shiftdown == 0, max shiftdown */
if (shiftdown > 0) {
h >>= shiftdown;
n += shiftdown;
}
/*
* enforce the 0 < h < 2^n rule
*/
if (h <= 0 || n <= 0) {
print "ERROR: reduced args violate the rule: 0 < h < 2^n";
print " ERROR: h=":oldh, "n=":oldn, "reduced h=":h, "n=":n;
ldebug("lucas", "unknown: h <= 0 || n <= 0");
return -1;
}
if (highbit(h) >= n) {
print "ERROR: reduced args violate the rule: h < 2^n";
print " ERROR: h=":oldh, "n=":oldn, "reduced h=":h, "n=":n;
ldebug("lucas", "unknown: highbit(h) >= n");
return -1;
}
/*
* catch the degenerate case of h*2^n-1 == 1
*/
if (h == 1 && n == 1) {
ldebug("lucas", "not prime: h == 1 && n == 1");
return 0; /* 1*2^1-1 == 1 is not prime */
}
/*
* catch the degenerate case of n==2
*
* n==2 and 0<h<2^n ==> 0<h<4
*
* Since h is now odd ==> h==1 or h==3
*/
if (h == 1 && n == 2) {
ldebug("lucas", "prime: h == 1 && n == 2");
return 1; /* 1*2^2-1 == 3 is prime */
}
if (h == 3 && n == 2) {
ldebug("lucas", "prime: h == 3 && n == 2");
return 1; /* 3*2^2-1 == 11 is prime */
}
/*
* catch small primes < 257
*
* We check for only a few primes because the other primes < 257
* violate the checks above.
*/
if (h == 1) {
if (n == 3 || n == 5 || n == 7) {
ldebug("lucas", "prime: 3, 7, 31, 127 are prime");
return 1; /* 3, 7, 31, 127 are prime */
}
}
if (h == 3) {
if (n == 2 || n == 3 || n == 4 || n == 6) {
ldebug("lucas", "prime: 11, 23, 47, 191 are prime");
return 1; /* 11, 23, 47, 191 are prime */
}
}
if (h == 5 && n == 4) {
ldebug("lucas", "prime: 79 is prime");
return 1; /* 79 is prime */
}
if (h == 7 && n == 5) {
ldebug("lucas", "prime: 223 is prime");
return 1; /* 223 is prime */
}
if (h == 15 && n == 4) {
ldebug("lucas", "prime: 239 is prime");
return 1; /* 239 is prime */
}
/*
* Verify that h*2^n-1 is not a multiple of 3
*
* The case for h*2^n-1 == 3 is handled above.
*/
if (((h % 3 == 1) && (n % 2 == 0)) || ((h % 3 == 2) && (n % 2 == 1))) {
/* no need to test h*2^n-1, it is a multiple of 3 */
ldebug("lucas","not-prime: != 3 and is a multiple of 3");
return 0;
}
/*
* Avoid any numbers divisible by small primes
*/
/*
* check for 5 <= prime factors < 31
* pfact(30)/6 = 1078282205
*/
testval = h*2^n - 1;
if (gcd(testval, 1078282205) > 1) {
/* a small 5 <= prime < 31 divides h*2^n-1 */
ldebug("lucas",\
"not-prime: a small 5<=prime<31 divides h*2^n-1");
return 0;
}
/*
* check for 31 <= prime factors < 53
* pfact(52)/pfact(30) = 95041567
*/
if (gcd(testval, 95041567) > 1) {
/* a small 31 <= prime < 53 divides h*2^n-1 */
ldebug("lucas","not-prime: 31<=prime<53 divides h*2^n-1");
return 0;
}
/*
* check for prime 53 <= factors < 257, if h*2^n-1 is large
* 2^276 > pfact(256)/pfact(52) > 2^275
*/
bits = highbit(testval);
if (bits >= 275) {
if (pprod256 <= 0) {
pprod256 = pfact(256)/pfact(52);
}
if (gcd(testval, pprod256) > 1) {
/* a small 53 <= prime < 257 divides h*2^n-1 */
ldebug("lucas",\
"not-prime: 53<=prime<257 divides h*2^n-1");
return 0;
}
}
/*
* try to compute u(2) (NOTE: some call this u(0))
*
* We will use gen_v1() to give us a v(1) using the values
* of 'h' and 'n'. We will then use gen_u2() to convert
* the v(1) into u(2).
*
* If gen_v1() returns a negative value, then we failed to
* generate a test for h*2^n-1. The legacy function,
* legacy_gen_v1() used by the Amdahl 6 could have returned
* -1. The new gen_v1() based on the method outlined in Ref4
* will never return -1 if h*2^n-1 is not a multiple of 3.
* Because the "multiple of 3" case is handled above, the
* call below to gen_v1() will never return -1.
*/
v1 = gen_v1(h, n);
if (v1 < 0) {
/* failure to test number */
print "unable to compute v(1) for", h : "*2^" : n : "-1";
ldebug("lucas", "unknown: no v(1)");
return -1;
}
u = gen_u2(h, n, v1);
/*
* compute u(n) (NOTE: some call this u(n-2))
*/
for (i=3; i <= n; ++i) {
/* u = (u^2 - 2) % testval; */
u = hnrmod(u^2 - 2, h, n, -1);
}
/*
* return 1 if prime, 0 is not prime
*/
if (u == 0) {
ldebug("lucas", "prime: end of test");
return 1;
} else {
ldebug("lucas", "not-prime: end of test");
return 0;
}
}
/*
* gen_u2 - determine the initial Lucas sequence for h*2^n-1
*
* Historically many start the Lucas sequence with u(0).
* Some, like the author of this code, prefer to start
* with U(2). This is so one may say:
*
* 2^p-1 is prime if u(p) = 0 mod 2^p-1
* or:
* h*2^p-1 is prime if u(p) = 0 mod h*2^p-1
*
* According to Ref1, Theorem 5:
*
* u(2) = alpha^h + alpha^(-h) (NOTE: Ref1 calls it u(0))
*
* Now:
*
* v(x) = alpha^x + alpha^(-x) (Ref1, bottom of page 872)
*
* Therefore:
*
* u(2) = v(h) (NOTE: Ref1 calls it u(0))
*
* We calculate v(h) as follows: (Ref1, top of page 873)
*
* v(0) = alpha^0 + alpha^(-0) = 2
* v(1) = alpha^1 + alpha^(-1) = gen_v1(h,n)
* v(n+2) = v(1)*v(n+1) - v(n)
*
* This function does not concern itself with the value of 'alpha'.
* The gen_v1() function is used to compute v(1), and identity
* functions take it from there.
*
* It can be shown that the following are true:
*
* v(2*n) = v(n)^2 - 2
* v(2*n+1) = v(n+1)*v(n) - v(1)
*
* To prevent v(x) from growing too large, one may replace v(x) with
* `v(x) mod h*2^n-1' at any time.
*
* See the function gen_v1() for details on the value of v(1).
*
* input:
* h - h as in h*2^n-1 (must be >= 1)
* n - n as in h*2^n-1 (must be >= 1)
* v1 - gen_v1(h,n) (must be >= 3) (see function below)
*
* returns:
* u(2) - initial value for Lucas test on h*2^n-1
* -1 - failed to generate u(2)
*/
define
gen_u2(h, n, v1)
{
local shiftdown; /* the power of 2 that divides h */
local r; /* low value: v(n) */
local s; /* high value: v(n+1) */
local hbits; /* highest bit set in h */
local oldh; /* pre-reduced h */
local oldn; /* pre-reduced n */
local i;
/*
* check arg types
*/
if (!isint(h)) {
quit "bad args: h must be an integer";
}
if (h < 0) {
quit "bad args: h must be an integer >= 1";
}
if (!isint(n)) {
quit "bad args: n must be an integer";
}
if (n < 1) {
quit "bad args: n must be an integer >= 1";
}
if (!isint(v1)) {
quit "bad args: v1 must be an integer";
}
if (v1 < 3) {
quit "bogus arg: v1 must be an integer >= 3";
}
/*
* reduce h if even
*
* we will force h to be odd by moving powers of two over to 2^n
*/
oldh = h;
oldn = n;
shiftdown = fcnt(h,2); /* h % 2^shiftdown == 0, max shiftdown */
if (shiftdown > 0) {
h >>= shiftdown;
n += shiftdown;
}
/*
* enforce the h > 0 and n >= 2 rules
*/
if (h <= 0 || n < 1) {
print " ERROR: h=":oldh, "n=":oldn, "reduced h=":h, "n=":n;
quit "reduced args violate the rule: 0 < h < 2^n";
}
hbits = highbit(h);
if (hbits >= n) {
print " ERROR: h=":oldh, "n=":oldn, "reduced h=":h, "n=":n;
quit "reduced args violate the rule: 0 < h < 2^n";
}
/*
* build up u2 based on the reversed bits of h
*/
/* setup for bit loop */
r = v1;
s = (r^2 - 2);
/*
* deal with small h as a special case
*
* The h value is odd > 0, and it needs to be
* at least 2 bits long for the loop below to work.
*/
if (h == 1) {
ldebug("gen_u2", "quick h == 1 case");
/* return r%(h*2^n-1); */
return hnrmod(r, h, n, -1);
}
/* cycle from second highest bit to second lowest bit of h */
for (i=hbits-1; i > 0; --i) {
/* bit(i) is 1 */
if (bit(h,i)) {
/* compute v(2n+1) = v(r+1)*v(r)-v1 */
/* r = (r*s - v1) % (h*2^n-1); */
r = hnrmod((r*s - v1), h, n, -1);
/* compute v(2n+2) = v(r+1)^2-2 */
/* s = (s^2 - 2) % (h*2^n-1); */
s = hnrmod((s^2 - 2), h, n, -1);
/* bit(i) is 0 */
} else {
/* compute v(2n+1) = v(r+1)*v(r)-v1 */
/* s = (r*s - v1) % (h*2^n-1); */
s = hnrmod((r*s - v1), h, n, -1);
/* compute v(2n) = v(r)^-2 */
/* r = (r^2 - 2) % (h*2^n-1); */
r = hnrmod((r^2 - 2), h, n, -1);
}
}
/* we know that h is odd, so the final bit(0) is 1 */
/* r = (r*s - v1) % (h*2^n-1); */
r = hnrmod((r*s - v1), h, n, -1);
/* compute the final u2 return value */
return r;
}
/*
* gen_u0 - determine the initial Lucas sequence for h*2^n-1
*
* Historically many start the Lucas sequence with u(0).
* Some, like the author of this code, prefer to start
* with u(2). This is so one may say:
*
* 2^p-1 is prime if u(p) = 0 mod 2^p-1
* or:
* h*2^n-1 is prime if U(n) = 0 mod h*2^n-1
*
* For those using the old code with gen_u0(), we
* simply call gen_u2() instead.
*
* See the function gen_u2() for details.
*
* input:
* h - h as in h*2^n-1 (must be >= 1)
* n - n as in h*2^n-1 (must be >= 1)
* v1 - gen_v1(h,n) (see function below)
*
* returns:
* u(2) - initial value for Lucas test on h*2^n-1
* -1 - failed to generate u(2)
*/
define
gen_u0(h, n, v1)
{
return gen_u2(h, n, v1);
}
/*
* rodseth_xhn - determine if v(1) == x for h*2^n-1
*
* For a given h*2^n-1, v(1) == x if:
*
* jacobi(x-2, h*2^n-1) == 1 (Ref4, condition 1) part 1
* jacobi(x+2, h*2^n-1) == -1 (Ref4, condition 1) part 2
*
* Now when x-2 <= 0:
*
* jacobi(x-2, h*2^n-1) == 0
*
* because:
*
* jacobi(x,y) == 0 if x <= 0
*
* So for (Ref4, condition 1) part 1 to be true:
*
* x-2 > 0
*
* And therefore:
*
* x > 2
*
* input:
* x potential v(1) value
* h h as in h*2^n-1 (h must be odd >= 1)
* n n as in h*2^n-1 (must be >= 1)
*
* returns:
* 1 if v(1) == x for h*2^n-1
* 0 otherwise
*/
define
rodseth_xhn(x, h, n)
{
local testval; /* h*2^n-1 */
/*
* check arg types
*/
if (!isint(h)) {
quit "bad args: h must be an integer";
}
if (iseven(h)) {
quit "bad args: h must be an odd integer";
}
if (h < 1) {
quit "bad args: h must be an integer >= 1";
}
if (!isint(n)) {
quit "bad args: n must be an integer";
}
if (n < 1) {
quit "bad args: n must be an integer >= 1";
}
if (!isint(x)) {
quit "bad args: x must be an integer";
}
/*
* firewall
*/
if (x <= 2) {
return 0;
}
/*
* Check for jacobi(x-2, h*2^n-1) == 1 (Ref4, condition 1) part 1
*/
testval = h*2^n-1;
if (jacobi(x-2, testval) != 1) {
return 0;
}
/*
* Check for jacobi(x+2, h*2^n-1) == -1 (Ref4, condition 1) part 2
*/
if (jacobi(x+2, testval) != -1) {
return 0;
}
/*
* v(1) == x for this h*2^n-1
*/
return 1;
}
/*
* Trial tables used by gen_v1()
*
* When h mod 3 == 0, according to Ref4 we need to find the first value X where:
*
* jacobi(X-2, h*2^n-1) == 1 (Ref4, condition 1) part 1
* jacobi(X+2, h*2^n-1) == -1 (Ref4, condition 1) part 2
*
* We can show that X > 2. See the comments in the rodseth_xhn(x,h,n) above.
*
* Some values of X satisfy more often than others. For example a large sample
* of h*2^n-1, h odd multiple of 3, and large n (some around 1e4, some near 1e6,
* others near 3e7) where the sample size was 66 973 365, here is the count of
* the smallest value of X that satisfies conditions in Ref4, condition 1:
*
* count X
* ----------
* 26791345 3
* 17223016 5
* 7829600 9
* 6988774 11
* 3301093 15
* 1517149 17
* 910346 21
* 711791 29
* 573403 20
* 390395 27
* 288637 35
* 149751 36
* 107733 39
* 58743 41
* 35619 45
* 25052 32
* 17775 51
* 13031 44
* 7563 56
* 7540 49
* 7060 59
* 4407 57
* 2948 65
* 2502 55
* 2388 69
* 2094 71
* 689 77
* 626 81
* 491 66
* 426 95
* 219 80
* 203 67
* 185 84
* 152 99
* 127 72
* 102 74
* 98 87
* 67 90
* 55 104
* 48 101
* 32 105
* 17 109
* 16 116
* 15 111
* 13 92
* 12 125
* 7 129
* 3 146
* 2 140
* 2 120
* 1 165
* 1 161
* 1 155
*
* It is important that we select the smallest possible v(1). While testing
* various values of X for V(1) is fast, using larger than necessary values
* of V(1) of can slow down calculating V(h).
*
* The above distribution was found to hold fairly well over many values of
* odd h that are also a multiple of 3 and for many values of n where h < 2^n.
*
* For example for in a sample size of 1000000 numbers of the form h*2^n-1
* where h is an odd multiple of 3, 12996351 <= h <= 13002351,
* 4331116 <= n <= 4332116, these are the smallest v(1) values that were found:
*
* smallest percentage
* v(1) used
* -------- ---------
* 3 40.0000 %
* 5 25.6833 %
* 9 11.6924 %
* 11 10.4528 %
* 15 4.8048 %
* 17 2.3458 %
* 21 1.3734 %
* 29 1.0527 %
* 20 0.8595 %
* 27 0.5758 %
* 35 0.4420 %
* 36 0.2433 %
* 39 0.1779 %
* 41 0.0885 %
* 45 0.0571 %
* 32 0.0337 %
* 51 0.0289 %
* 44 0.0205 %
* 49 0.0176 %
* 56 0.0137 %
* 59 0.0108 %
* 57 0.0053 %
* 65 0.0047 %
* 55 0.0045 %
* 69 0.0031 %
* 71 0.0024 %
* 66 0.0011 %
* 95 0.0008 %
* 81 0.0008 %
* 77 0.0006 %
* 72 0.0005 %
* 99 0.0004 %
* 80 0.0003 %
* 74 0.0003 %
* 84 0.0002 %
* 67 0.0002 %
* 87 0.0001 %
* 104 0.0001 %
* 129 0.0001 %
*
* However, a case can be made for considering only odd values for v(1)
* candidates. When h * 2^n-1 is prime and h is an odd multiple of 3,
* a smallest v(1) that is even is extremely rate. Of the list of 146553
* known primes of the form h*2^n-1 when h is an odd a multiple of 3,
* none has an smallest v(1) that was even.
*
* See:
*
* https://github.com/arcetri/verified-prime
*
* for that list of 146553 known primes of the form h*2^n-1.
*
* That same example for in a sample size of 1000000 numbers of the
* form h*2^n-1 where h is an odd multiple of 3, 12996351 <= h <= 13002351,
* 4331116 <= n <= 4332116, these are the smallest odd v(1) values that were
* found:
*
* smallest percentage
* odd v(1) used
* -------- ---------
* 3 40.0000 %
* 5 25.6833 %
* 9 11.6924 %
* 11 10.4528 %
* 15 4.8048 %
* 17 2.3458 %
* 21 1.6568 %
* 29 1.6174 %
* 35 0.4529 %
* 27 0.3546 %
* 39 0.3470 %
* 41 0.2159 %
* 45 0.1173 %
* 31 0.0661 %
* 51 0.0619 %
* 55 0.0419 %
* 59 0.0250 %
* 49 0.0170 %
* 69 0.0110 %
* 65 0.0098 %
* 71 0.0078 %
* 85 0.0048 %
* 81 0.0044 %
* 95 0.0038 %
* 99 0.0021 %
* 125 0.0009 %
* 57 0.0007 %
* 111 0.0005 %
* 77 0.0003 %
* 165 0.0003 %
* 155 0.0002 %
* 129 0.0002 %
* 101 0.0002 %
* 53 0.0001 %
*
* Moreover when evaluating odd candidates for v(1), one may cache Jacobi
* symbol evaluations to reduce the number of Jacobi symbol evaluations to
* a minimum. For example, if one tests 5 and finds that the 2nd case fails:
*
* jacobi(5+2, h*2^n-1) != -1
*
* Then if one is later testing 9, the Jacobi symbol value for the first
* 1st case:
*
* jacobi(7-2, h*2^n-1)
*
* is already known.
*
* Without Jacobi symbol value caching, it requires on average
* 4.851377 Jacobi symbol evaluations. With Jacobi symbol value caching
* cacheing, an averare of 4.348820 Jacobi symbol evaluations is needed.
*
* Given this information, when odd h is a multiple of 3 we try, in order,
* these odd values of X:
*
* 3, 5, 9, 11, 15, 17, 21, 29, 27, 35, 39, 41, 31, 45, 51, 55, 49, 59,
* 69, 65, 71, 57, 85, 81, 95, 99, 77, 53, 67, 125, 111, 105, 87, 129,
* 101, 83, 165, 155, 149, 141, 121, 109
*
* And stop on the first value of X where:
*
* jacobi(X-2, h*2^n-1) == 1
* jacobi(X+2, h*2^n-1) == -1
*
* Less than 1 case out of 1000000 will not be satisifed by the above list.
* If no value in that list works, we start simple search starting with X = 167
* and incrementing by 2 until a value of X is found.
*
* The x_tbl[] matrix contains those values of X to try in order.
* If all x_tbl_len fail to satisfy Ref4 condition 1 (this happens less than
* 1 in 1000000 cases), then we begin a linear search of odd values starting at
* next_x until we find a proper X value.
*/
x_tbl_len = 42;
mat x_tbl[x_tbl_len];
x_tbl = {
3, 5, 9, 11, 15, 17, 21, 29, 27, 35, 39, 41, 31, 45, 51, 55, 49, 59,
69, 65, 71, 57, 85, 81, 95, 99, 77, 53, 67, 125, 111, 105, 87, 129,
101, 83, 165, 155, 149, 141, 121, 109
};
next_x = 167; /* must be 2 more than the largest value in x_tbl[] */
/*
* gen_v1 - compute the v(1) for a given h*2^n-1 if we can
*
* This function assumes:
*
* n > 2 (n==2 has already been eliminated)
* h mod 2 == 1
* h < 2^n
* h*2^n-1 mod 3 != 0 (h*2^n-1 has no small factors, such as 3)
*
* The generation of v(1) depends on the value of h. There are two cases
* to consider, h mod 3 != 0, and h mod 3 == 0.
*
***
*
* Case 1: (h mod 3 != 0)
*
* This case is easy.
*
* In Ref1, page 869, one finds that if: (or see Ref2, page 131-132)
*
* h mod 6 == +/-1
* h*2^n-1 mod 3 != 0
*
* which translates, gives the functions assumptions, into the condition:
*
* h mod 3 != 0
*
* If this case condition is true, then:
*
* u(2) = (2+sqrt(3))^h + (2-sqrt(3))^h (see Ref1, page 869)
* = (2+sqrt(3))^h + (2+sqrt(3))^(-h) (NOTE: some call this u(2))
*
* and since Ref1, Theorem 5 states:
*
* u(2) = alpha^h + alpha^(-h) (NOTE: some call this u(2))
* r = abs(2^2 - 1^2*3) = 1
*
* where these values work for Case 1: (h mod 3 != 0)
*
* a = 1
* b = 2
* D = 1
*
* Now at the bottom of Ref1, page 872 states:
*
* v(x) = alpha^x + alpha^(-x)
*
* If we let:
*
* alpha = (2+sqrt(3))
*
* then
*
* u(2) = v(h) (NOTE: some call this u(2))
*
* so we can always return
*
* v(1) = alpha^1 + alpha^(-1)
* = (2+sqrt(3)) + (2-sqrt(3))
* = 4
*
* In 40% of the cases when h is not a multiple of 3, 3 is a valid value
* for v(1). We can test if 3 is a valid value for v(1) in this case:
*
* if jacobi(1, h*2^n-1) == 1 and jacobi(5, h*2^n-1) == -1, then
* v(1) = 3
* else
* v(1) = 4
*
* HOTE: The above "if then else" works only of h is not a multiple of 3.
*
***
*
* Case 2: (h mod 3 == 0)
*
* For the case where h is a multiple of 3, we turn to Ref4.
*
* The central theorem on page 3 of that paper states that
* we may set v(1) to the first value X that satisfies:
*
* jacobi(X-2, h*2^n-1) == 1 (Ref4, condition 1)
* jacobi(X+2, h*2^n-1) == -1 (Ref4, condition 1)
*
* NOTE: Ref4 uses P, which we shall refer to as X.
* Ref4 uses N, which we shall refer to as h*2^n-1.
*
* NOTE: Ref4 uses the term Legendre-Jacobi symbol, which
* we shall refer to as the Jacobi symbol.
*
* Before we address the two conditions, we need some background information
* on two symbols, Legendre and Jacobi. In Ref 2, pp 278, 284-285, we find
* the following definitions of jacobi(a,b) and L(a,p):
*
* The Legendre symbol L(a,p) takes the value:
*
* L(a,p) == 1 => a is a quadratic residue of p
* L(a,p) == -1 => a is NOT a quadratic residue of p
*
* when:
*
* p is prime
* p mod 2 == 1
* gcd(a,p) == 1
*
* The value a is a quadratic residue of b if there exists some integer z
* such that:
*
* z^2 mod b == a
*
* The Jacobi symbol jacobi(a,b) takes the value:
*
* jacobi(a,b) == 1 => b is not prime,
* or a is a quadratic residue of b
* jacobi(a,b) == -1 => a is NOT a quadratic residue of b
*
* when
*
* b mod 2 == 1
* gcd(a,b) == 1
*
* It is worth noting for the Legendre symbol, in order for L(X+/-2,
* h*2^n-1) to be defined, we must ensure that neither X-2 nor X+2 are
* factors of h*2^n-1. This is done by pre-screening h*2^n-1 to not
* have small factors and keeping X+2 less than that small factor
* limit. It is worth noting that in lucas(h, n), we first verify
* that h*2^n-1 does not have a factor < 257 before performing the
* primality test. So while X+/-2 < 257, we know that
* gcd(X+/-2, h*2^n-1) == 1.
*
* Returning to the testing of conditions in Ref4, condition 1:
*
* jacobi(X-2, h*2^n-1) == 1
* jacobi(X+2, h*2^n-1) == -1
*
* When such an X is found, we set:
*
* v(1) = X
*
***
*
* In conclusion, we can compute v,(1) by attempting to do the following:
*
* h mod 3 != 0
*
* we return:
*
* v(1) == 4
*
* h mod 3 == 0
*
* we return:
*
* v(1) = X
*
* where X > 2 in a integer such that:
*
* jacobi(X-2, h*2^n-1) == 1
* jacobi(X+2, h*2^n-1) == -1
*
***
*
* input:
* h h as in h*2^n-1 (h must be odd >= 1)
* n n as in h*2^n-1 (must be >= 1)
*
* output:
* returns v(1), or
* -1 when h*2^n-1 is a multiple of 3
*/
define
gen_v1(h, n)
{
local x; /* potential v(1) to test */
local i; /* x_tbl index */
local v1m2; /* X-2 1st case */
local v1p2; /* X+2 2nd case */
local testval; /* h*2^n-1 - value we are testing if prime */
local mat cached_v1[next_x]; /* cached Jacobi symbol values or 0 */
/*
* check arg types
*/
if (!isint(h)) {
quit "bad args: h must be an integer";
}
if (iseven(h)) {
quit "bad args: h must be an odd integer";
}
if (h < 1) {
quit "bad args: h must be an integer >= 1";
}
if (!isint(n)) {
quit "bad args: n must be an integer";
}
if (n < 1) {
quit "bad args: n must be an integer >= 1";
}
/*
* pretest: Verify that h*2^n-1 is not a multiple of 3
*/
if (((h % 3 == 1) && (n % 2 == 0)) || ((h % 3 == 2) && (n % 2 == 1))) {
/* no need to test h*2^n-1, it is not prime */
return -1;
}
/*
* Common Mersenne number case:
*
* For Mersenne numbers:
*
* 2^n-1
*
* we can use, 40% of the time, v(1) == 3. However nearly all code that
* implements the Lucas-Lehmer test uses v(1) == 4. Whenever for
* h != 0 mod 3, and particular the Mersenne number case of when h == 1:
*
* 1*2^n-1
*
* v(1) == 4 always works. For this reason, we return 4 when h == 1.
*/
if (h == 1) {
/* v(1) == 4 always works for the Mersenne number case */
return 4;
}
/*
* check for Case 1: (h mod 3 != 0)
*/
if (h % 3 != 0) {
if (rodseth_xhn(3, h, n) == 1) {
/* 40% of the time, 3 works when h mod 3 != 0 */
return 3;
} else {
/* otherwise 4 always works when h mod 3 != 0 */
return 4;
}
}
/*
* What follow is Case 2: (h mod 3 == 0)
*/
/*
* clear cache
*/
matfill(cached_v1, 0);
/*
* We will look for x that satisfies conditions in Ref4, condition 1:
*
* jacobi(X-2, h*2^n-1) == 1 part 1
* jacobi(X+2, h*2^n-1) == -1 part 2
*
* NOTE: If we wanted to be super optimial, we would cache
* jacobi(X+2, h*2^n-1) that that when we increment X
* to the next odd value, the now jacobi(X-2, h*2^n-1)
* does not need to be re-evaluted.
*/
testval = h*2^n-1;
for (i=0; i < x_tbl_len; ++i) {
/*
* obtain the next test candidate
*/
x = x_tbl[i];
/*
* Check x for condition 1 part 1
*
* jacobi(x-2, h*2^n-1) == 1
*/
v1m2 = x-2;
if (cached_v1[v1m2] == 0) {
cached_v1[v1m2] = jacobi(v1m2, testval);
}
if (cached_v1[v1m2] != 1) {
continue;
}
/*
* Check x for condition 1 part 2
*
* jacobi(x+2, h*2^n-1) == -1
*/
v1p2 = x+2;
if (cached_v1[v1p2] == 0) {
cached_v1[v1p2] = jacobi(v1p2, testval);
}
if (cached_v1[v1p2] != -1) {
continue;
}
/*
* found a x that satisfies Ref4 condition 1
*/
ldebug("gen_v1", "h= " + str(h) + " n= " + str(n) +
" v1= " + str(x) + " using tbl[ " +
str(i) + " ]");
return x;
}
/*
* We are in that rare case (less than 1 in 1 000 000) where none of the
* common X values satisfy Ref4 condition 1. We start a linear search
* of odd vules at next_x from here on.
*/
x = next_x;
while (rodseth_xhn(x, h, n) != 1) {
x += 2;
}
/* finally found a v(1) value */
ldebug("gen_v1", "h= " + str(h) + " n= " + str(n) +
" v1= " + str(x) + " beyond tbl");
return x;
}
/*
* ldebug - print a debug statement
*
* input:
* funct name of calling function
* str string to print
*/
define
ldebug(funct, str)
{
if (config("resource_debug") & 8) {
print "DEBUG:", funct:":", str;
}
return;
}
/*
************************
* Amdahl 6 legacy code *
************************
*
* NOTE: What follows is legacy code based on the method used by the
* Amdahl 6 group:
*
* John Brown, Landon Curt Noll, Bodo Parady, Gene Smith,
* Joel Smith and Sergio Zarantonello
*
* This method generated v(1) for nearly all values, except for a
* few rare cases when h mod 3 == 0. The code is NOT used by lucas.cal
* above. The gen_v1() function above is based on an improved method
* outlined in Ref4. That method generated v(1) for all h.
*
* The code below is kept for historical purposes only. The functions
* and global variables of the Amdahl 6 legacy code all begin with legacy_.
*/
/*
* Trial tables used by legacy_gen_v1()
*
* When h mod 3 == 0, one needs particular values of D, a and b (see
* legacy_gen_v1 documentation) in order to find a value of v(1).
*
* This table defines 'legacy_quickmax' possible tests to be taken in ascending
* order. The legacy_v1_qval[x] refers to a v(1) value from Ref1, Table 1. A
* related D value is found in legacy_d_qval[x]. All D values expect
* legacy_d_qval[1] are also taken from Ref1, Table 1. The case of D == 21 as
* listed in Ref1, Table 1 can be changed to D == 7 for the sake of the test
* because of {note 6}.
*
* It should be noted that the D values all satisfy the selection values
* as outlined in the legacy_gen_v1() function comments. That is:
*
* D == P*(2^f)*(3^g)
*
* where f == 0 and g == 0, P == D. So we simply need to check that
* one of the following two cases are true:
*
* P mod 4 == 1 and J(h*2^n-1 mod P, P) == -1
* P mod 4 == -1 and J(h*2^n-1 mod P, P) == 1
*
* In all cases, the value of r is:
*
* r == Q*(2^j)*(3^k)*(z^2)
*
* where Q == 1. No further processing is needed to compute v(1) when r
* is of this form.
*/
legacy_quickmax = 8;
mat legacy_d_qval[legacy_quickmax];
mat legacy_v1_qval[legacy_quickmax];
legacy_d_qval[0] = 5; legacy_v1_qval[0] = 3; /* a=1 b=1 r=4 */
legacy_d_qval[1] = 7; legacy_v1_qval[1] = 5; /* a=3 b=1 r=12 D=21 */
legacy_d_qval[2] = 13; legacy_v1_qval[2] = 11; /* a=3 b=1 r=4 */
legacy_d_qval[3] = 11; legacy_v1_qval[3] = 20; /* a=3 b=1 r=2 */
legacy_d_qval[4] = 29; legacy_v1_qval[4] = 27; /* a=5 b=1 r=4 */
legacy_d_qval[5] = 53; legacy_v1_qval[5] = 51; /* a=53 b=1 r=4 */
legacy_d_qval[6] = 17; legacy_v1_qval[6] = 66; /* a=17 b=1 r=1 */
legacy_d_qval[7] = 19; legacy_v1_qval[7] = 74; /* a=38 b=1 r=2 */
/*
* legacy_gen_v1 - compute the v(1) for a given h*2^n-1 if we can
*
* This function assumes:
*
* n > 2 (n==2 has already been eliminated)
* h mod 2 == 1
* h < 2^n
* h*2^n-1 mod 3 != 0 (h*2^n-1 has no small factors, such as 3)
*
* The generation of v(1) depends on the value of h. There are two cases
* to consider, h mod 3 != 0, and h mod 3 == 0.
*
***
*
* Case 1: (h mod 3 != 0)
*
* This case is easy and always finds v(1).
*
* In Ref1, page 869, one finds that if: (or see Ref2, page 131-132)
*
* h mod 6 == +/-1
* h*2^n-1 mod 3 != 0
*
* which translates, gives the functions assumptions, into the condition:
*
* h mod 3 != 0
*
* If this case condition is true, then:
*
* u(2) = (2+sqrt(3))^h + (2-sqrt(3))^h (see Ref1, page 869)
* = (2+sqrt(3))^h + (2+sqrt(3))^(-h) (some call this u(0))
*
* and since Ref1, Theorem 5 states:
*
* u(2) = alpha^h + alpha^(-h)
* r = abs(2^2 - 1^2*3) = 1
*
* where these values work for Case 1: (h mod 3 != 0)
*
* a = 1
* b = 2
* D = 1
*
* Now at the bottom of Ref1, page 872 states:
*
* v(x) = alpha^x + alpha^(-x)
*
* If we let:
*
* alpha = (2+sqrt(3))
*
* then
*
* u(2) = v(h)
*
* so we simply return
*
* v(1) = alpha^1 + alpha^(-1)
* = (2+sqrt(3)) + (2-sqrt(3))
* = 4
*
***
*
* Case 2: (h mod 3 == 0)
*
* This case is not so easy and finds v(1) in most all cases. In this
* version of this program, we will simply return -1 (failure) if we
* hit one of the cases that fall thru the cracks. This does not happen
* often, so this is not too bad.
*
* Ref1, Theorem 5 contains the following definitions:
*
* r = abs(a^2 - b^2*D)
* alpha = (a + b*sqrt(D))^2/r
*
* where D is 'square free', and 'alpha = epsilon^s' (for some s>0) are units
* in the quadratic field K(sqrt(D)).
*
* One can find possible values for a, b and D in Ref1, Table 1 (page 872).
* (see the file lucas_tbl.cal)
*
* Now Ref1, Theorem 5 states that if:
*
* L(D, h*2^n-1) = -1 [condition 1]
* L(r, h*2^n-1) * (a^2 - b^2*D)/r = -1 [condition 2]
*
* where L(x,y) is the Legendre symbol (see below), then:
*
* u(2) = alpha^h + alpha^(-h)
*
* The bottom of Ref1, page 872 states:
*
* v(x) = alpha^x + alpha^(-x)
*
* thus since:
*
* u(2) = v(h)
*
* so we want to return:
*
* v(1) = alpha^1 + alpha^(-1)
*
* Therefore we need to take a given (D,a,b), determine if the two conditions
* are true, and return the related v(1).
*
* Before we address the two conditions, we need some background information
* on two symbols, Legendre and Jacobi. In Ref 2, pp 278, 284-285, we find
* the following definitions of J(a,p) and L(a,n):
*
* The Legendre symbol L(a,p) takes the value:
*
* L(a,p) == 1 => a is a quadratic residue of p
* L(a,p) == -1 => a is NOT a quadratic residue of p
*
* when
*
* p is prime
* p mod 2 == 1
* gcd(a,p) == 1
*
* The value x is a quadratic residue of y if there exists some integer z
* such that:
*
* z^2 mod y == x
*
* The Jacobi symbol J(x,y) takes the value:
*
* J(x,y) == 1 => y is not prime, or x is a quadratic residue of y
* J(x,y) == -1 => x is NOT a quadratic residue of y
*
* when
*
* y mod 2 == 1
* gcd(x,y) == 1
*
* In the following comments on Legendre and Jacobi identities, we shall
* assume that the arguments to the symbolic are valid over the symbol
* definitions as stated above.
*
* In Ref2, pp 280-284, we find that:
*
* L(a,p)*L(b,p) == L(a*b,p) {A3.5}
* J(x,y)*J(z,y) == J(x*z,y) {A3.14}
* L(a,p) == L(p,a) * (-1)^((a-1)*(p-1)/4) {A3.8}
* J(x,y) == J(y,x) * (-1)^((x-1)*(y-1)/4) {A3.17}
*
* The equality L(a,p) == J(a,p) when: {note 0}
*
* p is prime
* p mod 2 == 1
* gcd(a,p) == 1
*
* It can be shown that (see Ref3):
*
* L(a,p) == L(a mod p, p) {note 1}
* L(z^2, p) == 1 {note 2}
*
* From Ref2, table 32:
*
* p mod 8 == +/-1 implies L(2,p) == 1 {note 3}
* p mod 12 == +/-1 implies L(3,p) == 1 {note 4}
*
* Since h*2^n-1 mod 8 == -1, for n>2, note 3 implies:
*
* L(2, h*2^n-1) == 1 (n>2) {note 5}
*
* Since h=3*A, h*2^n-1 mod 12 == -1, for A>0, note 4 implies:
*
* L(3, h*2^n-1) == 1 {note 6}
*
* By use of {A3.5}, {note 2}, {note 5} and {note 6}, one can show:
*
* L((2^g)*(3^l)*(z^2), h*2^n-1) == 1 (g>=0,l>=0,z>0,n>2) {note 7}
*
* Returning to the testing of conditions, take condition 1:
*
* L(D, h*2^n-1) == -1 [condition 1]
*
* In order for J(D, h*2^n-1) to be defined, we must ensure that D
* is not a factor of h*2^n-1. This is done by pre-screening h*2^n-1 to
* not have small factors and selecting D less than that factor check limit.
*
* By use of {note 7}, we can show that when we choose D to be:
*
* D is square free
* D = P*(2^f)*(3^g) (P is prime>2)
*
* The square free condition implies f = 0 or 1, g = 0 or 1. If f and g
* are both 1, P must be a prime > 3.
*
* So given such a D value:
*
* L(D, h*2^n-1) == L(P*(2^g)*(3^l), h*2^n-1)
* == L(P, h*2^n-1) * L((2^g)*(3^l), h*2^n-1) {A3.5}
* == L(P, h*2^n-1) * 1 {note 7}
* == L(h*2^n-1, P)*(-1)^((h*2^n-2)*(P-1)/4) {A3.8}
* == L(h*2^n-1 mod P, P)*(-1)^((h*2^n-2)*(P-1)/4) {note 1}
* == J(h*2^n-1 mod P, P)*(-1)^((h*2^n-2)*(P-1)/4) {note 0}
*
* When does J(h*2^n-1 mod P, P)*(-1)^((h*2^n-2)*(P-1)/4) take the value of -1,
* thus satisfy [condition 1]? The answer depends on P. Now P is a prime>2,
* thus P mod 4 == 1 or -1.
*
* Take P mod 4 == 1:
*
* P mod 4 == 1 implies (-1)^((h*2^n-2)*(P-1)/4) == 1
*
* Thus:
*
* L(D, h*2^n-1) == L(h*2^n-1 mod P, P) * (-1)^((h*2^n-2)*(P-1)/4)
* == L(h*2^n-1 mod P, P)
* == J(h*2^n-1 mod P, P)
*
* Take P mod 4 == -1:
*
* P mod 4 == -1 implies (-1)^((h*2^n-2)*(P-1)/4) == -1
*
* Thus:
*
* L(D, h*2^n-1) == L(h*2^n-1 mod P, P) * (-1)^((h*2^n-2)*(P-1)/4)
* == L(h*2^n-1 mod P, P) * -1
* == -J(h*2^n-1 mod P, P)
*
* Therefore [condition 1] is met if, and only if, one of the following
* to cases are true:
*
* P mod 4 == 1 and J(h*2^n-1 mod P, P) == -1
* P mod 4 == -1 and J(h*2^n-1 mod P, P) == 1
*
* Now consider [condition 2]:
*
* L(r, h*2^n-1) * (a^2 - b^2*D)/r == -1 [condition 2]
*
* We select only a, b, r and D values where:
*
* (a^2 - b^2*D)/r == -1
*
* Therefore in order for [condition 2] to be met, we must show that:
*
* L(r, h*2^n-1) == 1
*
* If we select r to be of the form:
*
* r == Q*(2^j)*(3^k)*(z^2) (Q == 1, j>=0, k>=0, z>0)
*
* then by use of {note 7}:
*
* L(r, h*2^n-1) == L(Q*(2^j)*(3^k)*(z^2), h*2^n-1)
* == L((2^j)*(3^k)*(z^2), h*2^n-1)
* == 1 {note 2}
*
* and thus, [condition 2] is met.
*
* If we select r to be of the form:
*
* r == Q*(2^j)*(3^k)*(z^2) (Q is prime>2, j>=0, k>=0, z>0)
*
* then by use of {note 7}:
*
* L(r, h*2^n-1) == L(Q*(2^j)*(3^k)*(z^2), h*2^n-1)
* == L(Q, h*2^n-1) * L((2^j)*(3^k)*(z^2), h*2^n-1) {A3.5}
* == L(Q, h*2^n-1) * 1 {note 2}
* == L(h*2^n-1, Q) * (-1)^((h*2^n-2)*(Q-1)/4) {A3.8}
* == L(h*2^n-1 mod Q, Q)*(-1)^((h*2^n-2)*(Q-1)/4) {note 1}
* == J(h*2^n-1 mod Q, Q)*(-1)^((h*2^n-2)*(Q-1)/4) {note 0}
*
* When does J(h*2^n-1 mod Q, Q)*(-1)^((h*2^n-2)*(Q-1)/4) take the value of 1,
* thus satisfy [condition 2]? The answer depends on Q. Now Q is a prime>2,
* thus Q mod 4 == 1 or -1.
*
* Take Q mod 4 == 1:
*
* Q mod 4 == 1 implies (-1)^((h*2^n-2)*(Q-1)/4) == 1
*
* Thus:
*
* L(D, h*2^n-1) == L(h*2^n-1 mod Q, Q) * (-1)^((h*2^n-2)*(Q-1)/4)
* == L(h*2^n-1 mod Q, Q)
* == J(h*2^n-1 mod Q, Q)
*
* Take Q mod 4 == -1:
*
* Q mod 4 == -1 implies (-1)^((h*2^n-2)*(Q-1)/4) == -1
*
* Thus:
*
* L(D, h*2^n-1) == L(h*2^n-1 mod Q, Q) * (-1)^((h*2^n-2)*(Q-1)/4)
* == L(h*2^n-1 mod Q, Q) * -1
* == -J(h*2^n-1 mod Q, Q)
*
* Therefore [condition 2] is met by selecting D = Q*(2^j)*(3^k)*(z^2),
* where Q is prime>2, j>=0, k>=0, z>0; if and only if one of the following
* to cases are true:
*
* Q mod 4 == 1 and J(h*2^n-1 mod Q, Q) == 1
* Q mod 4 == -1 and J(h*2^n-1 mod Q, Q) == -1
*
***
*
* In conclusion, we can compute v(1) by attempting to do the following:
*
* h mod 3 != 0
*
* we return:
*
* v(1) == 4
*
* h mod 3 == 0
*
* define:
*
* r == abs(a^2 - b^2*D)
* alpha == (a + b*sqrt(D))^2/r
*
* we return:
*
* v(1) = alpha^1 + alpha^(-1)
*
* if and only if we can find a given a, b, D that obey all the
* following selection rules:
*
* D is square free
*
* D == P*(2^f)*(3^g) (P is prime>2, f,g == 0 or 1)
*
* (a^2 - b^2*D)/r == -1
*
* r == Q*(2^j)*(3^k)*(z^2) (Q==1 or Q is prime>2, j>=0, k>=0, z>0)
*
* one of the following is true:
* P mod 4 == 1 and J(h*2^n-1 mod P, P) == -1
* P mod 4 == -1 and J(h*2^n-1 mod P, P) == 1
*
* if Q is prime, then one of the following is true:
* Q mod 4 == 1 and J(h*2^n-1 mod Q, Q) == 1
* Q mod 4 == -1 and J(h*2^n-1 mod Q, Q) == -1
*
* If we cannot find a v(1) quickly enough, then we will give up
* testing h*2^n-1. This does not happen too often, so this hack
* is not too bad.
*
***
*
* input:
* h h as in h*2^n-1 (must be >= 1)
* n n as in h*2^n-1 (must be >= 1)
*
* output:
* returns v(1), or -1 is there is no quick way
*/
define
legacy_gen_v1(h, n)
{
local d; /* the 'D' value to try */
local val_mod; /* h*2^n-1 mod 'D' */
local i;
/*
* check arg types
*/
if (!isint(h)) {
quit "bad args: h must be an integer";
}
if (h < 1) {
quit "bad args: h must be an integer >= 1";
}
if (!isint(n)) {
quit "bad args: n must be an integer";
}
if (n < 1) {
quit "bad args: n must be an integer >= 1";
}
/*
* check for case 1
*/
if (h % 3 != 0) {
/* v(1) is easy to compute */
return 4;
}
/*
* We will try all 'D' values until we find a proper v(1)
* or run out of 'D' values.
*/
for (i=0; i < legacy_quickmax; ++i) {
/* grab our 'D' value */
d = legacy_d_qval[i];
/* compute h*2^n-1 mod 'D' quickly */
val_mod = (h*pmod(2,n%(d-1),d)-1) % d;
/*
* if 'D' mod 4 == 1, then
* (h*2^n-1) mod 'D' can not be a quadratic residue of 'D'
* else
* (h*2^n-1) mod 'D' must be a quadratic residue of 'D'
*/
if (d%4 == 1) {
/* D mod 4 == 1, so check for J(D, h*2^n-1) == -1 */
if (jacobi(val_mod, d) == -1) {
/* it worked, return the related v(1) value */
return legacy_v1_qval[i];
}
} else {
/* D mod 4 == -1, so check for J(D, h*2^n-1) == 1 */
if (jacobi(val_mod, d) == 1) {
/* it worked, return the related v(1) value */
return legacy_v1_qval[i];
}
}
}
/*
* This is an example of a more complex proof construction.
* The code above will not be able to find the v(1) for:
*
* 81*2^81-1
*
* We will check with:
*
* v(1)=81 D=6557 a=79 b=1 r=316
*
* Now, D==79*83 and r=79*2^2. If we show that:
*
* J(h*2^n-1 mod 79, 79) == -1
* J(h*2^n-1 mod 83, 83) == 1
*
* then we will satisfy [condition 1]. Observe:
*
* 79 mod 4 == -1 implies (-1)^((h*2^n-2)*(79-1)/4) == -1
* 83 mod 4 == -1 implies (-1)^((h*2^n-2)*(83-1)/4) == -1
*
* J(D, h*2^n-1) == J(83, h*2^n-1) * J(79, h*2^n-1)
* == J(h*2^n-1, 83) * (-1)^((h*2^n-2)*(83-1)/4) *
* J(h*2^n-1, 79) * (-1)^((h*2^n-2)*(79-1)/4)
* == J(h*2^n-1 mod 83, 83) * -1 *
* J(h*2^n-1 mod 79, 79) * -1
* == 1 * -1 *
* -1 * -1
* == -1
*
* We will also satisfy [condition 2]. Observe:
*
* (a^2 - b^2*D)/r == (79^2 - 1^1*6557)/316
* == -1
*
* L(r, h*2^n-1) == L(Q*(2^j)*(3^k)*(z^2), h*2^n-1)
* == L(79, h*2^n-1) * L(2^2, h*2^n-1)
* == L(79, h*2^n-1) * 1
* == L(h*2^n-1, 79) * (-1)^((h*2^n-2)*(79-1)/4)
* == L(h*2^n-1, 79) * -1
* == L(h*2^n-1 mod 79, 79) * -1
* == J(h*2^n-1 mod 79, 79) * -1
* == -1 * -1
* == 1
*/
if (jacobi( ((h*pmod(2,n%(79-1),79)-1)%79), 79 ) == -1 &&
jacobi( ((h*pmod(2,n%(83-1),83)-1)%83), 83 ) == 1) {
/* return the associated v(1)=81 */
return 81;
}
/* no quick and dirty v(1), so return -1 */
return -1;
}