1/* 2 * lucas - perform a Lucas primality test on h*2^n-1 3 * 4 * Copyright (C) 1999,2017,2018,2021 Landon Curt Noll 5 * 6 * Calc is open software; you can redistribute it and/or modify it under 7 * the terms of the version 2.1 of the GNU Lesser General Public License 8 * as published by the Free Software Foundation. 9 * 10 * Calc is distributed in the hope that it will be useful, but WITHOUT 11 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 12 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General 13 * Public License for more details. 14 * 15 * A copy of version 2.1 of the GNU Lesser General Public License is 16 * distributed with calc under the filename COPYING-LGPL. You should have 17 * received a copy with calc; if not, write to Free Software Foundation, Inc. 18 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. 19 * 20 * Under source code control: 1990/05/03 16:49:51 21 * File existed as early as: 1990 22 * 23 * chongo <was here> /\oo/\ http://www.isthe.com/chongo/ 24 * Share and enjoy! :-) http://www.isthe.com/chongo/tech/comp/calc/ 25 */ 26 27/* 28 * For a general tutorial on how to find a new largest known prime, see: 29 * 30 * http://www.isthe.com/chongo/tech/math/prime/prime-tutorial.pdf 31 * 32 * Also see the reference code, available both C and go: 33 * 34 * https://github.com/arcetri/goprime 35 */ 36 37/* 38 * NOTE: This is a standard calc resource file. For information on calc see: 39 * 40 * http://www.isthe.com/chongo/tech/comp/calc/index.html 41 * 42 * To obtain your own copy of calc, see: 43 * 44 * http://www.isthe.com/chongo/tech/comp/calc/calc-download.html 45 */ 46 47/* 48 * HISTORICAL NOTE: 49 * 50 * On 6 August 1989 at 00:53 PDT, the 'Amdahl 6', a team consisting of 51 * John Brown, Landon Curt Noll, Bodo Parady, Gene Smith, Joel Smith and 52 * Sergio Zarantonello proved the following 65087 digit number to be prime: 53 * 54 * 216193 55 * 391581 * 2 -1 56 * 57 * At the time of discovery, this number was the largest known prime. 58 * The primality was demonstrated by a program implementing the test 59 * found in these routines. An Amdahl 1200 takes 1987 seconds to test 60 * the primality of this number. A Cray 2 took several hours to 61 * confirm this prime. As of 31 Dec 1995, this prime was the 3rd 62 * largest known prime and the largest known non-Mersenne prime. 63 * 64 * The same team also discovered the following twin prime pair: 65 * 66 * 11235 11235 67 * 1706595 * 2 -1 1706595 * 2 +1 68 * 69 * At the time of discovery, this was the largest known twin prime pair. 70 * 71 * See: 72 * 73 * http://www.isthe.com/chongo/tech/math/prime/amdahl6.html 74 * 75 * for more information on the Amdahl 6 group. 76 * 77 * NOTE: Both largest known and largest known twin prime records have been 78 * broken. Rather than update this file each time, I'll just 79 * congratulate the finders and encourage others to try for 80 * larger finds. Records were made to be broken after all! 81 */ 82 83/* 84 * ON GAINING A WORLD RECORD: 85 * 86 * For a general tutorial on how to find a new largest known prime, see: 87 * 88 * http://www.isthe.com/chongo/tech/math/prime/prime-tutorial.pdf 89 * 90 * The routines in calc were designed to be portable, and to work on 91 * numbers of 'sane' size. The Amdahl 6 team used a 'ultra-high speed 92 * multi-precision' package that a machine dependent collection of routines 93 * tuned for a long trace vector processor to work with very large numbers. 94 * The heart of the package was a multiplication and square routine that 95 * was based on the PFA Fast Fourier Transform and on Winograd's radix FFTs. 96 * 97 * NOTE: While the PFA Fast Fourier Transform and Winograd's radix FFTs 98 * might have been optimal for the Amdahl 6 team at the time, 99 * they might not be optimal for your CPU architecture. See 100 * the above mentioned tutorial for information on better 101 * methods of performing multiplications and squares of very 102 * large numbers. 103 * 104 * Having a fast computer, and a good multi-precision package are 105 * critical, but one also needs to know where to look in order to have 106 * a good chance at a record. Knowing what to test is beyond the scope 107 * of this routine. However the following observations are noted: 108 * 109 * test numbers of the form h*2^n-1 110 * fix a value of n and vary the value h 111 * n mod 2^x == 0 for some value of x, say > 7 or more 112 * h*2^n-1 is not divisible by any small prime < 2^40 113 * 0 < h < 2^39 114 * h*2^n+1 is not divisible by any small prime < 2^40 115 * 116 * The Mersenne test for '2^n-1' is the fastest known primality test 117 * for a given large numbers. However, it is faster to search for 118 * primes of the form 'h*2^n-1'. When n is around 200000, one can find 119 * a prime of the form 'h*2^n-1' in about 1/2 the time. 120 * 121 * Critical to understanding why 'h*2^n-1' is to observe that primes of 122 * the form '2^n-1' seem to bunch around "islands". Such "islands" 123 * seem to be getting fewer and farther in-between, forcing the time 124 * for each test to grow longer and longer (worse then O(n^2 log n)). 125 * On the other hand, when one tests 'h*2^n-1', fixes 'n' and varies 126 * 'h', the time to test each number remains relatively constant. 127 * 128 * It is clearly a win to eliminate potential test candidates by 129 * rejecting numbers that that are divisible by 'small' primes. We 130 * (the "Amdahl 6") rejected all numbers that were divisible by primes 131 * less than '2^40'. We stopped looking for small factors at '2^40' 132 * when the rate of candidates being eliminated was slowed down to 133 * just a trickle. 134 * 135 * The 'n mod 128 == 0' restriction allows one to test for divisibility 136 * of small primes more quickly. To test of 'q' is a factor of 'k*2^n-1', 137 * one check to see if 'k*2^n mod q' == 1, which is the same a checking 138 * if 'h*(2^n mod q) mod q' == 1. One can compute '2^n mod q' by making 139 * use of the following: 140 * 141 * if 142 * y = 2^x mod q 143 * then 144 * 2^(2x) mod q == y^2 mod q 0 bit 145 * 2^(2x+1) mod q == 2*y^2 mod q 1 bit 146 * 147 * The choice of which expression depends on the binary pattern of 'n'. 148 * Since '1' bits require an extra step (multiply by 2), one should 149 * select value of 'n' that contain mostly '0' bits. The restriction 150 * of 'n mod 128 == 0' ensures that the bottom 7 bits of 'n' are 0. 151 * 152 * By limiting 'h' to '2^39' and eliminating all values divisible by 153 * small primes < twice the 'h' limit (2^40), one knows that all 154 * remaining candidates are relatively prime. Thus, when a candidate 155 * is proven to be composite (not prime) by the big test, one knows 156 * that the factors for that number (whatever they may be) will not 157 * be the factors of another candidate. 158 * 159 * Finally, one should eliminate all values of 'h*2^n-1' where 160 * 'h*2^n+1' is divisible by a small primes. 161 * 162 * NOTE: Today, for world record sized h*2^n-1 primes, one might 163 * search for factors < 2^46 or more. By excluding h*2^n-1 164 * with prime factors < 2^46, where h*2^n-1 is a bit larger 165 * than the largest known prime, one may exclude about 96.5% 166 * of candidates that have "small" prime factors. 167 */ 168 169pprod256 = 0; /* product of "primes up to 256" / "primes up to 46" */ 170 171/* 172 * lucas - lucas primality test on h*2^n-1 173 * 174 * ABOUT THE TEST: 175 * 176 * This routine will perform a primality test on h*2^n-1 based on 177 * the mathematics of Lucas, Lehmer and Riesel. One should read 178 * the following article: 179 * 180 * Ref1: 181 * "Lucasian Criteria for the Primality of N=h*2^n-1", by Hans Riesel, 182 * Mathematics of Computation, Vol 23 #108, pp. 869-875, Oct 1969 183 * 184 * http://www.ams.org/journals/mcom/1969-23-108/S0025-5718-1969-0262163-1/ 185 * S0025-5718-1969-0262163-1.pdf 186 * 187 * NOTE: Join the above two lines for the complete URL of the paper. 188 * 189 * The following book is also useful: 190 * 191 * Ref2: 192 * "Prime numbers and Computer Methods for Factorization", by Hans Riesel, 193 * Birkhauser, 1985, pp 131-134, 278-285, 438-444 194 * 195 * A few useful Legendre identities may be found in: 196 * 197 * Ref3: 198 * "Introduction to Analytic Number Theory", by Tom A. Apostol, 199 * Springer-Verlag, 1984, p 188. 200 * 201 * An excellent 5-page paper by Oystein J. Rodseth (we apologize that the 202 * ASCII character set does not allow us to spell his name with the 203 * umlaut marks on the O's): 204 * 205 * NOTE: The original Amdahl 6 method predates the publication of Ref4. 206 * The gen_v1() function used by lucas() uses the Ref4 method. 207 * See the 'Amdahl 6 legacy code' section below for the original 208 * method of generating v(1). 209 * 210 * Ref4: 211 * 212 * "A note on primality tests for N = h*2^n-1", by Oystein J. Rodseth, 213 * Department of Mathematics, University of Bergen, BIT Numerical 214 * Mathematics. 34 (3): pp 451-454. 215 * 216 * http://folk.uib.no/nmaoy/papers/luc.pdf 217 * 218 * This test is performed as follows: (see Ref1, Theorem 5) 219 * 220 * a) generate u(2) (see the function gen_u2() below) 221 * (NOTE: some call this u(0)) 222 * 223 * b) generate u(n) according to the rule: 224 * 225 * u(i+1) = u(i)^2-2 mod h*2^n-1 226 * 227 * c) h*2^n-1 is prime if and only if u(n) == 0 Q.E.D. :-) 228 * 229 * Now the following conditions must be true for the test to work: 230 * 231 * n >= 2 232 * h >= 1 233 * h < 2^n 234 * h mod 2 == 1 235 * 236 * A few miscellaneous notes: 237 * 238 * In order to reduce the number of tests, as attempt to eliminate 239 * any number that is divisible by a prime less than 257. Valid prime 240 * candidates less than 257 are declared prime as a special case. 241 * 242 * In real life, you would eliminate candidates by checking for 243 * divisibility by a prime much larger than 257 (perhaps as high 244 * as 2^39). 245 * 246 * The condition 'h mod 2 == 1' is not a problem. Say one is testing 247 * 'j*2^m-1', where j is even. If we note that: 248 * 249 * j mod 2^x == 0 for x>0 implies j*2^m-1 == ((j/2^x)*2^(m+x))-1, 250 * 251 * then we can let h=j/2^x and n=m+x and test 'h*2^n-1' which is the value. 252 * We need only consider odd values of h because we can rewrite our numbers 253 * do make this so. 254 * 255 * input: 256 * h h as in h*2^n-1 (must be >= 1) 257 * n n as in h*2^n-1 (must be >= 1) 258 * 259 * returns: 260 * 1 => h*2^n-1 is prime 261 * 0 => h*2^n-1 is not prime 262 * -1 => a test could not be formed, or h >= 2^n, h <= 0, n <= 0 263 */ 264define 265lucas(h, n) 266{ 267 local testval; /* h*2^n-1 */ 268 local shiftdown; /* the power of 2 that divides h */ 269 local u; /* the u(i) sequence value */ 270 local v1; /* the v(1) generator of u(2) */ 271 local i; /* u sequence cycle number */ 272 local oldh; /* pre-reduced h */ 273 local oldn; /* pre-reduced n */ 274 local bits; /* highbit of h*2^n-1 */ 275 276 /* 277 * check arg types 278 */ 279 if (!isint(h)) { 280 quit "FATAL: bad args: h must be an integer"; 281 } 282 if (h < 1) { 283 quit "FATAL: bad args: h must be an integer >= 1"; 284 } 285 if (!isint(n)) { 286 quit "FATAL: bad args: n must be an integer"; 287 } 288 if (n < 1) { 289 quit "FATAL: bad args: n must be an integer >= 1"; 290 } 291 292 /* 293 * reduce h if even 294 * 295 * we will force h to be odd by moving powers of two over to 2^n 296 */ 297 oldh = h; 298 oldn = n; 299 shiftdown = fcnt(h,2); /* h % 2^shiftdown == 0, max shiftdown */ 300 if (shiftdown > 0) { 301 h >>= shiftdown; 302 n += shiftdown; 303 } 304 305 /* 306 * enforce the 0 < h < 2^n rule 307 */ 308 if (h <= 0 || n <= 0) { 309 print "ERROR: reduced args violate the rule: 0 < h < 2^n"; 310 print " ERROR: h=":oldh, "n=":oldn, "reduced h=":h, "n=":n; 311 ldebug("lucas", "unknown: h <= 0 || n <= 0"); 312 return -1; 313 } 314 if (highbit(h) >= n) { 315 print "ERROR: reduced args violate the rule: h < 2^n"; 316 print " ERROR: h=":oldh, "n=":oldn, "reduced h=":h, "n=":n; 317 ldebug("lucas", "unknown: highbit(h) >= n"); 318 return -1; 319 } 320 321 /* 322 * catch the degenerate case of h*2^n-1 == 1 323 */ 324 if (h == 1 && n == 1) { 325 ldebug("lucas", "not prime: h == 1 && n == 1"); 326 return 0; /* 1*2^1-1 == 1 is not prime */ 327 } 328 329 /* 330 * catch the degenerate case of n==2 331 * 332 * n==2 and 0<h<2^n ==> 0<h<4 333 * 334 * Since h is now odd ==> h==1 or h==3 335 */ 336 if (h == 1 && n == 2) { 337 ldebug("lucas", "prime: h == 1 && n == 2"); 338 return 1; /* 1*2^2-1 == 3 is prime */ 339 } 340 if (h == 3 && n == 2) { 341 ldebug("lucas", "prime: h == 3 && n == 2"); 342 return 1; /* 3*2^2-1 == 11 is prime */ 343 } 344 345 /* 346 * catch small primes < 257 347 * 348 * We check for only a few primes because the other primes < 257 349 * violate the checks above. 350 */ 351 if (h == 1) { 352 if (n == 3 || n == 5 || n == 7) { 353 ldebug("lucas", "prime: 3, 7, 31, 127 are prime"); 354 return 1; /* 3, 7, 31, 127 are prime */ 355 } 356 } 357 if (h == 3) { 358 if (n == 2 || n == 3 || n == 4 || n == 6) { 359 ldebug("lucas", "prime: 11, 23, 47, 191 are prime"); 360 return 1; /* 11, 23, 47, 191 are prime */ 361 } 362 } 363 if (h == 5 && n == 4) { 364 ldebug("lucas", "prime: 79 is prime"); 365 return 1; /* 79 is prime */ 366 } 367 if (h == 7 && n == 5) { 368 ldebug("lucas", "prime: 223 is prime"); 369 return 1; /* 223 is prime */ 370 } 371 if (h == 15 && n == 4) { 372 ldebug("lucas", "prime: 239 is prime"); 373 return 1; /* 239 is prime */ 374 } 375 376 /* 377 * Verify that h*2^n-1 is not a multiple of 3 378 * 379 * The case for h*2^n-1 == 3 is handled above. 380 */ 381 if (((h % 3 == 1) && (n % 2 == 0)) || ((h % 3 == 2) && (n % 2 == 1))) { 382 /* no need to test h*2^n-1, it is a multiple of 3 */ 383 ldebug("lucas","not-prime: != 3 and is a multiple of 3"); 384 return 0; 385 } 386 387 /* 388 * Avoid any numbers divisible by small primes 389 */ 390 /* 391 * check for 5 <= prime factors < 31 392 * pfact(30)/6 = 1078282205 393 */ 394 testval = h*2^n - 1; 395 if (gcd(testval, 1078282205) > 1) { 396 /* a small 5 <= prime < 31 divides h*2^n-1 */ 397 ldebug("lucas",\ 398 "not-prime: a small 5<=prime<31 divides h*2^n-1"); 399 return 0; 400 } 401 /* 402 * check for 31 <= prime factors < 53 403 * pfact(52)/pfact(30) = 95041567 404 */ 405 if (gcd(testval, 95041567) > 1) { 406 /* a small 31 <= prime < 53 divides h*2^n-1 */ 407 ldebug("lucas","not-prime: 31<=prime<53 divides h*2^n-1"); 408 return 0; 409 } 410 /* 411 * check for prime 53 <= factors < 257, if h*2^n-1 is large 412 * 2^276 > pfact(256)/pfact(52) > 2^275 413 */ 414 bits = highbit(testval); 415 if (bits >= 275) { 416 if (pprod256 <= 0) { 417 pprod256 = pfact(256)/pfact(52); 418 } 419 if (gcd(testval, pprod256) > 1) { 420 /* a small 53 <= prime < 257 divides h*2^n-1 */ 421 ldebug("lucas",\ 422 "not-prime: 53<=prime<257 divides h*2^n-1"); 423 return 0; 424 } 425 } 426 427 /* 428 * try to compute u(2) (NOTE: some call this u(0)) 429 * 430 * We will use gen_v1() to give us a v(1) using the values 431 * of 'h' and 'n'. We will then use gen_u2() to convert 432 * the v(1) into u(2). 433 * 434 * If gen_v1() returns a negative value, then we failed to 435 * generate a test for h*2^n-1. The legacy function, 436 * legacy_gen_v1() used by the Amdahl 6 could have returned 437 * -1. The new gen_v1() based on the method outlined in Ref4 438 * will never return -1 if h*2^n-1 is not a multiple of 3. 439 * Because the "multiple of 3" case is handled above, the 440 * call below to gen_v1() will never return -1. 441 */ 442 v1 = gen_v1(h, n); 443 if (v1 < 0) { 444 /* failure to test number */ 445 print "unable to compute v(1) for", h : "*2^" : n : "-1"; 446 ldebug("lucas", "unknown: no v(1)"); 447 return -1; 448 } 449 u = gen_u2(h, n, v1); 450 451 /* 452 * compute u(n) (NOTE: some call this u(n-2)) 453 */ 454 for (i=3; i <= n; ++i) { 455 /* u = (u^2 - 2) % testval; */ 456 u = hnrmod(u^2 - 2, h, n, -1); 457 } 458 459 /* 460 * return 1 if prime, 0 is not prime 461 */ 462 if (u == 0) { 463 ldebug("lucas", "prime: end of test"); 464 return 1; 465 } else { 466 ldebug("lucas", "not-prime: end of test"); 467 return 0; 468 } 469} 470 471/* 472 * gen_u2 - determine the initial Lucas sequence for h*2^n-1 473 * 474 * Historically many start the Lucas sequence with u(0). 475 * Some, like the author of this code, prefer to start 476 * with U(2). This is so one may say: 477 * 478 * 2^p-1 is prime if u(p) = 0 mod 2^p-1 479 * or: 480 * h*2^p-1 is prime if u(p) = 0 mod h*2^p-1 481 * 482 * According to Ref1, Theorem 5: 483 * 484 * u(2) = alpha^h + alpha^(-h) (NOTE: Ref1 calls it u(0)) 485 * 486 * Now: 487 * 488 * v(x) = alpha^x + alpha^(-x) (Ref1, bottom of page 872) 489 * 490 * Therefore: 491 * 492 * u(2) = v(h) (NOTE: Ref1 calls it u(0)) 493 * 494 * We calculate v(h) as follows: (Ref1, top of page 873) 495 * 496 * v(0) = alpha^0 + alpha^(-0) = 2 497 * v(1) = alpha^1 + alpha^(-1) = gen_v1(h,n) 498 * v(n+2) = v(1)*v(n+1) - v(n) 499 * 500 * This function does not concern itself with the value of 'alpha'. 501 * The gen_v1() function is used to compute v(1), and identity 502 * functions take it from there. 503 * 504 * It can be shown that the following are true: 505 * 506 * v(2*n) = v(n)^2 - 2 507 * v(2*n+1) = v(n+1)*v(n) - v(1) 508 * 509 * To prevent v(x) from growing too large, one may replace v(x) with 510 * `v(x) mod h*2^n-1' at any time. 511 * 512 * See the function gen_v1() for details on the value of v(1). 513 * 514 * input: 515 * h - h as in h*2^n-1 (must be >= 1) 516 * n - n as in h*2^n-1 (must be >= 1) 517 * v1 - gen_v1(h,n) (must be >= 3) (see function below) 518 * 519 * returns: 520 * u(2) - initial value for Lucas test on h*2^n-1 521 * -1 - failed to generate u(2) 522 */ 523define 524gen_u2(h, n, v1) 525{ 526 local shiftdown; /* the power of 2 that divides h */ 527 local r; /* low value: v(n) */ 528 local s; /* high value: v(n+1) */ 529 local hbits; /* highest bit set in h */ 530 local oldh; /* pre-reduced h */ 531 local oldn; /* pre-reduced n */ 532 local i; 533 534 /* 535 * check arg types 536 */ 537 if (!isint(h)) { 538 quit "bad args: h must be an integer"; 539 } 540 if (h < 0) { 541 quit "bad args: h must be an integer >= 1"; 542 } 543 if (!isint(n)) { 544 quit "bad args: n must be an integer"; 545 } 546 if (n < 1) { 547 quit "bad args: n must be an integer >= 1"; 548 } 549 if (!isint(v1)) { 550 quit "bad args: v1 must be an integer"; 551 } 552 if (v1 < 3) { 553 quit "bogus arg: v1 must be an integer >= 3"; 554 } 555 556 /* 557 * reduce h if even 558 * 559 * we will force h to be odd by moving powers of two over to 2^n 560 */ 561 oldh = h; 562 oldn = n; 563 shiftdown = fcnt(h,2); /* h % 2^shiftdown == 0, max shiftdown */ 564 if (shiftdown > 0) { 565 h >>= shiftdown; 566 n += shiftdown; 567 } 568 569 /* 570 * enforce the h > 0 and n >= 2 rules 571 */ 572 if (h <= 0 || n < 1) { 573 print " ERROR: h=":oldh, "n=":oldn, "reduced h=":h, "n=":n; 574 quit "reduced args violate the rule: 0 < h < 2^n"; 575 } 576 hbits = highbit(h); 577 if (hbits >= n) { 578 print " ERROR: h=":oldh, "n=":oldn, "reduced h=":h, "n=":n; 579 quit "reduced args violate the rule: 0 < h < 2^n"; 580 } 581 582 /* 583 * build up u2 based on the reversed bits of h 584 */ 585 /* setup for bit loop */ 586 r = v1; 587 s = (r^2 - 2); 588 589 /* 590 * deal with small h as a special case 591 * 592 * The h value is odd > 0, and it needs to be 593 * at least 2 bits long for the loop below to work. 594 */ 595 if (h == 1) { 596 ldebug("gen_u2", "quick h == 1 case"); 597 /* return r%(h*2^n-1); */ 598 return hnrmod(r, h, n, -1); 599 } 600 601 /* cycle from second highest bit to second lowest bit of h */ 602 for (i=hbits-1; i > 0; --i) { 603 604 /* bit(i) is 1 */ 605 if (bit(h,i)) { 606 607 /* compute v(2n+1) = v(r+1)*v(r)-v1 */ 608 /* r = (r*s - v1) % (h*2^n-1); */ 609 r = hnrmod((r*s - v1), h, n, -1); 610 611 /* compute v(2n+2) = v(r+1)^2-2 */ 612 /* s = (s^2 - 2) % (h*2^n-1); */ 613 s = hnrmod((s^2 - 2), h, n, -1); 614 615 /* bit(i) is 0 */ 616 } else { 617 618 /* compute v(2n+1) = v(r+1)*v(r)-v1 */ 619 /* s = (r*s - v1) % (h*2^n-1); */ 620 s = hnrmod((r*s - v1), h, n, -1); 621 622 /* compute v(2n) = v(r)^-2 */ 623 /* r = (r^2 - 2) % (h*2^n-1); */ 624 r = hnrmod((r^2 - 2), h, n, -1); 625 } 626 } 627 628 /* we know that h is odd, so the final bit(0) is 1 */ 629 /* r = (r*s - v1) % (h*2^n-1); */ 630 r = hnrmod((r*s - v1), h, n, -1); 631 632 /* compute the final u2 return value */ 633 return r; 634} 635 636/* 637 * gen_u0 - determine the initial Lucas sequence for h*2^n-1 638 * 639 * Historically many start the Lucas sequence with u(0). 640 * Some, like the author of this code, prefer to start 641 * with u(2). This is so one may say: 642 * 643 * 2^p-1 is prime if u(p) = 0 mod 2^p-1 644 * or: 645 * h*2^n-1 is prime if U(n) = 0 mod h*2^n-1 646 * 647 * For those using the old code with gen_u0(), we 648 * simply call gen_u2() instead. 649 * 650 * See the function gen_u2() for details. 651 * 652 * input: 653 * h - h as in h*2^n-1 (must be >= 1) 654 * n - n as in h*2^n-1 (must be >= 1) 655 * v1 - gen_v1(h,n) (see function below) 656 * 657 * returns: 658 * u(2) - initial value for Lucas test on h*2^n-1 659 * -1 - failed to generate u(2) 660 */ 661define 662gen_u0(h, n, v1) 663{ 664 return gen_u2(h, n, v1); 665} 666 667/* 668 * rodseth_xhn - determine if v(1) == x for h*2^n-1 669 * 670 * For a given h*2^n-1, v(1) == x if: 671 * 672 * jacobi(x-2, h*2^n-1) == 1 (Ref4, condition 1) part 1 673 * jacobi(x+2, h*2^n-1) == -1 (Ref4, condition 1) part 2 674 * 675 * Now when x-2 <= 0: 676 * 677 * jacobi(x-2, h*2^n-1) == 0 678 * 679 * because: 680 * 681 * jacobi(x,y) == 0 if x <= 0 682 * 683 * So for (Ref4, condition 1) part 1 to be true: 684 * 685 * x-2 > 0 686 * 687 * And therefore: 688 * 689 * x > 2 690 * 691 * input: 692 * x potential v(1) value 693 * h h as in h*2^n-1 (h must be odd >= 1) 694 * n n as in h*2^n-1 (must be >= 1) 695 * 696 * returns: 697 * 1 if v(1) == x for h*2^n-1 698 * 0 otherwise 699 */ 700define 701rodseth_xhn(x, h, n) 702{ 703 local testval; /* h*2^n-1 */ 704 705 /* 706 * check arg types 707 */ 708 if (!isint(h)) { 709 quit "bad args: h must be an integer"; 710 } 711 if (iseven(h)) { 712 quit "bad args: h must be an odd integer"; 713 } 714 if (h < 1) { 715 quit "bad args: h must be an integer >= 1"; 716 } 717 if (!isint(n)) { 718 quit "bad args: n must be an integer"; 719 } 720 if (n < 1) { 721 quit "bad args: n must be an integer >= 1"; 722 } 723 if (!isint(x)) { 724 quit "bad args: x must be an integer"; 725 } 726 727 /* 728 * firewall 729 */ 730 if (x <= 2) { 731 return 0; 732 } 733 734 /* 735 * Check for jacobi(x-2, h*2^n-1) == 1 (Ref4, condition 1) part 1 736 */ 737 testval = h*2^n-1; 738 if (jacobi(x-2, testval) != 1) { 739 return 0; 740 } 741 742 /* 743 * Check for jacobi(x+2, h*2^n-1) == -1 (Ref4, condition 1) part 2 744 */ 745 if (jacobi(x+2, testval) != -1) { 746 return 0; 747 } 748 749 /* 750 * v(1) == x for this h*2^n-1 751 */ 752 return 1; 753} 754 755/* 756 * Trial tables used by gen_v1() 757 * 758 * When h mod 3 == 0, according to Ref4 we need to find the first value X where: 759 * 760 * jacobi(X-2, h*2^n-1) == 1 (Ref4, condition 1) part 1 761 * jacobi(X+2, h*2^n-1) == -1 (Ref4, condition 1) part 2 762 * 763 * We can show that X > 2. See the comments in the rodseth_xhn(x,h,n) above. 764 * 765 * Some values of X satisfy more often than others. For example a large sample 766 * of h*2^n-1, h odd multiple of 3, and large n (some around 1e4, some near 1e6, 767 * others near 3e7) where the sample size was 66 973 365, here is the count of 768 * the smallest value of X that satisfies conditions in Ref4, condition 1: 769 * 770 * count X 771 * ---------- 772 * 26791345 3 773 * 17223016 5 774 * 7829600 9 775 * 6988774 11 776 * 3301093 15 777 * 1517149 17 778 * 910346 21 779 * 711791 29 780 * 573403 20 781 * 390395 27 782 * 288637 35 783 * 149751 36 784 * 107733 39 785 * 58743 41 786 * 35619 45 787 * 25052 32 788 * 17775 51 789 * 13031 44 790 * 7563 56 791 * 7540 49 792 * 7060 59 793 * 4407 57 794 * 2948 65 795 * 2502 55 796 * 2388 69 797 * 2094 71 798 * 689 77 799 * 626 81 800 * 491 66 801 * 426 95 802 * 219 80 803 * 203 67 804 * 185 84 805 * 152 99 806 * 127 72 807 * 102 74 808 * 98 87 809 * 67 90 810 * 55 104 811 * 48 101 812 * 32 105 813 * 17 109 814 * 16 116 815 * 15 111 816 * 13 92 817 * 12 125 818 * 7 129 819 * 3 146 820 * 2 140 821 * 2 120 822 * 1 165 823 * 1 161 824 * 1 155 825 * 826 * It is important that we select the smallest possible v(1). While testing 827 * various values of X for V(1) is fast, using larger than necessary values 828 * of V(1) of can slow down calculating V(h). 829 * 830 * The above distribution was found to hold fairly well over many values of 831 * odd h that are also a multiple of 3 and for many values of n where h < 2^n. 832 * 833 * For example for in a sample size of 1000000 numbers of the form h*2^n-1 834 * where h is an odd multiple of 3, 12996351 <= h <= 13002351, 835 * 4331116 <= n <= 4332116, these are the smallest v(1) values that were found: 836 * 837 * smallest percentage 838 * v(1) used 839 * -------- --------- 840 * 3 40.0000 % 841 * 5 25.6833 % 842 * 9 11.6924 % 843 * 11 10.4528 % 844 * 15 4.8048 % 845 * 17 2.3458 % 846 * 21 1.3734 % 847 * 29 1.0527 % 848 * 20 0.8595 % 849 * 27 0.5758 % 850 * 35 0.4420 % 851 * 36 0.2433 % 852 * 39 0.1779 % 853 * 41 0.0885 % 854 * 45 0.0571 % 855 * 32 0.0337 % 856 * 51 0.0289 % 857 * 44 0.0205 % 858 * 49 0.0176 % 859 * 56 0.0137 % 860 * 59 0.0108 % 861 * 57 0.0053 % 862 * 65 0.0047 % 863 * 55 0.0045 % 864 * 69 0.0031 % 865 * 71 0.0024 % 866 * 66 0.0011 % 867 * 95 0.0008 % 868 * 81 0.0008 % 869 * 77 0.0006 % 870 * 72 0.0005 % 871 * 99 0.0004 % 872 * 80 0.0003 % 873 * 74 0.0003 % 874 * 84 0.0002 % 875 * 67 0.0002 % 876 * 87 0.0001 % 877 * 104 0.0001 % 878 * 129 0.0001 % 879 * 880 * However, a case can be made for considering only odd values for v(1) 881 * candidates. When h * 2^n-1 is prime and h is an odd multiple of 3, 882 * a smallest v(1) that is even is extremely rate. Of the list of 146553 883 * known primes of the form h*2^n-1 when h is an odd a multiple of 3, 884 * none has an smallest v(1) that was even. 885 * 886 * See: 887 * 888 * https://github.com/arcetri/verified-prime 889 * 890 * for that list of 146553 known primes of the form h*2^n-1. 891 * 892 * That same example for in a sample size of 1000000 numbers of the 893 * form h*2^n-1 where h is an odd multiple of 3, 12996351 <= h <= 13002351, 894 * 4331116 <= n <= 4332116, these are the smallest odd v(1) values that were 895 * found: 896 * 897 * smallest percentage 898 * odd v(1) used 899 * -------- --------- 900 * 3 40.0000 % 901 * 5 25.6833 % 902 * 9 11.6924 % 903 * 11 10.4528 % 904 * 15 4.8048 % 905 * 17 2.3458 % 906 * 21 1.6568 % 907 * 29 1.6174 % 908 * 35 0.4529 % 909 * 27 0.3546 % 910 * 39 0.3470 % 911 * 41 0.2159 % 912 * 45 0.1173 % 913 * 31 0.0661 % 914 * 51 0.0619 % 915 * 55 0.0419 % 916 * 59 0.0250 % 917 * 49 0.0170 % 918 * 69 0.0110 % 919 * 65 0.0098 % 920 * 71 0.0078 % 921 * 85 0.0048 % 922 * 81 0.0044 % 923 * 95 0.0038 % 924 * 99 0.0021 % 925 * 125 0.0009 % 926 * 57 0.0007 % 927 * 111 0.0005 % 928 * 77 0.0003 % 929 * 165 0.0003 % 930 * 155 0.0002 % 931 * 129 0.0002 % 932 * 101 0.0002 % 933 * 53 0.0001 % 934 * 935 * Moreover when evaluating odd candidates for v(1), one may cache Jacobi 936 * symbol evaluations to reduce the number of Jacobi symbol evaluations to 937 * a minimum. For example, if one tests 5 and finds that the 2nd case fails: 938 * 939 * jacobi(5+2, h*2^n-1) != -1 940 * 941 * Then if one is later testing 9, the Jacobi symbol value for the first 942 * 1st case: 943 * 944 * jacobi(7-2, h*2^n-1) 945 * 946 * is already known. 947 * 948 * Without Jacobi symbol value caching, it requires on average 949 * 4.851377 Jacobi symbol evaluations. With Jacobi symbol value caching 950 * cacheing, an average of 4.348820 Jacobi symbol evaluations is needed. 951 * 952 * Given this information, when odd h is a multiple of 3 we try, in order, 953 * these odd values of X: 954 * 955 * 3, 5, 9, 11, 15, 17, 21, 29, 27, 35, 39, 41, 31, 45, 51, 55, 49, 59, 956 * 69, 65, 71, 57, 85, 81, 95, 99, 77, 53, 67, 125, 111, 105, 87, 129, 957 * 101, 83, 165, 155, 149, 141, 121, 109 958 * 959 * And stop on the first value of X where: 960 * 961 * jacobi(X-2, h*2^n-1) == 1 962 * jacobi(X+2, h*2^n-1) == -1 963 * 964 * Less than 1 case out of 1000000 will not be satisfied by the above list. 965 * If no value in that list works, we start simple search starting with X = 167 966 * and incrementing by 2 until a value of X is found. 967 * 968 * The x_tbl[] matrix contains those values of X to try in order. 969 * If all x_tbl_len fail to satisfy Ref4 condition 1 (this happens less than 970 * 1 in 1000000 cases), then we begin a linear search of odd values starting at 971 * next_x until we find a proper X value. 972 */ 973x_tbl_len = 42; 974mat x_tbl[x_tbl_len]; 975x_tbl = { 976 3, 5, 9, 11, 15, 17, 21, 29, 27, 35, 39, 41, 31, 45, 51, 55, 49, 59, 977 69, 65, 71, 57, 85, 81, 95, 99, 77, 53, 67, 125, 111, 105, 87, 129, 978 101, 83, 165, 155, 149, 141, 121, 109 979}; 980next_x = 167; /* must be 2 more than the largest value in x_tbl[] */ 981 982/* 983 * gen_v1 - compute the v(1) for a given h*2^n-1 if we can 984 * 985 * This function assumes: 986 * 987 * n > 2 (n==2 has already been eliminated) 988 * h mod 2 == 1 989 * h < 2^n 990 * h*2^n-1 mod 3 != 0 (h*2^n-1 has no small factors, such as 3) 991 * 992 * The generation of v(1) depends on the value of h. There are two cases 993 * to consider, h mod 3 != 0, and h mod 3 == 0. 994 * 995 *** 996 * 997 * Case 1: (h mod 3 != 0) 998 * 999 * This case is easy. 1000 * 1001 * In Ref1, page 869, one finds that if: (or see Ref2, page 131-132) 1002 * 1003 * h mod 6 == +/-1 1004 * h*2^n-1 mod 3 != 0 1005 * 1006 * which translates, gives the functions assumptions, into the condition: 1007 * 1008 * h mod 3 != 0 1009 * 1010 * If this case condition is true, then: 1011 * 1012 * u(2) = (2+sqrt(3))^h + (2-sqrt(3))^h (see Ref1, page 869) 1013 * = (2+sqrt(3))^h + (2+sqrt(3))^(-h) (NOTE: some call this u(2)) 1014 * 1015 * and since Ref1, Theorem 5 states: 1016 * 1017 * u(2) = alpha^h + alpha^(-h) (NOTE: some call this u(2)) 1018 * r = abs(2^2 - 1^2*3) = 1 1019 * 1020 * where these values work for Case 1: (h mod 3 != 0) 1021 * 1022 * a = 1 1023 * b = 2 1024 * D = 1 1025 * 1026 * Now at the bottom of Ref1, page 872 states: 1027 * 1028 * v(x) = alpha^x + alpha^(-x) 1029 * 1030 * If we let: 1031 * 1032 * alpha = (2+sqrt(3)) 1033 * 1034 * then 1035 * 1036 * u(2) = v(h) (NOTE: some call this u(2)) 1037 * 1038 * so we can always return 1039 * 1040 * v(1) = alpha^1 + alpha^(-1) 1041 * = (2+sqrt(3)) + (2-sqrt(3)) 1042 * = 4 1043 * 1044 * In 40% of the cases when h is not a multiple of 3, 3 is a valid value 1045 * for v(1). We can test if 3 is a valid value for v(1) in this case: 1046 * 1047 * if jacobi(1, h*2^n-1) == 1 and jacobi(5, h*2^n-1) == -1, then 1048 * v(1) = 3 1049 * else 1050 * v(1) = 4 1051 * 1052 * NOTE: The above "if then else" works only of h is not a multiple of 3. 1053 * 1054 *** 1055 * 1056 * Case 2: (h mod 3 == 0) 1057 * 1058 * For the case where h is a multiple of 3, we turn to Ref4. 1059 * 1060 * The central theorem on page 3 of that paper states that 1061 * we may set v(1) to the first value X that satisfies: 1062 * 1063 * jacobi(X-2, h*2^n-1) == 1 (Ref4, condition 1) 1064 * jacobi(X+2, h*2^n-1) == -1 (Ref4, condition 1) 1065 * 1066 * NOTE: Ref4 uses P, which we shall refer to as X. 1067 * Ref4 uses N, which we shall refer to as h*2^n-1. 1068 * 1069 * NOTE: Ref4 uses the term Legendre-Jacobi symbol, which 1070 * we shall refer to as the Jacobi symbol. 1071 * 1072 * Before we address the two conditions, we need some background information 1073 * on two symbols, Legendre and Jacobi. In Ref 2, pp 278, 284-285, we find 1074 * the following definitions of jacobi(a,b) and L(a,p): 1075 * 1076 * The Legendre symbol L(a,p) takes the value: 1077 * 1078 * L(a,p) == 1 => a is a quadratic residue of p 1079 * L(a,p) == -1 => a is NOT a quadratic residue of p 1080 * 1081 * when: 1082 * 1083 * p is prime 1084 * p mod 2 == 1 1085 * gcd(a,p) == 1 1086 * 1087 * The value a is a quadratic residue of b if there exists some integer z 1088 * such that: 1089 * 1090 * z^2 mod b == a 1091 * 1092 * The Jacobi symbol jacobi(a,b) takes the value: 1093 * 1094 * jacobi(a,b) == 1 => b is not prime, 1095 * or a is a quadratic residue of b 1096 * jacobi(a,b) == -1 => a is NOT a quadratic residue of b 1097 * 1098 * when 1099 * 1100 * b mod 2 == 1 1101 * gcd(a,b) == 1 1102 * 1103 * It is worth noting for the Legendre symbol, in order for L(X+/-2, 1104 * h*2^n-1) to be defined, we must ensure that neither X-2 nor X+2 are 1105 * factors of h*2^n-1. This is done by pre-screening h*2^n-1 to not 1106 * have small factors and keeping X+2 less than that small factor 1107 * limit. It is worth noting that in lucas(h, n), we first verify 1108 * that h*2^n-1 does not have a factor < 257 before performing the 1109 * primality test. So while X+/-2 < 257, we know that 1110 * gcd(X+/-2, h*2^n-1) == 1. 1111 * 1112 * Returning to the testing of conditions in Ref4, condition 1: 1113 * 1114 * jacobi(X-2, h*2^n-1) == 1 1115 * jacobi(X+2, h*2^n-1) == -1 1116 * 1117 * When such an X is found, we set: 1118 * 1119 * v(1) = X 1120 * 1121 *** 1122 * 1123 * In conclusion, we can compute v,(1) by attempting to do the following: 1124 * 1125 * h mod 3 != 0 1126 * 1127 * we return: 1128 * 1129 * v(1) == 4 1130 * 1131 * h mod 3 == 0 1132 * 1133 * we return: 1134 * 1135 * v(1) = X 1136 * 1137 * where X > 2 in a integer such that: 1138 * 1139 * jacobi(X-2, h*2^n-1) == 1 1140 * jacobi(X+2, h*2^n-1) == -1 1141 * 1142 *** 1143 * 1144 * input: 1145 * h h as in h*2^n-1 (h must be odd >= 1) 1146 * n n as in h*2^n-1 (must be >= 1) 1147 * 1148 * output: 1149 * returns v(1), or 1150 * -1 when h*2^n-1 is a multiple of 3 1151 */ 1152define 1153gen_v1(h, n) 1154{ 1155 local x; /* potential v(1) to test */ 1156 local i; /* x_tbl index */ 1157 local v1m2; /* X-2 1st case */ 1158 local v1p2; /* X+2 2nd case */ 1159 local testval; /* h*2^n-1 - value we are testing if prime */ 1160 local mat cached_v1[next_x]; /* cached Jacobi symbol values or 0 */ 1161 1162 /* 1163 * check arg types 1164 */ 1165 if (!isint(h)) { 1166 quit "bad args: h must be an integer"; 1167 } 1168 if (iseven(h)) { 1169 quit "bad args: h must be an odd integer"; 1170 } 1171 if (h < 1) { 1172 quit "bad args: h must be an integer >= 1"; 1173 } 1174 if (!isint(n)) { 1175 quit "bad args: n must be an integer"; 1176 } 1177 if (n < 1) { 1178 quit "bad args: n must be an integer >= 1"; 1179 } 1180 1181 /* 1182 * pretest: Verify that h*2^n-1 is not a multiple of 3 1183 */ 1184 if (((h % 3 == 1) && (n % 2 == 0)) || ((h % 3 == 2) && (n % 2 == 1))) { 1185 /* no need to test h*2^n-1, it is not prime */ 1186 return -1; 1187 } 1188 1189 /* 1190 * Common Mersenne number case: 1191 * 1192 * For Mersenne numbers: 1193 * 1194 * 2^n-1 1195 * 1196 * we can use, 40% of the time, v(1) == 3. However nearly all code that 1197 * implements the Lucas-Lehmer test uses v(1) == 4. Whenever for 1198 * h != 0 mod 3, and particular the Mersenne number case of when h == 1: 1199 * 1200 * 1*2^n-1 1201 * 1202 * v(1) == 4 always works. For this reason, we return 4 when h == 1. 1203 */ 1204 if (h == 1) { 1205 /* v(1) == 4 always works for the Mersenne number case */ 1206 return 4; 1207 } 1208 1209 /* 1210 * check for Case 1: (h mod 3 != 0) 1211 */ 1212 if (h % 3 != 0) { 1213 if (rodseth_xhn(3, h, n) == 1) { 1214 /* 40% of the time, 3 works when h mod 3 != 0 */ 1215 return 3; 1216 } else { 1217 /* otherwise 4 always works when h mod 3 != 0 */ 1218 return 4; 1219 } 1220 } 1221 1222 /* 1223 * What follow is Case 2: (h mod 3 == 0) 1224 */ 1225 1226 /* 1227 * clear cache 1228 */ 1229 matfill(cached_v1, 0); 1230 1231 /* 1232 * We will look for x that satisfies conditions in Ref4, condition 1: 1233 * 1234 * jacobi(X-2, h*2^n-1) == 1 part 1 1235 * jacobi(X+2, h*2^n-1) == -1 part 2 1236 * 1237 * NOTE: If we wanted to be super optimal, we would cache 1238 * jacobi(X+2, h*2^n-1) that that when we increment X 1239 * to the next odd value, the now jacobi(X-2, h*2^n-1) 1240 * does not need to be re-evaluated. 1241 */ 1242 testval = h*2^n-1; 1243 for (i=0; i < x_tbl_len; ++i) { 1244 1245 /* 1246 * obtain the next test candidate 1247 */ 1248 x = x_tbl[i]; 1249 1250 /* 1251 * Check x for condition 1 part 1 1252 * 1253 * jacobi(x-2, h*2^n-1) == 1 1254 */ 1255 v1m2 = x-2; 1256 if (cached_v1[v1m2] == 0) { 1257 cached_v1[v1m2] = jacobi(v1m2, testval); 1258 } 1259 if (cached_v1[v1m2] != 1) { 1260 continue; 1261 } 1262 1263 /* 1264 * Check x for condition 1 part 2 1265 * 1266 * jacobi(x+2, h*2^n-1) == -1 1267 */ 1268 v1p2 = x+2; 1269 if (cached_v1[v1p2] == 0) { 1270 cached_v1[v1p2] = jacobi(v1p2, testval); 1271 } 1272 if (cached_v1[v1p2] != -1) { 1273 continue; 1274 } 1275 1276 /* 1277 * found a x that satisfies Ref4 condition 1 1278 */ 1279 ldebug("gen_v1", "h= " + str(h) + " n= " + str(n) + 1280 " v1= " + str(x) + " using tbl[ " + 1281 str(i) + " ]"); 1282 return x; 1283 } 1284 1285 /* 1286 * We are in that rare case (less than 1 in 1 000 000) where none of the 1287 * common X values satisfy Ref4 condition 1. We start a linear search 1288 * of odd values at next_x from here on. 1289 */ 1290 x = next_x; 1291 while (rodseth_xhn(x, h, n) != 1) { 1292 x += 2; 1293 } 1294 /* finally found a v(1) value */ 1295 ldebug("gen_v1", "h= " + str(h) + " n= " + str(n) + 1296 " v1= " + str(x) + " beyond tbl"); 1297 return x; 1298} 1299 1300/* 1301 * ldebug - print a debug statement 1302 * 1303 * input: 1304 * funct name of calling function 1305 * str string to print 1306 */ 1307define 1308ldebug(funct, str) 1309{ 1310 if (config("resource_debug") & 8) { 1311 print "DEBUG:", funct:":", str; 1312 } 1313 return; 1314} 1315 1316/* 1317 ************************ 1318 * Amdahl 6 legacy code * 1319 ************************ 1320 * 1321 * NOTE: What follows is legacy code based on the method used by the 1322 * Amdahl 6 group: 1323 * 1324 * John Brown, Landon Curt Noll, Bodo Parady, Gene Smith, 1325 * Joel Smith and Sergio Zarantonello 1326 * 1327 * This method generated v(1) for nearly all values, except for a 1328 * few rare cases when h mod 3 == 0. The code is NOT used by lucas.cal 1329 * above. The gen_v1() function above is based on an improved method 1330 * outlined in Ref4. That method generated v(1) for all h. 1331 * 1332 * The code below is kept for historical purposes only. The functions 1333 * and global variables of the Amdahl 6 legacy code all begin with legacy_. 1334 */ 1335 1336/* 1337 * Trial tables used by legacy_gen_v1() 1338 * 1339 * When h mod 3 == 0, one needs particular values of D, a and b (see 1340 * legacy_gen_v1 documentation) in order to find a value of v(1). 1341 * 1342 * This table defines 'legacy_quickmax' possible tests to be taken in ascending 1343 * order. The legacy_v1_qval[x] refers to a v(1) value from Ref1, Table 1. A 1344 * related D value is found in legacy_d_qval[x]. All D values expect 1345 * legacy_d_qval[1] are also taken from Ref1, Table 1. The case of D == 21 as 1346 * listed in Ref1, Table 1 can be changed to D == 7 for the sake of the test 1347 * because of {note 6}. 1348 * 1349 * It should be noted that the D values all satisfy the selection values 1350 * as outlined in the legacy_gen_v1() function comments. That is: 1351 * 1352 * D == P*(2^f)*(3^g) 1353 * 1354 * where f == 0 and g == 0, P == D. So we simply need to check that 1355 * one of the following two cases are true: 1356 * 1357 * P mod 4 == 1 and J(h*2^n-1 mod P, P) == -1 1358 * P mod 4 == -1 and J(h*2^n-1 mod P, P) == 1 1359 * 1360 * In all cases, the value of r is: 1361 * 1362 * r == Q*(2^j)*(3^k)*(z^2) 1363 * 1364 * where Q == 1. No further processing is needed to compute v(1) when r 1365 * is of this form. 1366 */ 1367legacy_quickmax = 8; 1368mat legacy_d_qval[legacy_quickmax]; 1369mat legacy_v1_qval[legacy_quickmax]; 1370legacy_d_qval[0] = 5; legacy_v1_qval[0] = 3; /* a=1 b=1 r=4 */ 1371legacy_d_qval[1] = 7; legacy_v1_qval[1] = 5; /* a=3 b=1 r=12 D=21 */ 1372legacy_d_qval[2] = 13; legacy_v1_qval[2] = 11; /* a=3 b=1 r=4 */ 1373legacy_d_qval[3] = 11; legacy_v1_qval[3] = 20; /* a=3 b=1 r=2 */ 1374legacy_d_qval[4] = 29; legacy_v1_qval[4] = 27; /* a=5 b=1 r=4 */ 1375legacy_d_qval[5] = 53; legacy_v1_qval[5] = 51; /* a=53 b=1 r=4 */ 1376legacy_d_qval[6] = 17; legacy_v1_qval[6] = 66; /* a=17 b=1 r=1 */ 1377legacy_d_qval[7] = 19; legacy_v1_qval[7] = 74; /* a=38 b=1 r=2 */ 1378 1379/* 1380 * legacy_gen_v1 - compute the v(1) for a given h*2^n-1 if we can 1381 * 1382 * This function assumes: 1383 * 1384 * n > 2 (n==2 has already been eliminated) 1385 * h mod 2 == 1 1386 * h < 2^n 1387 * h*2^n-1 mod 3 != 0 (h*2^n-1 has no small factors, such as 3) 1388 * 1389 * The generation of v(1) depends on the value of h. There are two cases 1390 * to consider, h mod 3 != 0, and h mod 3 == 0. 1391 * 1392 *** 1393 * 1394 * Case 1: (h mod 3 != 0) 1395 * 1396 * This case is easy and always finds v(1). 1397 * 1398 * In Ref1, page 869, one finds that if: (or see Ref2, page 131-132) 1399 * 1400 * h mod 6 == +/-1 1401 * h*2^n-1 mod 3 != 0 1402 * 1403 * which translates, gives the functions assumptions, into the condition: 1404 * 1405 * h mod 3 != 0 1406 * 1407 * If this case condition is true, then: 1408 * 1409 * u(2) = (2+sqrt(3))^h + (2-sqrt(3))^h (see Ref1, page 869) 1410 * = (2+sqrt(3))^h + (2+sqrt(3))^(-h) (some call this u(0)) 1411 * 1412 * and since Ref1, Theorem 5 states: 1413 * 1414 * u(2) = alpha^h + alpha^(-h) 1415 * r = abs(2^2 - 1^2*3) = 1 1416 * 1417 * where these values work for Case 1: (h mod 3 != 0) 1418 * 1419 * a = 1 1420 * b = 2 1421 * D = 1 1422 * 1423 * Now at the bottom of Ref1, page 872 states: 1424 * 1425 * v(x) = alpha^x + alpha^(-x) 1426 * 1427 * If we let: 1428 * 1429 * alpha = (2+sqrt(3)) 1430 * 1431 * then 1432 * 1433 * u(2) = v(h) 1434 * 1435 * so we simply return 1436 * 1437 * v(1) = alpha^1 + alpha^(-1) 1438 * = (2+sqrt(3)) + (2-sqrt(3)) 1439 * = 4 1440 * 1441 *** 1442 * 1443 * Case 2: (h mod 3 == 0) 1444 * 1445 * This case is not so easy and finds v(1) in most all cases. In this 1446 * version of this program, we will simply return -1 (failure) if we 1447 * hit one of the cases that fall thru the cracks. This does not happen 1448 * often, so this is not too bad. 1449 * 1450 * Ref1, Theorem 5 contains the following definitions: 1451 * 1452 * r = abs(a^2 - b^2*D) 1453 * alpha = (a + b*sqrt(D))^2/r 1454 * 1455 * where D is 'square free', and 'alpha = epsilon^s' (for some s>0) are units 1456 * in the quadratic field K(sqrt(D)). 1457 * 1458 * One can find possible values for a, b and D in Ref1, Table 1 (page 872). 1459 * (see the file lucas_tbl.cal) 1460 * 1461 * Now Ref1, Theorem 5 states that if: 1462 * 1463 * L(D, h*2^n-1) = -1 [condition 1] 1464 * L(r, h*2^n-1) * (a^2 - b^2*D)/r = -1 [condition 2] 1465 * 1466 * where L(x,y) is the Legendre symbol (see below), then: 1467 * 1468 * u(2) = alpha^h + alpha^(-h) 1469 * 1470 * The bottom of Ref1, page 872 states: 1471 * 1472 * v(x) = alpha^x + alpha^(-x) 1473 * 1474 * thus since: 1475 * 1476 * u(2) = v(h) 1477 * 1478 * so we want to return: 1479 * 1480 * v(1) = alpha^1 + alpha^(-1) 1481 * 1482 * Therefore we need to take a given (D,a,b), determine if the two conditions 1483 * are true, and return the related v(1). 1484 * 1485 * Before we address the two conditions, we need some background information 1486 * on two symbols, Legendre and Jacobi. In Ref 2, pp 278, 284-285, we find 1487 * the following definitions of J(a,p) and L(a,n): 1488 * 1489 * The Legendre symbol L(a,p) takes the value: 1490 * 1491 * L(a,p) == 1 => a is a quadratic residue of p 1492 * L(a,p) == -1 => a is NOT a quadratic residue of p 1493 * 1494 * when 1495 * 1496 * p is prime 1497 * p mod 2 == 1 1498 * gcd(a,p) == 1 1499 * 1500 * The value x is a quadratic residue of y if there exists some integer z 1501 * such that: 1502 * 1503 * z^2 mod y == x 1504 * 1505 * The Jacobi symbol J(x,y) takes the value: 1506 * 1507 * J(x,y) == 1 => y is not prime, or x is a quadratic residue of y 1508 * J(x,y) == -1 => x is NOT a quadratic residue of y 1509 * 1510 * when 1511 * 1512 * y mod 2 == 1 1513 * gcd(x,y) == 1 1514 * 1515 * In the following comments on Legendre and Jacobi identities, we shall 1516 * assume that the arguments to the symbolic are valid over the symbol 1517 * definitions as stated above. 1518 * 1519 * In Ref2, pp 280-284, we find that: 1520 * 1521 * L(a,p)*L(b,p) == L(a*b,p) {A3.5} 1522 * J(x,y)*J(z,y) == J(x*z,y) {A3.14} 1523 * L(a,p) == L(p,a) * (-1)^((a-1)*(p-1)/4) {A3.8} 1524 * J(x,y) == J(y,x) * (-1)^((x-1)*(y-1)/4) {A3.17} 1525 * 1526 * The equality L(a,p) == J(a,p) when: {note 0} 1527 * 1528 * p is prime 1529 * p mod 2 == 1 1530 * gcd(a,p) == 1 1531 * 1532 * It can be shown that (see Ref3): 1533 * 1534 * L(a,p) == L(a mod p, p) {note 1} 1535 * L(z^2, p) == 1 {note 2} 1536 * 1537 * From Ref2, table 32: 1538 * 1539 * p mod 8 == +/-1 implies L(2,p) == 1 {note 3} 1540 * p mod 12 == +/-1 implies L(3,p) == 1 {note 4} 1541 * 1542 * Since h*2^n-1 mod 8 == -1, for n>2, note 3 implies: 1543 * 1544 * L(2, h*2^n-1) == 1 (n>2) {note 5} 1545 * 1546 * Since h=3*A, h*2^n-1 mod 12 == -1, for A>0, note 4 implies: 1547 * 1548 * L(3, h*2^n-1) == 1 {note 6} 1549 * 1550 * By use of {A3.5}, {note 2}, {note 5} and {note 6}, one can show: 1551 * 1552 * L((2^g)*(3^l)*(z^2), h*2^n-1) == 1 (g>=0,l>=0,z>0,n>2) {note 7} 1553 * 1554 * Returning to the testing of conditions, take condition 1: 1555 * 1556 * L(D, h*2^n-1) == -1 [condition 1] 1557 * 1558 * In order for J(D, h*2^n-1) to be defined, we must ensure that D 1559 * is not a factor of h*2^n-1. This is done by pre-screening h*2^n-1 to 1560 * not have small factors and selecting D less than that factor check limit. 1561 * 1562 * By use of {note 7}, we can show that when we choose D to be: 1563 * 1564 * D is square free 1565 * D = P*(2^f)*(3^g) (P is prime>2) 1566 * 1567 * The square free condition implies f = 0 or 1, g = 0 or 1. If f and g 1568 * are both 1, P must be a prime > 3. 1569 * 1570 * So given such a D value: 1571 * 1572 * L(D, h*2^n-1) == L(P*(2^g)*(3^l), h*2^n-1) 1573 * == L(P, h*2^n-1) * L((2^g)*(3^l), h*2^n-1) {A3.5} 1574 * == L(P, h*2^n-1) * 1 {note 7} 1575 * == L(h*2^n-1, P)*(-1)^((h*2^n-2)*(P-1)/4) {A3.8} 1576 * == L(h*2^n-1 mod P, P)*(-1)^((h*2^n-2)*(P-1)/4) {note 1} 1577 * == J(h*2^n-1 mod P, P)*(-1)^((h*2^n-2)*(P-1)/4) {note 0} 1578 * 1579 * When does J(h*2^n-1 mod P, P)*(-1)^((h*2^n-2)*(P-1)/4) take the value of -1, 1580 * thus satisfy [condition 1]? The answer depends on P. Now P is a prime>2, 1581 * thus P mod 4 == 1 or -1. 1582 * 1583 * Take P mod 4 == 1: 1584 * 1585 * P mod 4 == 1 implies (-1)^((h*2^n-2)*(P-1)/4) == 1 1586 * 1587 * Thus: 1588 * 1589 * L(D, h*2^n-1) == L(h*2^n-1 mod P, P) * (-1)^((h*2^n-2)*(P-1)/4) 1590 * == L(h*2^n-1 mod P, P) 1591 * == J(h*2^n-1 mod P, P) 1592 * 1593 * Take P mod 4 == -1: 1594 * 1595 * P mod 4 == -1 implies (-1)^((h*2^n-2)*(P-1)/4) == -1 1596 * 1597 * Thus: 1598 * 1599 * L(D, h*2^n-1) == L(h*2^n-1 mod P, P) * (-1)^((h*2^n-2)*(P-1)/4) 1600 * == L(h*2^n-1 mod P, P) * -1 1601 * == -J(h*2^n-1 mod P, P) 1602 * 1603 * Therefore [condition 1] is met if, and only if, one of the following 1604 * to cases are true: 1605 * 1606 * P mod 4 == 1 and J(h*2^n-1 mod P, P) == -1 1607 * P mod 4 == -1 and J(h*2^n-1 mod P, P) == 1 1608 * 1609 * Now consider [condition 2]: 1610 * 1611 * L(r, h*2^n-1) * (a^2 - b^2*D)/r == -1 [condition 2] 1612 * 1613 * We select only a, b, r and D values where: 1614 * 1615 * (a^2 - b^2*D)/r == -1 1616 * 1617 * Therefore in order for [condition 2] to be met, we must show that: 1618 * 1619 * L(r, h*2^n-1) == 1 1620 * 1621 * If we select r to be of the form: 1622 * 1623 * r == Q*(2^j)*(3^k)*(z^2) (Q == 1, j>=0, k>=0, z>0) 1624 * 1625 * then by use of {note 7}: 1626 * 1627 * L(r, h*2^n-1) == L(Q*(2^j)*(3^k)*(z^2), h*2^n-1) 1628 * == L((2^j)*(3^k)*(z^2), h*2^n-1) 1629 * == 1 {note 2} 1630 * 1631 * and thus, [condition 2] is met. 1632 * 1633 * If we select r to be of the form: 1634 * 1635 * r == Q*(2^j)*(3^k)*(z^2) (Q is prime>2, j>=0, k>=0, z>0) 1636 * 1637 * then by use of {note 7}: 1638 * 1639 * L(r, h*2^n-1) == L(Q*(2^j)*(3^k)*(z^2), h*2^n-1) 1640 * == L(Q, h*2^n-1) * L((2^j)*(3^k)*(z^2), h*2^n-1) {A3.5} 1641 * == L(Q, h*2^n-1) * 1 {note 2} 1642 * == L(h*2^n-1, Q) * (-1)^((h*2^n-2)*(Q-1)/4) {A3.8} 1643 * == L(h*2^n-1 mod Q, Q)*(-1)^((h*2^n-2)*(Q-1)/4) {note 1} 1644 * == J(h*2^n-1 mod Q, Q)*(-1)^((h*2^n-2)*(Q-1)/4) {note 0} 1645 * 1646 * When does J(h*2^n-1 mod Q, Q)*(-1)^((h*2^n-2)*(Q-1)/4) take the value of 1, 1647 * thus satisfy [condition 2]? The answer depends on Q. Now Q is a prime>2, 1648 * thus Q mod 4 == 1 or -1. 1649 * 1650 * Take Q mod 4 == 1: 1651 * 1652 * Q mod 4 == 1 implies (-1)^((h*2^n-2)*(Q-1)/4) == 1 1653 * 1654 * Thus: 1655 * 1656 * L(D, h*2^n-1) == L(h*2^n-1 mod Q, Q) * (-1)^((h*2^n-2)*(Q-1)/4) 1657 * == L(h*2^n-1 mod Q, Q) 1658 * == J(h*2^n-1 mod Q, Q) 1659 * 1660 * Take Q mod 4 == -1: 1661 * 1662 * Q mod 4 == -1 implies (-1)^((h*2^n-2)*(Q-1)/4) == -1 1663 * 1664 * Thus: 1665 * 1666 * L(D, h*2^n-1) == L(h*2^n-1 mod Q, Q) * (-1)^((h*2^n-2)*(Q-1)/4) 1667 * == L(h*2^n-1 mod Q, Q) * -1 1668 * == -J(h*2^n-1 mod Q, Q) 1669 * 1670 * Therefore [condition 2] is met by selecting D = Q*(2^j)*(3^k)*(z^2), 1671 * where Q is prime>2, j>=0, k>=0, z>0; if and only if one of the following 1672 * to cases are true: 1673 * 1674 * Q mod 4 == 1 and J(h*2^n-1 mod Q, Q) == 1 1675 * Q mod 4 == -1 and J(h*2^n-1 mod Q, Q) == -1 1676 * 1677 *** 1678 * 1679 * In conclusion, we can compute v(1) by attempting to do the following: 1680 * 1681 * h mod 3 != 0 1682 * 1683 * we return: 1684 * 1685 * v(1) == 4 1686 * 1687 * h mod 3 == 0 1688 * 1689 * define: 1690 * 1691 * r == abs(a^2 - b^2*D) 1692 * alpha == (a + b*sqrt(D))^2/r 1693 * 1694 * we return: 1695 * 1696 * v(1) = alpha^1 + alpha^(-1) 1697 * 1698 * if and only if we can find a given a, b, D that obey all the 1699 * following selection rules: 1700 * 1701 * D is square free 1702 * 1703 * D == P*(2^f)*(3^g) (P is prime>2, f,g == 0 or 1) 1704 * 1705 * (a^2 - b^2*D)/r == -1 1706 * 1707 * r == Q*(2^j)*(3^k)*(z^2) (Q==1 or Q is prime>2, j>=0, k>=0, z>0) 1708 * 1709 * one of the following is true: 1710 * P mod 4 == 1 and J(h*2^n-1 mod P, P) == -1 1711 * P mod 4 == -1 and J(h*2^n-1 mod P, P) == 1 1712 * 1713 * if Q is prime, then one of the following is true: 1714 * Q mod 4 == 1 and J(h*2^n-1 mod Q, Q) == 1 1715 * Q mod 4 == -1 and J(h*2^n-1 mod Q, Q) == -1 1716 * 1717 * If we cannot find a v(1) quickly enough, then we will give up 1718 * testing h*2^n-1. This does not happen too often, so this hack 1719 * is not too bad. 1720 * 1721 *** 1722 * 1723 * input: 1724 * h h as in h*2^n-1 (must be >= 1) 1725 * n n as in h*2^n-1 (must be >= 1) 1726 * 1727 * output: 1728 * returns v(1), or -1 is there is no quick way 1729 */ 1730define 1731legacy_gen_v1(h, n) 1732{ 1733 local d; /* the 'D' value to try */ 1734 local val_mod; /* h*2^n-1 mod 'D' */ 1735 local i; 1736 1737 /* 1738 * check arg types 1739 */ 1740 if (!isint(h)) { 1741 quit "bad args: h must be an integer"; 1742 } 1743 if (h < 1) { 1744 quit "bad args: h must be an integer >= 1"; 1745 } 1746 if (!isint(n)) { 1747 quit "bad args: n must be an integer"; 1748 } 1749 if (n < 1) { 1750 quit "bad args: n must be an integer >= 1"; 1751 } 1752 1753 /* 1754 * check for case 1 1755 */ 1756 if (h % 3 != 0) { 1757 /* v(1) is easy to compute */ 1758 return 4; 1759 } 1760 1761 /* 1762 * We will try all 'D' values until we find a proper v(1) 1763 * or run out of 'D' values. 1764 */ 1765 for (i=0; i < legacy_quickmax; ++i) { 1766 1767 /* grab our 'D' value */ 1768 d = legacy_d_qval[i]; 1769 1770 /* compute h*2^n-1 mod 'D' quickly */ 1771 val_mod = (h*pmod(2,n%(d-1),d)-1) % d; 1772 1773 /* 1774 * if 'D' mod 4 == 1, then 1775 * (h*2^n-1) mod 'D' can not be a quadratic residue of 'D' 1776 * else 1777 * (h*2^n-1) mod 'D' must be a quadratic residue of 'D' 1778 */ 1779 if (d%4 == 1) { 1780 /* D mod 4 == 1, so check for J(D, h*2^n-1) == -1 */ 1781 if (jacobi(val_mod, d) == -1) { 1782 /* it worked, return the related v(1) value */ 1783 return legacy_v1_qval[i]; 1784 } 1785 } else { 1786 /* D mod 4 == -1, so check for J(D, h*2^n-1) == 1 */ 1787 if (jacobi(val_mod, d) == 1) { 1788 /* it worked, return the related v(1) value */ 1789 return legacy_v1_qval[i]; 1790 } 1791 } 1792 } 1793 1794 /* 1795 * This is an example of a more complex proof construction. 1796 * The code above will not be able to find the v(1) for: 1797 * 1798 * 81*2^81-1 1799 * 1800 * We will check with: 1801 * 1802 * v(1)=81 D=6557 a=79 b=1 r=316 1803 * 1804 * Now, D==79*83 and r=79*2^2. If we show that: 1805 * 1806 * J(h*2^n-1 mod 79, 79) == -1 1807 * J(h*2^n-1 mod 83, 83) == 1 1808 * 1809 * then we will satisfy [condition 1]. Observe: 1810 * 1811 * 79 mod 4 == -1 implies (-1)^((h*2^n-2)*(79-1)/4) == -1 1812 * 83 mod 4 == -1 implies (-1)^((h*2^n-2)*(83-1)/4) == -1 1813 * 1814 * J(D, h*2^n-1) == J(83, h*2^n-1) * J(79, h*2^n-1) 1815 * == J(h*2^n-1, 83) * (-1)^((h*2^n-2)*(83-1)/4) * 1816 * J(h*2^n-1, 79) * (-1)^((h*2^n-2)*(79-1)/4) 1817 * == J(h*2^n-1 mod 83, 83) * -1 * 1818 * J(h*2^n-1 mod 79, 79) * -1 1819 * == 1 * -1 * 1820 * -1 * -1 1821 * == -1 1822 * 1823 * We will also satisfy [condition 2]. Observe: 1824 * 1825 * (a^2 - b^2*D)/r == (79^2 - 1^1*6557)/316 1826 * == -1 1827 * 1828 * L(r, h*2^n-1) == L(Q*(2^j)*(3^k)*(z^2), h*2^n-1) 1829 * == L(79, h*2^n-1) * L(2^2, h*2^n-1) 1830 * == L(79, h*2^n-1) * 1 1831 * == L(h*2^n-1, 79) * (-1)^((h*2^n-2)*(79-1)/4) 1832 * == L(h*2^n-1, 79) * -1 1833 * == L(h*2^n-1 mod 79, 79) * -1 1834 * == J(h*2^n-1 mod 79, 79) * -1 1835 * == -1 * -1 1836 * == 1 1837 */ 1838 if (jacobi( ((h*pmod(2,n%(79-1),79)-1)%79), 79 ) == -1 && 1839 jacobi( ((h*pmod(2,n%(83-1),83)-1)%83), 83 ) == 1) { 1840 /* return the associated v(1)=81 */ 1841 return 81; 1842 } 1843 1844 /* no quick and dirty v(1), so return -1 */ 1845 return -1; 1846} 1847