/* * lucas - perform a Lucas primality test on h*2^n-1 * * Copyright (C) 1999,2017,2018,2021 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 /\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 0 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 average 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 satisfied 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 * * NOTE: 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 optimal, 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-evaluated. */ 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 values 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; }