Annotation of OpenXM_contrib/gmp/tests/mpz/t-jac.c, Revision 1.1.1.1
1.1 ohara 1: /* Exercise mpz_*_kronecker_*() and mpz_jacobi() functions. */
2:
3: /*
4: Copyright 1999, 2000, 2001, 2002 Free Software Foundation, Inc.
5:
6: This file is part of the GNU MP Library.
7:
8: The GNU MP Library is free software; you can redistribute it and/or modify
9: it under the terms of the GNU Lesser General Public License as published by
10: the Free Software Foundation; either version 2.1 of the License, or (at your
11: option) any later version.
12:
13: The GNU MP Library is distributed in the hope that it will be useful, but
14: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16: License for more details.
17:
18: You should have received a copy of the GNU Lesser General Public License
19: along with the GNU MP Library; see the file COPYING.LIB. If not, write to
20: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
21: MA 02111-1307, USA.
22: */
23:
24:
25: /* With no arguments the various Kronecker/Jacobi symbol routines are
26: checked against some test data and a lot of derived data.
27:
28: To check the test data against PARI-GP, run
29:
30: t-jac -p | gp -q
31:
32: It takes a while because the output from "t-jac -p" is big.
33:
34:
35: Enhancements:
36:
37: More big test cases than those given by check_squares_zi would be good. */
38:
39:
40: #include <stdio.h>
41: #include <stdlib.h>
42: #include <string.h>
43:
44: #include "gmp.h"
45: #include "gmp-impl.h"
46: #include "tests.h"
47:
48:
49: #ifdef _LONG_LONG_LIMB
50: #define LL(l,ll) ll
51: #else
52: #define LL(l,ll) l
53: #endif
54:
55:
56: int option_pari = 0;
57:
58:
59: unsigned long
60: mpz_mod4 (mpz_srcptr z)
61: {
62: mpz_t m;
63: unsigned long ret;
64:
65: mpz_init (m);
66: mpz_fdiv_r_2exp (m, z, 2);
67: ret = mpz_get_ui (m);
68: mpz_clear (m);
69: return ret;
70: }
71:
72: int
73: mpz_fits_ulimb_p (mpz_srcptr z)
74: {
75: return (SIZ(z) == 1 || SIZ(z) == 0);
76: }
77:
78: mp_limb_t
79: mpz_get_ulimb (mpz_srcptr z)
80: {
81: if (SIZ(z) == 0)
82: return 0;
83: else
84: return PTR(z)[0];
85: }
86:
87:
88: void
89: try_base (mp_limb_t a, mp_limb_t b, int answer)
90: {
91: int got;
92:
93: if ((b & 1) == 0 || b == 1 || a > b)
94: return;
95:
96: got = mpn_jacobi_base (a, b, 0);
97: if (got != answer)
98: {
99: printf (LL("mpn_jacobi_base (%lu, %lu) is %d should be %d\n",
100: "mpn_jacobi_base (%llu, %llu) is %d should be %d\n"),
101: a, b, got, answer);
102: abort ();
103: }
104: }
105:
106:
107: void
108: try_zi_ui (mpz_srcptr a, unsigned long b, int answer)
109: {
110: int got;
111:
112: got = mpz_kronecker_ui (a, b);
113: if (got != answer)
114: {
115: printf ("mpz_kronecker_ui (");
116: mpz_out_str (stdout, 10, a);
117: printf (", %lu) is %d should be %d\n", b, got, answer);
118: abort ();
119: }
120: }
121:
122:
123: void
124: try_zi_si (mpz_srcptr a, long b, int answer)
125: {
126: int got;
127:
128: got = mpz_kronecker_si (a, b);
129: if (got != answer)
130: {
131: printf ("mpz_kronecker_si (");
132: mpz_out_str (stdout, 10, a);
133: printf (", %ld) is %d should be %d\n", b, got, answer);
134: abort ();
135: }
136: }
137:
138:
139: void
140: try_ui_zi (unsigned long a, mpz_srcptr b, int answer)
141: {
142: int got;
143:
144: got = mpz_ui_kronecker (a, b);
145: if (got != answer)
146: {
147: printf ("mpz_ui_kronecker (%lu, ", a);
148: mpz_out_str (stdout, 10, b);
149: printf (") is %d should be %d\n", got, answer);
150: abort ();
151: }
152: }
153:
154:
155: void
156: try_si_zi (long a, mpz_srcptr b, int answer)
157: {
158: int got;
159:
160: got = mpz_si_kronecker (a, b);
161: if (got != answer)
162: {
163: printf ("mpz_si_kronecker (%d, ", a);
164: mpz_out_str (stdout, 10, b);
165: printf (") is %d should be %d\n", got, answer);
166: abort ();
167: }
168: }
169:
170:
171: /* Don't bother checking mpz_jacobi, since it only differs for b even, and
172: we don't have an actual expected answer for it. tests/devel/try.c does
173: some checks though. */
174: void
175: try_zi_zi (mpz_srcptr a, mpz_srcptr b, int answer)
176: {
177: int got;
178:
179: got = mpz_kronecker (a, b);
180: if (got != answer)
181: {
182: printf ("mpz_kronecker (");
183: mpz_out_str (stdout, 10, a);
184: printf (", ");
185: mpz_out_str (stdout, 10, b);
186: printf (") is %d should be %d\n", got, answer);
187: abort ();
188: }
189: }
190:
191:
192: void
193: try_pari (mpz_srcptr a, mpz_srcptr b, int answer)
194: {
195: printf ("try(");
196: mpz_out_str (stdout, 10, a);
197: printf (",");
198: mpz_out_str (stdout, 10, b);
199: printf (",%d)\n", answer);
200: }
201:
202:
203: void
204: try_each (mpz_srcptr a, mpz_srcptr b, int answer)
205: {
206: if (option_pari)
207: {
208: try_pari (a, b, answer);
209: return;
210: }
211:
212: if (mpz_fits_ulimb_p (a) && mpz_fits_ulimb_p (b))
213: try_base (mpz_get_ulimb (a), mpz_get_ulimb (b), answer);
214:
215: if (mpz_fits_ulong_p (b))
216: try_zi_ui (a, mpz_get_ui (b), answer);
217:
218: if (mpz_fits_slong_p (b))
219: try_zi_si (a, mpz_get_si (b), answer);
220:
221: if (mpz_fits_ulong_p (a))
222: try_ui_zi (mpz_get_ui (a), b, answer);
223:
224: if (mpz_fits_sint_p (a))
225: try_si_zi (mpz_get_si (a), b, answer);
226:
227: try_zi_zi (a, b, answer);
228: }
229:
230:
231: /* Try (a/b) and (a/-b). */
232: void
233: try_pn (mpz_srcptr a, mpz_srcptr b_orig, int answer)
234: {
235: mpz_t b;
236:
237: mpz_init_set (b, b_orig);
238: try_each (a, b, answer);
239:
240: mpz_neg (b, b);
241: if (mpz_sgn (a) < 0)
242: answer = -answer;
243:
244: try_each (a, b, answer);
245:
246: mpz_clear (b);
247: }
248:
249:
250: /* Try (a+k*p/b) for various k, using the fact (a/b) is periodic in a with
251: period p. For b>0, p=b if b!=2mod4 or p=4*b if b==2mod4. */
252:
253: void
254: try_periodic_num (mpz_srcptr a_orig, mpz_srcptr b, int answer)
255: {
256: mpz_t a, a_period;
257: int i;
258:
259: if (mpz_sgn (b) <= 0)
260: return;
261:
262: mpz_init_set (a, a_orig);
263: mpz_init_set (a_period, b);
264: if (mpz_mod4 (b) == 2)
265: mpz_mul_ui (a_period, a_period, 4);
266:
267: /* don't bother with these tests if they're only going to produce
268: even/even */
269: if (mpz_even_p (a) && mpz_even_p (b) && mpz_even_p (a_period))
270: goto done;
271:
272: for (i = 0; i < 10; i++)
273: {
274: mpz_add (a, a, a_period);
275: try_pn (a, b, answer);
276: }
277:
278: mpz_set (a, a_orig);
279: for (i = 0; i < 10; i++)
280: {
281: mpz_sub (a, a, a_period);
282: try_pn (a, b, answer);
283: }
284:
285: done:
286: mpz_clear (a);
287: mpz_clear (a_period);
288: }
289:
290:
291: /* Try (a/b+k*p) for various k, using the fact (a/b) is periodic in b of
292: period p.
293:
294: period p
295: a==0,1mod4 a
296: a==2mod4 4*a
297: a==3mod4 and b odd 4*a
298: a==3mod4 and b even 8*a
299:
300: In Henri Cohen's book the period is given as 4*a for all a==2,3mod4, but
301: a counterexample would seem to be (3/2)=-1 which with (3/14)=+1 doesn't
302: have period 4*a (but rather 8*a with (3/26)=-1). Maybe the plain 4*a is
303: to be read as applying to a plain Jacobi symbol with b odd, rather than
304: the Kronecker extension to b even. */
305:
306: void
307: try_periodic_den (mpz_srcptr a, mpz_srcptr b_orig, int answer)
308: {
309: mpz_t b, b_period;
310: int i;
311:
312: if (mpz_sgn (a) == 0 || mpz_sgn (b_orig) == 0)
313: return;
314:
315: mpz_init_set (b, b_orig);
316:
317: mpz_init_set (b_period, a);
318: if (mpz_mod4 (a) == 3 && mpz_even_p (b))
319: mpz_mul_ui (b_period, b_period, 8);
320: else if (mpz_mod4 (a) >= 2)
321: mpz_mul_ui (b_period, b_period, 4);
322:
323: /* don't bother with these tests if they're only going to produce
324: even/even */
325: if (mpz_even_p (a) && mpz_even_p (b) && mpz_even_p (b_period))
326: goto done;
327:
328: for (i = 0; i < 10; i++)
329: {
330: mpz_add (b, b, b_period);
331: try_pn (a, b, answer);
332: }
333:
334: mpz_set (b, b_orig);
335: for (i = 0; i < 10; i++)
336: {
337: mpz_sub (b, b, b_period);
338: try_pn (a, b, answer);
339: }
340:
341: done:
342: mpz_clear (b);
343: mpz_clear (b_period);
344: }
345:
346:
347: /* Try (a/b*2^k) for various k. */
348: void
349: try_2den (mpz_srcptr a, mpz_srcptr b_orig, int answer)
350: {
351: mpz_t b;
352: int i;
353: int answer_two;
354:
355: /* don't bother when b==0 */
356: if (mpz_sgn (b_orig) == 0)
357: return;
358:
359: mpz_init_set (b, b_orig);
360:
361: /* answer_two = (a/2) */
362: switch (mpz_fdiv_ui (a, 8)) {
363: case 1:
364: case 7:
365: answer_two = 1;
366: break;
367: case 3:
368: case 5:
369: answer_two = -1;
370: break;
371: default:
372: /* 0, 2, 4, 6 */
373: answer_two = 0;
374: break;
375: }
376:
377: for (i = 0; i < 3 * BITS_PER_MP_LIMB; i++)
378: {
379: answer *= answer_two;
380: mpz_mul_2exp (b, b, 1);
381: try_pn (a, b, answer);
382: }
383:
384: mpz_clear (b);
385: }
386:
387:
388: /* Try (a*2^k/b) for various k. If it happens mpz_ui_kronecker() gets (2/b)
389: wrong it will show up as wrong answers demanded. */
390: void
391: try_2num (mpz_srcptr a_orig, mpz_srcptr b, int answer)
392: {
393: mpz_t a;
394: int i;
395: int answer_twos;
396:
397: /* don't bother when a==0 */
398: if (mpz_sgn (a_orig) == 0)
399: return;
400:
401: mpz_init_set (a, a_orig);
402: answer_twos = mpz_ui_kronecker (2, b);
403:
404: for (i = 0; i < 3 * BITS_PER_MP_LIMB; i++)
405: {
406: answer *= answer_twos;
407: mpz_mul_2exp (a, a, 1);
408: try_pn (a, b, answer);
409: }
410:
411: mpz_clear (a);
412: }
413:
414:
415: /* The try_2num() and try_2den() routines don't in turn call
416: try_periodic_num() and try_periodic_den() because it hugely increases the
417: number of tests performed, without obviously increasing coverage.
418:
419: Useful extra derived cases can be added here. */
420:
421: void
422: try_all (mpz_t a, mpz_t b, int answer)
423: {
424: try_pn (a, b, answer);
425: try_periodic_num (a, b, answer);
426: try_periodic_den (a, b, answer);
427: try_2num (a, b, answer);
428: try_2den (a, b, answer);
429: }
430:
431:
432: void
433: check_data (void)
434: {
435: static const struct {
436: const char *a;
437: const char *b;
438: int answer;
439:
440: } data[] = {
441:
442: /* Note that the various derived checks in try_all() reduce the cases
443: that need to be given here. */
444:
445: /* some zeros */
446: { "0", "0", 0 },
447: { "0", "2", 0 },
448: { "0", "6", 0 },
449: { "5", "0", 0 },
450: { "24", "60", 0 },
451:
452: /* (a/1) = 1, any a
453: In particular note (0/1)=1 so that (a/b)=(a mod b/b). */
454: { "0", "1", 1 },
455: { "1", "1", 1 },
456: { "2", "1", 1 },
457: { "3", "1", 1 },
458: { "4", "1", 1 },
459: { "5", "1", 1 },
460:
461: /* (0/b) = 0, b != 1 */
462: { "0", "3", 0 },
463: { "0", "5", 0 },
464: { "0", "7", 0 },
465: { "0", "9", 0 },
466: { "0", "11", 0 },
467: { "0", "13", 0 },
468: { "0", "15", 0 },
469:
470: /* (1/b) = 1 */
471: { "1", "1", 1 },
472: { "1", "3", 1 },
473: { "1", "5", 1 },
474: { "1", "7", 1 },
475: { "1", "9", 1 },
476: { "1", "11", 1 },
477:
478: /* (-1/b) = (-1)^((b-1)/2) which is -1 for b==3 mod 4 */
479: { "-1", "1", 1 },
480: { "-1", "3", -1 },
481: { "-1", "5", 1 },
482: { "-1", "7", -1 },
483: { "-1", "9", 1 },
484: { "-1", "11", -1 },
485: { "-1", "13", 1 },
486: { "-1", "15", -1 },
487: { "-1", "17", 1 },
488: { "-1", "19", -1 },
489:
490: /* (2/b) = (-1)^((b^2-1)/8) which is -1 for b==3,5 mod 8.
491: try_2num() will exercise multiple powers of 2 in the numerator. */
492: { "2", "1", 1 },
493: { "2", "3", -1 },
494: { "2", "5", -1 },
495: { "2", "7", 1 },
496: { "2", "9", 1 },
497: { "2", "11", -1 },
498: { "2", "13", -1 },
499: { "2", "15", 1 },
500: { "2", "17", 1 },
501:
502: /* (-2/b) = (-1)^((b^2-1)/8)*(-1)^((b-1)/2) which is -1 for b==5,7mod8.
503: try_2num() will exercise multiple powers of 2 in the numerator, which
504: will test that the shift in mpz_si_kronecker() uses unsigned not
505: signed. */
506: { "-2", "1", 1 },
507: { "-2", "3", 1 },
508: { "-2", "5", -1 },
509: { "-2", "7", -1 },
510: { "-2", "9", 1 },
511: { "-2", "11", 1 },
512: { "-2", "13", -1 },
513: { "-2", "15", -1 },
514: { "-2", "17", 1 },
515:
516: /* (a/2)=(2/a).
517: try_2den() will exercise multiple powers of 2 in the denominator. */
518: { "3", "2", -1 },
519: { "5", "2", -1 },
520: { "7", "2", 1 },
521: { "9", "2", 1 },
522: { "11", "2", -1 },
523:
524: /* Harriet Griffin, "Elementary Theory of Numbers", page 155, various
525: examples. */
526: { "2", "135", 1 },
527: { "135", "19", -1 },
528: { "2", "19", -1 },
529: { "19", "135", 1 },
530: { "173", "135", 1 },
531: { "38", "135", 1 },
532: { "135", "173", 1 },
533: { "173", "5", -1 },
534: { "3", "5", -1 },
535: { "5", "173", -1 },
536: { "173", "3", -1 },
537: { "2", "3", -1 },
538: { "3", "173", -1 },
539: { "253", "21", 1 },
540: { "1", "21", 1 },
541: { "21", "253", 1 },
542: { "21", "11", -1 },
543: { "-1", "11", -1 },
544:
545: /* Griffin page 147 */
546: { "-1", "17", 1 },
547: { "2", "17", 1 },
548: { "-2", "17", 1 },
549: { "-1", "89", 1 },
550: { "2", "89", 1 },
551:
552: /* Griffin page 148 */
553: { "89", "11", 1 },
554: { "1", "11", 1 },
555: { "89", "3", -1 },
556: { "2", "3", -1 },
557: { "3", "89", -1 },
558: { "11", "89", 1 },
559: { "33", "89", -1 },
560:
561: /* H. Davenport, "The Higher Arithmetic", page 65, the quadratic
562: residues and non-residues mod 19. */
563: { "1", "19", 1 },
564: { "4", "19", 1 },
565: { "5", "19", 1 },
566: { "6", "19", 1 },
567: { "7", "19", 1 },
568: { "9", "19", 1 },
569: { "11", "19", 1 },
570: { "16", "19", 1 },
571: { "17", "19", 1 },
572: { "2", "19", -1 },
573: { "3", "19", -1 },
574: { "8", "19", -1 },
575: { "10", "19", -1 },
576: { "12", "19", -1 },
577: { "13", "19", -1 },
578: { "14", "19", -1 },
579: { "15", "19", -1 },
580: { "18", "19", -1 },
581:
582: /* Residues and non-residues mod 13 */
583: { "0", "13", 0 },
584: { "1", "13", 1 },
585: { "2", "13", -1 },
586: { "3", "13", 1 },
587: { "4", "13", 1 },
588: { "5", "13", -1 },
589: { "6", "13", -1 },
590: { "7", "13", -1 },
591: { "8", "13", -1 },
592: { "9", "13", 1 },
593: { "10", "13", 1 },
594: { "11", "13", -1 },
595: { "12", "13", 1 },
596:
597: /* various */
598: { "5", "7", -1 },
599: { "15", "17", 1 },
600: { "67", "89", 1 },
601:
602: /* special values inducing a==b==1 at the end of jac_or_kron() */
603: { "0x10000000000000000000000000000000000000000000000001",
604: "0x10000000000000000000000000000000000000000000000003", 1 },
605: };
606:
607: int i, answer;
608: mpz_t a, b;
609:
610: mpz_init (a);
611: mpz_init (b);
612:
613: for (i = 0; i < numberof (data); i++)
614: {
615: mpz_set_str_or_abort (a, data[i].a, 0);
616: mpz_set_str_or_abort (b, data[i].b, 0);
617:
618: answer = data[i].answer;
619: try_all (a, b, data[i].answer);
620: }
621:
622: mpz_clear (a);
623: mpz_clear (b);
624: }
625:
626:
627: /* (a^2/b)=1 if gcd(a,b)=1, or (a^2/b)=0 if gcd(a,b)!=1.
628: This includes when a=0 or b=0. */
629: void
630: check_squares_zi (void)
631: {
632: gmp_randstate_ptr rands = RANDS;
633: mpz_t a, b, g;
634: int i, answer;
635: mp_size_t size_range, an, bn;
636: mpz_t bs;
637:
638: mpz_init (bs);
639: mpz_init (a);
640: mpz_init (b);
641: mpz_init (g);
642:
643: for (i = 0; i < 200; i++)
644: {
645: mpz_urandomb (bs, rands, 32);
646: size_range = mpz_get_ui (bs) % 10 + 2;
647:
648: mpz_urandomb (bs, rands, size_range);
649: an = mpz_get_ui (bs);
650: mpz_rrandomb (a, rands, an);
651:
652: mpz_urandomb (bs, rands, size_range);
653: bn = mpz_get_ui (bs);
654: mpz_rrandomb (b, rands, bn);
655:
656: mpz_gcd (g, a, b);
657: if (mpz_cmp_ui (g, 1L) == 0)
658: answer = 1;
659: else
660: answer = 0;
661:
662: mpz_mul (a, a, a);
663:
664: try_all (a, b, answer);
665: }
666:
667: mpz_clear (bs);
668: mpz_clear (a);
669: mpz_clear (b);
670: mpz_clear (g);
671: }
672:
673:
674: /* Check the handling of asize==0, make sure it isn't affected by the low
675: limb. */
676: void
677: check_a_zero (void)
678: {
679: mpz_t a, b;
680:
681: mpz_init_set_ui (a, 0);
682: mpz_init (b);
683:
684: mpz_set_ui (b, 1L);
685: PTR(a)[0] = 0;
686: try_all (a, b, 1); /* (0/1)=1 */
687: PTR(a)[0] = 1;
688: try_all (a, b, 1); /* (0/1)=1 */
689:
690: mpz_set_si (b, -1L);
691: PTR(a)[0] = 0;
692: try_all (a, b, 1); /* (0/-1)=1 */
693: PTR(a)[0] = 1;
694: try_all (a, b, 1); /* (0/-1)=1 */
695:
696: mpz_set_ui (b, 0);
697: PTR(a)[0] = 0;
698: try_all (a, b, 0); /* (0/0)=0 */
699: PTR(a)[0] = 1;
700: try_all (a, b, 0); /* (0/0)=0 */
701:
702: mpz_set_ui (b, 2);
703: PTR(a)[0] = 0;
704: try_all (a, b, 0); /* (0/2)=0 */
705: PTR(a)[0] = 1;
706: try_all (a, b, 0); /* (0/2)=0 */
707:
708: mpz_clear (a);
709: mpz_clear (b);
710: }
711:
712:
713: int
714: main (int argc, char *argv[])
715: {
716: tests_start ();
717:
718: if (argc >= 2 && strcmp (argv[1], "-p") == 0)
719: {
720: option_pari = 1;
721:
722: printf ("\
723: try(a,b,answer) =\n\
724: {\n\
725: if (kronecker(a,b) != answer,\n\
726: print(\"wrong at \", a, \",\", b,\n\
727: \" expected \", answer,\n\
728: \" pari says \", kronecker(a,b)))\n\
729: }\n");
730: }
731:
732: check_data ();
733: check_squares_zi ();
734: check_a_zero ();
735:
736: tests_end ();
737: exit (0);
738: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to t-jac.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / tests / mpz