| blob | 24ba911b3dfa30bde3fb0047cd3ce284f1eb8e17 |
1 /*
2 * lucas - perform a Lucas primality test on h*2^n-1
3 *
4 * Copyright (C) 1999 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 * @(#) $Revision: 30.1 $
21 * @(#) $Id: lucas.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
22 * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/lucas.cal,v $
23 *
24 * Under source code control: 1990/05/03 16:49:51
25 * File existed as early as: 1990
26 *
27 * chongo <was here> /\oo/\ http://www.isthe.com/chongo/
28 * Share and enjoy! :-) http://www.isthe.com/chongo/tech/comp/calc/
29 */
31 /*
32 * NOTE: This is a standard calc resource file. For information on calc see:
33 *
34 * http://www.isthe.com/chongo/tech/comp/calc/index.html
35 *
36 * To obtain your own copy of calc, see:
37 *
38 * http://www.isthe.com/chongo/tech/comp/calc/calc-download.html
39 */
41 /*
42 * HISTORICAL NOTE:
43 *
44 * On 6 August 1989 at 00:53 PDT, the 'Amdahl 6', a team consisting of
45 * John Brown, Landon Curt Noll, Bodo Parady, Gene Smith, Joel Smith and
46 * Sergio Zarantonello proved the following 65087 digit number to be prime:
47 *
48 * 216193
49 * 391581 * 2 -1
50 *
51 * At the time of discovery, this number was the largest known prime.
52 * The primality was demonstrated by a program implementing the test
53 * found in these routines. An Amdahl 1200 takes 1987 seconds to test
54 * the primality of this number. A Cray 2 took several hours to
55 * confirm this prime. As of 31 Dec 1995, this prime was the 3rd
56 * largest known prime and the largest known non-Mersenne prime.
57 *
58 * The same team also discovered the following twin prime pair:
59 *
60 * 11235 11235
61 * 1706595 * 2 -1 1706595 * 2 +1
62 *
63 * At the time of discovery, this was the largest known twin prime pair.
64 *
65 * See:
66 *
67 * http://www.isthe.com/chongo/tech/math/prime/amdahl6.html
68 *
69 * for more information on the Amdahl 6 group.
70 *
71 * NOTE: Both largest known and largest known twin prime records have been
72 * broken. Rather than update this file each time, I'll just
73 * congratulate the finders and encourage others to try for
74 * larger finds. Records were made to be broken afterall!
75 */
77 /* ON GAINING A WORLD RECORD:
78 *
79 * The routines in calc were designed to be portable, and to work on
80 * numbers of 'sane' size. The Amdahl 6 team used a 'ultra-high speed
81 * multi-precision' package that a machine dependent collection of routines
82 * tuned for a long trace vector processor to work with very large numbers.
83 * The heart of the package was a multiplication and square routine that
84 * was based on the PFA Fast Fourier Transform and on Winograd's radix FFTs.
85 *
86 * Having a fast computer, and a good multi-precision package are
87 * critical, but one also needs to know where to look in order to have
88 * a good chance at a record. Knowing what to test is beyond the scope
89 * of this routine. However the following observations are noted:
90 *
91 * test numbers of the form h*2^n-1
92 * fix a value of n and vary the value h
93 * n mod 2^x == 0 for some value of x, say > 7 or more
94 * h*2^n-1 is not divisible by any small prime < 2^40
95 * 0 < h < 2^39
96 * h*2^n+1 is not divisible by any small prime < 2^40
97 *
98 * The Mersenne test for '2^n-1' is the fastest known primality test
99 * for a given large numbers. However, it is faster to search for
100 * primes of the form 'h*2^n-1'. When n is around 200000, one can find
101 * a prime of the form 'h*2^n-1' in about 1/2 the time.
102 *
103 * Critical to understanding why 'h*2^n-1' is to observe that primes of
104 * the form '2^n-1' seem to bunch around "islands". Such "islands"
105 * seem to be getting fewer and farther in-between, forcing the time
106 * for each test to grow longer and longer (worse then O(n^2 log n)).
107 * On the other hand, when one tests 'h*2^n-1', fixes 'n' and varies
108 * 'h', the time to test each number remains relatively constant.
109 *
110 * It is clearly a win to eliminate potential test candidates by
111 * rejecting numbers that that are divisible by 'small' primes. We
112 * (the "Amdahl 6") rejected all numbers that were divisible by primes
113 * less than '2^40'. We stopped looking for small factors at '2^40'
114 * when the rate of candidates being eliminated was slowed down to
115 * just a trickle.
116 *
117 * The 'n mod 128 == 0' restriction allows one to test for divisibility
118 * of small primes more quickly. To test of 'q' is a factor of 'k*2^n-1',
119 * one check to see if 'k*2^n mod q' == 1, which is the same a checking
120 * if 'h*(2^n mod q) mod q' == 1. One can compute '2^n mod q' by making
121 * use of the following:
122 *
123 * if
124 * y = 2^x mod q
125 * then
126 * 2^(2x) mod q == y^2 mod q 0 bit
127 * 2^(2x+1) mod q == 2*y^2 mod q 1 bit
128 *
129 * The choice of which expression depends on the binary pattern of 'n'.
130 * Since '1' bits require an extra step (multiply by 2), one should
131 * select value of 'n' that contain mostly '0' bits. The restriction
132 * of 'n mod 128 == 0' ensures that the bottom 7 bits of 'n' are 0.
133 *
134 * By limiting 'h' to '2^39' and eliminating all values divisible by
135 * small primes < twice the 'h' limit (2^40), one knows that all
136 * remaining candidates are relatively prime. Thus, when a candidate
137 * is proven to be composite (not prime) by the big test, one knows
138 * that the factors for that number (whatever they may be) will not
139 * be the factors of another candidate.
140 *
141 * Finally, one should eliminate all values of 'h*2^n-1' where
142 * 'h*2^n+1' is divisible by a small primes. The ideas behind this
143 * point is beyond the scope of this program.
144 */
147 global pprod256; /* product of "primes up to 256" / "primes up to 46" */
149 /*
150 * lucas - lucas primality test on h*2^n-1
151 *
152 * ABOUT THE TEST:
153 *
154 * This routine will perform a primality test on h*2^n-1 based on
155 * the mathematics of Lucas, Lehmer and Riesel. One should read
156 * the following article:
157 *
158 * Ref1:
159 * "Lucasian Criteria for the Primality of N=h*2^n-1", by Hans Riesel,
160 * Mathematics of Computation, Vol 23 #108, pp. 869-875, Oct 1969
161 *
162 * The following book is also useful:
163 *
164 * Ref2:
165 * "Prime numbers and Computer Methods for Factorization", by Hans Riesel,
166 * Birkhauser, 1985, pp 131-134, 278-285, 438-444
167 *
168 * A few useful Legendre identities may be found in:
169 *
170 * Ref3:
171 * "Introduction to Analytic Number Theory", by Tom A. Apostol,
172 * Springer-Verlag, 1984, p 188.
173 *
174 * This test is performed as follows: (see Ref1, Theorem 5)
175 *
176 * a) generate u(0) (see the function gen_u0() below)
177 *
178 * b) generate u(n-2) according to the rule:
179 *
180 * u(i+1) = u(i)^2-2 mod h*2^n-1
181 *
182 * c) h*2^n-1 is prime if and only if u(n-2) == 0 Q.E.D. :-)
183 *
184 * Now the following conditions must be true for the test to work:
185 *
186 * n >= 2
187 * h >= 1
188 * h < 2^n
189 * h mod 2 == 1
190 *
191 * A few misc notes:
192 *
193 * In order to reduce the number of tests, as attempt to eliminate
194 * any number that is divisible by a prime less than 257. Valid prime
195 * candidates less than 257 are declared prime as a special case.
196 *
197 * In real life, you would eliminate candidates by checking for
198 * divisibility by a prime much larger than 257 (perhaps as high
199 * as 2^39).
200 *
201 * The condition 'h mod 2 == 1' is not a problem. Say one is testing
202 * 'j*2^m-1', where j is even. If we note that:
203 *
204 * j mod 2^x == 0 for x>0 implies j*2^m-1 == ((j/2^x)*2^(m+x))-1,
205 *
206 * then we can let h=j/2^x and n=m+x and test 'h*2^n-1' which is the value.
207 * We need only consider odd values of h because we can rewrite our numbers
208 * do make this so.
209 *
210 * input:
211 * h the h as in h*2^n-1
212 * n the n as in h*2^n-1
213 *
214 * returns:
215 * 1 => h*2^n-1 is prime
216 * 0 => h*2^n-1 is not prime
217 * -1 => a test could not be formed, or h >= 2^n, h <= 0, n <= 0
218 */
219 define
220 lucas(h, n)
221 {
222 local testval; /* h*2^n-1 */
223 local shiftdown; /* the power of 2 that divides h */
224 local u; /* the u(i) sequence value */
225 local v1; /* the v(1) generator of u(0) */
226 local i; /* u sequence cycle number */
227 local oldh; /* pre-reduced h */
228 local oldn; /* pre-reduced n */
229 local bits; /* highbit of h*2^n-1 */
231 /*
232 * check arg types
233 */
234 if (!isint(h)) {
235 ldebug("lucas", "h is non-int");
236 quit "FATAL: bad args: h must be an integer";
237 }
238 if (!isint(n)) {
239 ldebug("lucas", "n is non-int");
240 quit "FATAL: bad args: n must be an integer";
241 }
243 /*
244 * reduce h if even
245 *
246 * we will force h to be odd by moving powers of two over to 2^n
247 */
248 oldh = h;
249 oldn = n;
250 shiftdown = fcnt(h,2); /* h % 2^shiftdown == 0, max shiftdown */
251 if (shiftdown > 0) {
252 h >>= shiftdown;
253 n += shiftdown;
254 }
256 /*
257 * enforce the 0 < h < 2^n rule
258 */
259 if (h <= 0 || n <= 0) {
260 print "ERROR: reduced args violate the rule: 0 < h < 2^n";
261 print " ERROR: h=":oldh, "n=":oldn, "reduced h=":h, "n=":n;
262 ldebug("lucas", "unknown: h <= 0 || n <= 0");
263 return -1;
264 }
265 if (highbit(h) >= n) {
266 print "ERROR: reduced args violate the rule: h < 2^n";
267 print " ERROR: h=":oldh, "n=":oldn, "reduced h=":h, "n=":n;
268 ldebug("lucas", "unknown: highbit(h) >= n");
269 return -1;
270 }
272 /*
273 * catch the degenerate case of h*2^n-1 == 1
274 */
275 if (h == 1 && n == 1) {
276 ldebug("lucas", "not prime: h == 1 && n == 1");
277 return 0; /* 1*2^1-1 == 1 is not prime */
278 }
280 /*
281 * catch the degenerate case of n==2
282 *
283 * n==2 and 0<h<2^n ==> 0<h<4
284 *
285 * Since h is now odd ==> h==1 or h==3
286 */
287 if (h == 1 && n == 2) {
288 ldebug("lucas", "prime: h == 1 && n == 2");
289 return 1; /* 1*2^2-1 == 3 is prime */
290 }
291 if (h == 3 && n == 2) {
292 ldebug("lucas", "prime: h == 3 && n == 2");
293 return 1; /* 3*2^2-1 == 11 is prime */
294 }
296 /*
297 * catch small primes < 257
298 *
299 * We check for only a few primes because the other primes < 257
300 * violate the checks above.
301 */
302 if (h == 1) {
303 if (n == 3 || n == 5 || n == 7) {
304 ldebug("lucas", "prime: 3, 7, 31, 127 are prime");
305 return 1; /* 3, 7, 31, 127 are prime */
306 }
307 }
308 if (h == 3) {
309 if (n == 2 || n == 3 || n == 4 || n == 6) {
310 ldebug("lucas", "prime: 11, 23, 47, 191 are prime");
311 return 1; /* 11, 23, 47, 191 are prime */
312 }
313 }
314 if (h == 5 && n == 4) {
315 ldebug("lucas", "prime: 79 is prime");
316 return 1; /* 79 is prime */
317 }
318 if (h == 7 && n == 5) {
319 ldebug("lucas", "prime: 223 is prime");
320 return 1; /* 223 is prime */
321 }
322 if (h == 15 && n == 4) {
323 ldebug("lucas", "prime: 239 is prime");
324 return 1; /* 239 is prime */
325 }
327 /*
328 * Avoid any numbers divisible by small primes
329 */
330 /*
331 * check for 3 <= prime factors < 29
332 * pfact(28)/2 = 111546435
333 */
334 testval = h*2^n - 1;
335 if (gcd(testval, 111546435) > 1) {
336 /* a small 3 <= prime < 29 divides h*2^n-1 */
337 ldebug("lucas","not-prime: 3<=prime<29 divides h*2^n-1");
338 return 0;
339 }
340 /*
341 * check for 29 <= prime factors < 47
342 * pfact(46)/pfact(28) = 5864229
343 */
344 if (gcd(testval, 58642669) > 1) {
345 /* a small 29 <= prime < 47 divides h*2^n-1 */
346 ldebug("lucas","not-prime: 29<=prime<47 divides h*2^n-1");
347 return 0;
348 }
349 /*
350 * check for prime 47 <= factors < 257, if h*2^n-1 is large
351 * 2^282 > pfact(256)/pfact(46) > 2^281
352 */
353 bits = highbit(testval);
354 if (bits >= 281) {
355 if (pprod256 <= 0) {
356 pprod256 = pfact(256)/pfact(46);
357 }
358 if (gcd(testval, pprod256) > 1) {
359 /* a small 47 <= prime < 257 divides h*2^n-1 */
360 ldebug("lucas",\
361 "not-prime: 47<=prime<257 divides h*2^n-1");
362 return 0;
363 }
364 }
366 /*
367 * try to compute u(0)
368 *
369 * We will use gen_v1() to give us a v(1) using the values
370 * of 'h' and 'n'. We will then use gen_u0() to convert
371 * the v(1) into u(0).
372 *
373 * If gen_v1() returns a negative value, then we failed to
374 * generate a test for h*2^n-1. This is because h mod 3 == 0
375 * is hard to do, and in rare cases, exceed the tables found
376 * in this program. We will generate an message and assume
377 * the number is not prime, even though if we had a larger
378 * table, we might have been able to show that it is prime.
379 */
380 v1 = gen_v1(h, n);
381 if (v1 < 0) {
382 /* failure to test number */
383 print "unable to compute v(1) for", h : "*2^" : n : "-1";
384 ldebug("lucas", "unknown: no v(1)");
385 return -1;
386 }
387 u = gen_u0(h, n, v1);
389 /*
390 * compute u(n-2)
391 */
392 for (i=3; i <= n; ++i) {
393 /* u = (u^2 - 2) % testval; */
394 u = hnrmod(u^2 - 2, h, n, -1);
395 }
397 /*
398 * return 1 if prime, 0 is not prime
399 */
400 if (u == 0) {
401 ldebug("lucas", "prime: end of test");
402 return 1;
403 } else {
404 ldebug("lucas", "not-prime: end of test");
405 return 0;
406 }
407 }
409 /*
410 * gen_u0 - determine the initial Lucas sequence for h*2^n-1
411 *
412 * According to Ref1, Theorem 5:
413 *
414 * u(0) = alpha^h + alpha^(-h)
415 *
416 * Now:
417 *
418 * v(x) = alpha^x + alpha^(-x) (Ref1, bottom of page 872)
419 *
420 * Therefore:
421 *
422 * u(0) = v(h)
423 *
424 * We calculate v(h) as follows: (Ref1, top of page 873)
425 *
426 * v(0) = alpha^0 + alpha^(-0) = 2
427 * v(1) = alpha^1 + alpha^(-1) = gen_v1(h,n)
428 * v(n+2) = v(1)*v(n+1) - v(n)
429 *
430 * This function does not concern itself with the value of 'alpha'.
431 * The gen_v1() function is used to compute v(1), and identity
432 * functions take it from there.
433 *
434 * It can be shown that the following are true:
435 *
436 * v(2*n) = v(n)^2 - 2
437 * v(2*n+1) = v(n+1)*v(n) - v(1)
438 *
439 * To prevent v(x) from growing too large, one may replace v(x) with
440 * `v(x) mod h*2^n-1' at any time.
441 *
442 * See the function gen_v1() for details on the value of v(1).
443 *
444 * input:
445 * h - h as in h*2^n-1 (h mod 2 != 0)
446 * n - n as in h*2^n-1
447 * v1 - gen_v1(h,n) (see function below)
448 *
449 * returns:
450 * u(0) - initial value for Lucas test on h*2^n-1
451 * -1 - failed to generate u(0)
452 */
453 define
454 gen_u0(h, n, v1)
455 {
456 local shiftdown; /* the power of 2 that divides h */
457 local r; /* low value: v(n) */
458 local s; /* high value: v(n+1) */
459 local hbits; /* highest bit set in h */
460 local i;
462 /*
463 * check arg types
464 */
465 if (!isint(h)) {
466 quit "bad args: h must be an integer";
467 }
468 if (!isint(n)) {
469 quit "bad args: n must be an integer";
470 }
471 if (!isint(v1)) {
472 quit "bad args: v1 must be an integer";
473 }
474 if (v1 <= 0) {
475 quit "bogus arg: v1 is <= 0";
476 }
478 /*
479 * enforce the h mod rules
480 */
481 if (h%2 == 0) {
482 quit "h must not be even";
483 }
485 /*
486 * enforce the h > 0 and n >= 2 rules
487 */
488 if (h <= 0 || n < 1) {
489 quit "reduced args violate the rule: 0 < h < 2^n";
490 }
491 hbits = highbit(h);
492 if (hbits >= n) {
493 quit "reduced args violate the rule: 0 < h < 2^n";
494 }
496 /*
497 * build up u2 based on the reversed bits of h
498 */
499 /* setup for bit loop */
500 r = v1;
501 s = (r^2 - 2);
503 /*
504 * deal with small h as a special case
505 *
506 * The h value is odd > 0, and it needs to be
507 * at least 2 bits long for the loop below to work.
508 */
509 if (h == 1) {
510 ldebug("gen_u0", "quick h == 1 case");
511 /* return r%(h*2^n-1); */
512 return hnrmod(r, h, n, -1);
513 }
515 /* cycle from second highest bit to second lowest bit of h */
516 for (i=hbits-1; i > 0; --i) {
518 /* bit(i) is 1 */
519 if (bit(h,i)) {
521 /* compute v(2n+1) = v(r+1)*v(r)-v1 */
522 /* r = (r*s - v1) % (h*2^n-1); */
523 r = hnrmod((r*s - v1), h, n, -1);
525 /* compute v(2n+2) = v(r+1)^2-2 */
526 /* s = (s^2 - 2) % (h*2^n-1); */
527 s = hnrmod((s^2 - 2), h, n, -1);
529 /* bit(i) is 0 */
530 } else {
532 /* compute v(2n+1) = v(r+1)*v(r)-v1 */
533 /* s = (r*s - v1) % (h*2^n-1); */
534 s = hnrmod((r*s - v1), h, n, -1);
536 /* compute v(2n) = v(r)^-2 */
537 /* r = (r^2 - 2) % (h*2^n-1); */
538 r = hnrmod((r^2 - 2), h, n, -1);
539 }
540 }
542 /* we know that h is odd, so the final bit(0) is 1 */
543 /* r = (r*s - v1) % (h*2^n-1); */
544 r = hnrmod((r*s - v1), h, n, -1);
546 /* compute the final u2 return value */
547 return r;
548 }
550 /*
551 * Trial tables used by gen_v1()
552 *
553 * When h mod 3 == 0, one needs particular values of D, a and b (see gen_v1
554 * documentation) in order to find a value of v(1).
555 *
556 * This table defines 'quickmax' possible tests to be taken in ascending
557 * order. The v1_qval[x] refers to a v(1) value from Ref1, Table 1. A
558 * related D value is found in d_qval[x]. All D values expect d_qval[1]
559 * are also taken from Ref1, Table 1. The case of D == 21 as listed in
560 * Ref1, Table 1 can be changed to D == 7 for the sake of the test because
561 * of {note 6}.
562 *
563 * It should be noted that the D values all satisfy the selection values
564 * as outlined in the gen_v1() function comments. That is:
565 *
566 * D == P*(2^f)*(3^g)
567 *
568 * where f == 0 and g == 0, P == D. So we simply need to check that
569 * one of the following two cases are true:
570 *
571 * P mod 4 == 1 and J(h*2^n-1 mod P, P) == -1
572 * P mod 4 == -1 and J(h*2^n-1 mod P, P) == 1
573 *
574 * In all cases, the value of r is:
575 *
576 * r == Q*(2^j)*(3^k)*(z^2)
577 *
578 * where Q == 1. No further processing is needed to compute v(1) when r
579 * is of this form.
580 */
581 quickmax = 8;
582 mat d_qval[quickmax];
583 mat v1_qval[quickmax];
584 d_qval[0] = 5; v1_qval[0] = 3; /* a=1 b=1 r=4 */
585 d_qval[1] = 7; v1_qval[1] = 5; /* a=3 b=1 r=12 D=21 */
586 d_qval[2] = 13; v1_qval[2] = 11; /* a=3 b=1 r=4 */
587 d_qval[3] = 11; v1_qval[3] = 20; /* a=3 b=1 r=2 */
588 d_qval[4] = 29; v1_qval[4] = 27; /* a=5 b=1 r=4 */
589 d_qval[5] = 53; v1_qval[5] = 51; /* a=53 b=1 r=4 */
590 d_qval[6] = 17; v1_qval[6] = 66; /* a=17 b=1 r=1 */
591 d_qval[7] = 19; v1_qval[7] = 74; /* a=38 b=1 r=2 */
593 /*
594 * gen_v1 - compute the v(1) for a given h*2^n-1 if we can
595 *
596 * This function assumes:
597 *
598 * n > 2 (n==2 has already been eliminated)
599 * h mod 2 == 1
600 * h < 2^n
601 * h*2^n-1 mod 3 != 0 (h*2^n-1 has no small factors, such as 3)
602 *
603 * The generation of v(1) depends on the value of h. There are two cases
604 * to consider, h mod 3 != 0, and h mod 3 == 0.
605 *
606 ***
607 *
608 * Case 1: (h mod 3 != 0)
609 *
610 * This case is easy and always finds v(1).
611 *
612 * In Ref1, page 869, one finds that if: (or see Ref2, page 131-132)
613 *
614 * h mod 6 == +/-1
615 * h*2^n-1 mod 3 != 0
616 *
617 * which translates, gives the functions assumptions, into the condition:
618 *
619 * h mod 3 != 0
620 *
621 * If this case condition is true, then:
622 *
623 * u(0) = (2+sqrt(3))^h + (2-sqrt(3))^h (see Ref1, page 869)
624 * = (2+sqrt(3))^h + (2+sqrt(3))^(-h)
625 *
626 * and since Ref1, Theorem 5 states:
627 *
628 * u(0) = alpha^h + alpha^(-h)
629 * r = abs(2^2 - 1^2*3) = 1
630 *
631 * and the bottom of Ref1, page 872 states:
632 *
633 * v(x) = alpha^x + alpha^(-x)
634 *
635 * If we let:
636 *
637 * alpha = (2+sqrt(3))
638 *
639 * then
640 *
641 * u(0) = v(h)
642 *
643 * so we simply return
644 *
645 * v(1) = alpha^1 + alpha^(-1)
646 * = (2+sqrt(3)) + (2-sqrt(3))
647 * = 4
648 *
649 ***
650 *
651 * Case 2: (h mod 3 == 0)
652 *
653 * This case is not so easy and finds v(1) in most all cases. In this
654 * version of this program, we will simply return -1 (failure) if we
655 * hit one of the cases that fall thru the cracks. This does not happen
656 * often, so this is not too bad.
657 *
658 * Ref1, Theorem 5 contains the following definitions:
659 *
660 * r = abs(a^2 - b^2*D)
661 * alpha = (a + b*sqrt(D))^2/r
662 *
663 * where D is 'square free', and 'alpha = epsilon^s' (for some s>0) are units
664 * in the quadratic field K(sqrt(D)).
665 *
666 * One can find possible values for a, b and D in Ref1, Table 1 (page 872).
667 * (see the file lucas_tbl.cal)
668 *
669 * Now Ref1, Theorem 5 states that if:
670 *
671 * L(D, h*2^n-1) = -1 [condition 1]
672 * L(r, h*2^n-1) * (a^2 - b^2*D)/r = -1 [condition 2]
673 *
674 * where L(x,y) is the Legendre symbol (see below), then:
675 *
676 * u(0) = alpha^h + alpha^(-h)
677 *
678 * The bottom of Ref1, page 872 states:
679 *
680 * v(x) = alpha^x + alpha^(-x)
681 *
682 * thus since:
683 *
684 * u(0) = v(h)
685 *
686 * so we want to return:
687 *
688 * v(1) = alpha^1 + alpha^(-1)
689 *
690 * Therefore we need to take a given (D,a,b), determine if the two conditions
691 * are true, and return the related v(1).
692 *
693 * Before we address the two conditions, we need some background information
694 * on two symbols, Legendre and Jacobi. In Ref 2, pp 278, 284-285, we find
695 * the following definitions of J(a,p) and L(a,n):
696 *
697 * The Legendre symbol L(a,p) takes the value:
698 *
699 * L(a,p) == 1 => a is a quadratic residue of p
700 * L(a,p) == -1 => a is NOT a quadratic residue of p
701 *
702 * when
703 *
704 * p is prime
705 * p mod 2 == 1
706 * gcd(a,p) == 1
707 *
708 * The value x is a quadratic residue of y if there exists some integer z
709 * such that:
710 *
711 * z^2 mod y == x
712 *
713 * The Jacobi symbol J(x,y) takes the value:
714 *
715 * J(x,y) == 1 => y is not prime, or x is a quadratic residue of y
716 * J(x,y) == -1 => x is NOT a quadratic residue of y
717 *
718 * when
719 *
720 * y mod 2 == 1
721 * gcd(x,y) == 1
722 *
723 * In the following comments on Legendre and Jacobi identities, we shall
724 * assume that the arguments to the symbolic are valid over the symbol
725 * definitions as stated above.
726 *
727 * In Ref2, pp 280-284, we find that:
728 *
729 * L(a,p)*L(b,p) == L(a*b,p) {A3.5}
730 * J(x,y)*J(z,y) == J(x*z,y) {A3.14}
731 * L(a,p) == L(p,a) * (-1)^((a-1)*(p-1)/4) {A3.8}
732 * J(x,y) == J(y,x) * (-1)^((x-1)*(y-1)/4) {A3.17}
733 *
734 * The equality L(a,p) == J(a,p) when: {note 0}
735 *
736 * p is prime
737 * p mod 2 == 1
738 * gcd(a,p) == 1
739 *
740 * It can be shown that (see Ref3):
741 *
742 * L(a,p) == L(a mod p, p) {note 1}
743 * L(z^2, p) == 1 {note 2}
744 *
745 * From Ref2, table 32:
746 *
747 * p mod 8 == +/-1 implies L(2,p) == 1 {note 3}
748 * p mod 12 == +/-1 implies L(3,p) == 1 {note 4}
749 *
750 * Since h*2^n-1 mod 8 == -1, for n>2, note 3 implies:
751 *
752 * L(2, h*2^n-1) == 1 (n>2) {note 5}
753 *
754 * Since h=3*A, h*2^n-1 mod 12 == -1, for A>0, note 4 implies:
755 *
756 * L(3, h*2^n-1) == 1 {note 6}
757 *
758 * By use of {A3.5}, {note 2}, {note 5} and {note 6}, one can show:
759 *
760 * L((2^g)*(3^l)*(z^2), h*2^n-1) == 1 (g>=0,l>=0,z>0,n>2) {note 7}
761 *
762 * Returning to the testing of conditions, take condition 1:
763 *
764 * L(D, h*2^n-1) == -1 [condition 1]
765 *
766 * In order for J(D, h*2^n-1) to be defined, we must ensure that D
767 * is not a factor of h*2^n-1. This is done by pre-screening h*2^n-1 to
768 * not have small factors and selecting D less than that factor check limit.
769 *
770 * By use of {note 7}, we can show that when we choose D to be:
771 *
772 * D is square free
773 * D = P*(2^f)*(3^g) (P is prime>2)
774 *
775 * The square free condition implies f = 0 or 1, g = 0 or 1. If f and g
776 * are both 1, P must be a prime > 3.
777 *
778 * So given such a D value:
779 *
780 * L(D, h*2^n-1) == L(P*(2^g)*(3^l), h*2^n-1)
781 * == L(P, h*2^n-1) * L((2^g)*(3^l), h*2^n-1) {A3.5}
782 * == L(P, h*2^n-1) * 1 {note 7}
783 * == L(h*2^n-1, P)*(-1)^((h*2^n-2)*(P-1)/4) {A3.8}
784 * == L(h*2^n-1 mod P, P)*(-1)^((h*2^n-2)*(P-1)/4) {note 1}
785 * == J(h*2^n-1 mod P, P)*(-1)^((h*2^n-2)*(P-1)/4) {note 0}
786 *
787 * When does J(h*2^n-1 mod P, P)*(-1)^((h*2^n-2)*(P-1)/4) take the value of -1,
788 * thus satisfy [condition 1]? The answer depends on P. Now P is a prime>2,
789 * thus P mod 4 == 1 or -1.
790 *
791 * Take P mod 4 == 1:
792 *
793 * P mod 4 == 1 implies (-1)^((h*2^n-2)*(P-1)/4) == 1
794 *
795 * Thus:
796 *
797 * L(D, h*2^n-1) == L(h*2^n-1 mod P, P) * (-1)^((h*2^n-2)*(P-1)/4)
798 * == L(h*2^n-1 mod P, P)
799 * == J(h*2^n-1 mod P, P)
800 *
801 * Take P mod 4 == -1:
802 *
803 * P mod 4 == -1 implies (-1)^((h*2^n-2)*(P-1)/4) == -1
804 *
805 * Thus:
806 *
807 * L(D, h*2^n-1) == L(h*2^n-1 mod P, P) * (-1)^((h*2^n-2)*(P-1)/4)
808 * == L(h*2^n-1 mod P, P) * -1
809 * == -J(h*2^n-1 mod P, P)
810 *
811 * Therefore [condition 1] is met if, and only if, one of the following
812 * to cases are true:
813 *
814 * P mod 4 == 1 and J(h*2^n-1 mod P, P) == -1
815 * P mod 4 == -1 and J(h*2^n-1 mod P, P) == 1
816 *
817 * Now consider [condition 2]:
818 *
819 * L(r, h*2^n-1) * (a^2 - b^2*D)/r == -1 [condition 2]
820 *
821 * We select only a, b, r and D values where:
822 *
823 * (a^2 - b^2*D)/r == -1
824 *
825 * Therefore in order for [condition 2] to be met, we must show that:
826 *
827 * L(r, h*2^n-1) == 1
828 *
829 * If we select r to be of the form:
830 *
831 * r == Q*(2^j)*(3^k)*(z^2) (Q == 1, j>=0, k>=0, z>0)
832 *
833 * then by use of {note 7}:
834 *
835 * L(r, h*2^n-1) == L(Q*(2^j)*(3^k)*(z^2), h*2^n-1)
836 * == L((2^j)*(3^k)*(z^2), h*2^n-1)
837 * == 1 {note 2}
838 *
839 * and thus, [condition 2] is met.
840 *
841 * If we select r to be of the form:
842 *
843 * r == Q*(2^j)*(3^k)*(z^2) (Q is prime>2, j>=0, k>=0, z>0)
844 *
845 * then by use of {note 7}:
846 *
847 * L(r, h*2^n-1) == L(Q*(2^j)*(3^k)*(z^2), h*2^n-1)
848 * == L(Q, h*2^n-1) * L((2^j)*(3^k)*(z^2), h*2^n-1) {A3.5}
849 * == L(Q, h*2^n-1) * 1 {note 2}
850 * == L(h*2^n-1, Q) * (-1)^((h*2^n-2)*(Q-1)/4) {A3.8}
851 * == L(h*2^n-1 mod Q, Q)*(-1)^((h*2^n-2)*(Q-1)/4) {note 1}
852 * == J(h*2^n-1 mod Q, Q)*(-1)^((h*2^n-2)*(Q-1)/4) {note 0}
853 *
854 * When does J(h*2^n-1 mod Q, Q)*(-1)^((h*2^n-2)*(Q-1)/4) take the value of 1,
855 * thus satisfy [condition 2]? The answer depends on Q. Now Q is a prime>2,
856 * thus Q mod 4 == 1 or -1.
857 *
858 * Take Q mod 4 == 1:
859 *
860 * Q mod 4 == 1 implies (-1)^((h*2^n-2)*(Q-1)/4) == 1
861 *
862 * Thus:
863 *
864 * L(D, h*2^n-1) == L(h*2^n-1 mod Q, Q) * (-1)^((h*2^n-2)*(Q-1)/4)
865 * == L(h*2^n-1 mod Q, Q)
866 * == J(h*2^n-1 mod Q, Q)
867 *
868 * Take Q mod 4 == -1:
869 *
870 * Q mod 4 == -1 implies (-1)^((h*2^n-2)*(Q-1)/4) == -1
871 *
872 * Thus:
873 *
874 * L(D, h*2^n-1) == L(h*2^n-1 mod Q, Q) * (-1)^((h*2^n-2)*(Q-1)/4)
875 * == L(h*2^n-1 mod Q, Q) * -1
876 * == -J(h*2^n-1 mod Q, Q)
877 *
878 * Therefore [condition 2] is met by selecting D = Q*(2^j)*(3^k)*(z^2),
879 * where Q is prime>2, j>=0, k>=0, z>0; if and only if one of the following
880 * to cases are true:
881 *
882 * Q mod 4 == 1 and J(h*2^n-1 mod Q, Q) == 1
883 * Q mod 4 == -1 and J(h*2^n-1 mod Q, Q) == -1
884 *
885 ***
886 *
887 * In conclusion, we can compute v(1) by attempting to do the following:
888 *
889 * h mod 3 != 0
890 *
891 * we return:
892 *
893 * v(1) == 4
894 *
895 * h mod 3 == 0
896 *
897 * define:
898 *
899 * r == abs(a^2 - b^2*D)
900 * alpha == (a + b*sqrt(D))^2/r
901 *
902 * we return:
903 *
904 * v(1) = alpha^1 + alpha^(-1)
905 *
906 * if and only if we can find a given a, b, D that obey all the
907 * following selection rules:
908 *
909 * D is square free
910 *
911 * D == P*(2^f)*(3^g) (P is prime>2, f,g == 0 or 1)
912 *
913 * (a^2 - b^2*D)/r == -1
914 *
915 * r == Q*(2^j)*(3^k)*(z^2) (Q==1 or Q is prime>2, j>=0, k>=0, z>0)
916 *
917 * one of the following is true:
918 * P mod 4 == 1 and J(h*2^n-1 mod P, P) == -1
919 * P mod 4 == -1 and J(h*2^n-1 mod P, P) == 1
920 *
921 * if Q is prime, then one of the following is true:
922 * Q mod 4 == 1 and J(h*2^n-1 mod Q, Q) == 1
923 * Q mod 4 == -1 and J(h*2^n-1 mod Q, Q) == -1
924 *
925 * If we cannot find a v(1) quickly enough, then we will give up
926 * testing h*2^n-1. This does not happen too often, so this hack
927 * is not too bad.
928 *
929 ***
930 *
931 * input:
932 * h h as in h*2^n-1
933 * n n as in h*2^n-1
934 *
935 * output:
936 * returns v(1), or -1 is there is no quick way
937 */
938 define
939 gen_v1(h, n)
940 {
941 local d; /* the 'D' value to try */
942 local val_mod; /* h*2^n-1 mod 'D' */
943 local i;
945 /*
946 * check for case 1
947 */
948 if (h % 3 != 0) {
949 /* v(1) is easy to compute */
950 return 4;
951 }
953 /*
954 * We will try all 'D' values until we find a proper v(1)
955 * or run out of 'D' values.
956 */
957 for (i=0; i < quickmax; ++i) {
959 /* grab our 'D' value */
960 d = d_qval[i];
962 /* compute h*2^n-1 mod 'D' quickly */
963 val_mod = (h*pmod(2,n%(d-1),d)-1) % d;
965 /*
966 * if 'D' mod 4 == 1, then
967 * (h*2^n-1) mod 'D' can not be a quadratic residue of 'D'
968 * else
969 * (h*2^n-1) mod 'D' must be a quadratic residue of 'D'
970 */
971 if (d%4 == 1) {
972 /* D mod 4 == 1, so check for J(D, h*2^n-1) == -1 */
973 if (jacobi(val_mod, d) == -1) {
974 /* it worked, return the related v(1) value */
975 return v1_qval[i];
976 }
977 } else {
978 /* D mod 4 == -1, so check for J(D, h*2^n-1) == 1 */
979 if (jacobi(val_mod, d) == 1) {
980 /* it worked, return the related v(1) value */
981 return v1_qval[i];
982 }
983 }
984 }
986 /*
987 * This is an example of a more complex proof construction.
988 * The code above will not be able to find the v(1) for:
989 *
990 * 81*2^81-1
991 *
992 * We will check with:
993 *
994 * v(1)=81 D=6557 a=79 b=1 r=316
995 *
996 * Now, D==79*83 and r=79*2^2. If we show that:
997 *
998 * J(h*2^n-1 mod 79, 79) == -1
999 * J(h*2^n-1 mod 83, 83) == 1
1000 *
1001 * then we will satisfy [condition 1]. Observe:
1002 *
1003 * 79 mod 4 == -1 implies (-1)^((h*2^n-2)*(79-1)/4) == -1
1004 * 83 mod 4 == -1 implies (-1)^((h*2^n-2)*(83-1)/4) == -1
1005 *
1006 * J(D, h*2^n-1) == J(83, h*2^n-1) * J(79, h*2^n-1)
1007 * == J(h*2^n-1, 83) * (-1)^((h*2^n-2)*(83-1)/4) *
1008 * J(h*2^n-1, 79) * (-1)^((h*2^n-2)*(79-1)/4)
1009 * == J(h*2^n-1 mod 83, 83) * -1 *
1010 * J(h*2^n-1 mod 79, 79) * -1
1011 * == 1 * -1 *
1012 * -1 * -1
1013 * == -1
1014 *
1015 * We will also satisfy [condition 2]. Observe:
1016 *
1017 * (a^2 - b^2*D)/r == (79^2 - 1^1*6557)/316
1018 * == -1
1019 *
1020 * L(r, h*2^n-1) == L(Q*(2^j)*(3^k)*(z^2), h*2^n-1)
1021 * == L(79, h*2^n-1) * L(2^2, h*2^n-1)
1022 * == L(79, h*2^n-1) * 1
1023 * == L(h*2^n-1, 79) * (-1)^((h*2^n-2)*(79-1)/4)
1024 * == L(h*2^n-1, 79) * -1
1025 * == L(h*2^n-1 mod 79, 79) * -1
1026 * == J(h*2^n-1 mod 79, 79) * -1
1027 * == -1 * -1
1028 * == 1
1029 */
1030 if (jacobi( ((h*pmod(2,n%(79-1),79)-1)%79), 79 ) == -1 &&
1031 jacobi( ((h*pmod(2,n%(83-1),83)-1)%83), 83 ) == 1) {
1032 /* return the associated v(1)=81 */
1033 return 81;
1034 }
1036 /* no quick and dirty v(1), so return -1 */
1037 return -1;
1038 }
1040 /*
1041 * ldebug - print a debug statement
1042 *
1043 * input:
1044 * funct name of calling function
1045 * str string to print
1046 */
1047 define
1048 ldebug(funct, str)
1049 {
1050 if (config("resource_debug") & 8) {
1051 print "DEBUG:", funct:":", str;
1052 }
1053 return;
1054 }