Annotation of OpenXM_contrib/gmp/mpz/tests/t-jac.c, Revision 1.1.1.1
1.1 maekawa 1: /* Exercise mpz_*_kronecker_*() and mpz_jacobi() functions. */
2:
3: /*
4: Copyright (C) 1999, 2000 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 "gmp.h"
42: #include "gmp-impl.h"
43:
44:
45: int option_pari = 0;
46:
47:
48: unsigned long
49: mpz_mod4 (mpz_srcptr z)
50: {
51: mpz_t m;
52: unsigned long ret;
53:
54: mpz_init (m);
55: mpz_fdiv_r_2exp (m, z, 2);
56: ret = mpz_get_ui (m);
57: mpz_clear (m);
58: return ret;
59: }
60:
61: int
62: mpz_fits_ulimb_p (mpz_srcptr z)
63: {
64: return (SIZ(z) == 1 || SIZ(z) == 0);
65: }
66:
67: mp_limb_t
68: mpz_get_ulimb (mpz_srcptr z)
69: {
70: if (SIZ(z) == 0)
71: return 0;
72: else
73: return PTR(z)[0];
74: }
75:
76:
77: void
78: try_base (mp_limb_t a, mp_limb_t b, int answer)
79: {
80: int got;
81:
82: if ((b & 1) == 0 || b == 1 || a > b)
83: return;
84:
85: got = mpn_jacobi_base (a, b, 0);
86: if (got != answer)
87: {
88: printf ("mpn_jacobi_base (%lu, %lu) is %d should be %d\n",
89: a, b, got, answer);
90: exit (1);
91: }
92: }
93:
94:
95: void
96: try_zi_ui (mpz_srcptr a, unsigned long b, int answer)
97: {
98: int got;
99:
100: got = mpz_kronecker_ui (a, b);
101: if (got != answer)
102: {
103: printf ("mpz_kronecker_ui (");
104: mpz_out_str (stdout, 10, a);
105: printf (", %lu) is %d should be %d\n", b, got, answer);
106: exit (1);
107: }
108: }
109:
110:
111: void
112: try_zi_si (mpz_srcptr a, long b, int answer)
113: {
114: int got;
115:
116: got = mpz_kronecker_si (a, b);
117: if (got != answer)
118: {
119: printf ("mpz_kronecker_si (");
120: mpz_out_str (stdout, 10, a);
121: printf (", %ld) is %d should be %d\n", b, got, answer);
122: exit (1);
123: }
124: }
125:
126:
127: void
128: try_ui_zi (unsigned long a, mpz_srcptr b, int answer)
129: {
130: int got;
131:
132: got = mpz_ui_kronecker (a, b);
133: if (got != answer)
134: {
135: printf ("mpz_ui_kronecker (%lu, ", a);
136: mpz_out_str (stdout, 10, b);
137: printf (") is %d should be %d\n", got, answer);
138: exit (1);
139: }
140: }
141:
142:
143: void
144: try_si_zi (int a, mpz_srcptr b, int answer)
145: {
146: int got;
147:
148: got = mpz_si_kronecker (a, b);
149: if (got != answer)
150: {
151: printf ("mpz_si_kronecker (%d, ", a);
152: mpz_out_str (stdout, 10, b);
153: printf (") is %d should be %d\n", got, answer);
154: exit (1);
155: }
156: }
157:
158:
159: void
160: try_zi_zi (mpz_srcptr a, mpz_srcptr b, int answer)
161: {
162: int got;
163:
164: /* gmp 2.0.2 mpz_jacobi() doesn't handle negative or even b */
165: if (mpz_sgn (b) <= 0 || mpz_even_p (b))
166: return;
167:
168: got = mpz_jacobi (a, b);
169: if (got != answer)
170: {
171: printf ("mpz_jacobi (");
172: mpz_out_str (stdout, 10, a);
173: printf (", ");
174: mpz_out_str (stdout, 10, b);
175: printf (") is %d should be %d\n", got, answer);
176: exit (1);
177: }
178: }
179:
180:
181: void
182: try_pari (mpz_srcptr a, mpz_srcptr b, int answer)
183: {
184: printf ("try(");
185: mpz_out_str (stdout, 10, a);
186: printf (",");
187: mpz_out_str (stdout, 10, b);
188: printf (",%d)\n", answer);
189: }
190:
191:
192: void
193: try_each (mpz_srcptr a, mpz_srcptr b, int answer)
194: {
195: if (option_pari)
196: {
197: try_pari (a, b, answer);
198: return;
199: }
200:
201: if (mpz_fits_ulimb_p (a) && mpz_fits_ulimb_p (b))
202: try_base (mpz_get_ulimb (a), mpz_get_ulimb (b), answer);
203:
204: if (mpz_fits_ulong_p (b))
205: try_zi_ui (a, mpz_get_ui (b), answer);
206:
207: if (mpz_fits_slong_p (b))
208: try_zi_si (a, mpz_get_si (b), answer);
209:
210: if (mpz_fits_ulong_p (a))
211: try_ui_zi (mpz_get_ui (a), b, answer);
212:
213: if (mpz_fits_sint_p (a))
214: try_si_zi (mpz_get_si (a), b, answer);
215:
216: try_zi_zi (a, b, answer);
217: }
218:
219:
220: /* Try (a/b) and (a/-b). */
221: void
222: try_pn (mpz_srcptr a, mpz_srcptr b_orig, int answer)
223: {
224: mpz_t b;
225:
226: mpz_init_set (b, b_orig);
227: try_each (a, b, answer);
228:
229: mpz_neg (b, b);
230: if (mpz_sgn (a) < 0)
231: answer = -answer;
232:
233: try_each (a, b, answer);
234:
235: mpz_clear (b);
236: }
237:
238:
239: /* Try (a+k*p/b) for various k, using the fact (a/b) is periodic in a with
240: period p. For b>0, p=b if b!=2mod4 or p=4*b if b==2mod4. */
241:
242: void
243: try_periodic_num (mpz_srcptr a_orig, mpz_srcptr b, int answer)
244: {
245: mpz_t a, a_period;
246: int i;
247:
248: if (mpz_sgn (b) <= 0)
249: return;
250:
251: mpz_init_set (a, a_orig);
252: mpz_init_set (a_period, b);
253: if (mpz_mod4 (b) == 2)
254: mpz_mul_ui (a_period, a_period, 4);
255:
256: /* don't bother with these tests if they're only going to produce
257: even/even */
258: if (mpz_even_p (a) && mpz_even_p (b) && mpz_even_p (a_period))
259: goto done;
260:
261: for (i = 0; i < 10; i++)
262: {
263: mpz_add (a, a, a_period);
264: try_pn (a, b, answer);
265: }
266:
267: mpz_set (a, a_orig);
268: for (i = 0; i < 10; i++)
269: {
270: mpz_sub (a, a, a_period);
271: try_pn (a, b, answer);
272: }
273:
274: done:
275: mpz_clear (a);
276: mpz_clear (a_period);
277: }
278:
279:
280: /* Try (a/b+k*p) for various k, using the fact (a/b) is periodic in b of
281: period p.
282:
283: period p
284: a==0,1mod4 a
285: a==2mod4 4*a
286: a==3mod4 and b odd 4*a
287: a==3mod4 and b even 8*a
288:
289: In Henri Cohen's book the period is given as 4*a for all a==2,3mod4, but
290: a counterexample would seem to be (3/2)=-1 which with (3/14)=+1 doesn't
291: have period 4*a (but rather 8*a with (3/26)=-1). Maybe the plain 4*a is
292: to be read as applying to a plain Jacobi symbol with b odd, rather than
293: the Kronecker extension to b even. */
294:
295: void
296: try_periodic_den (mpz_srcptr a, mpz_srcptr b_orig, int answer)
297: {
298: mpz_t b, b_period;
299: int i;
300:
301: if (mpz_sgn (a) == 0 || mpz_sgn (b_orig) == 0)
302: return;
303:
304: mpz_init_set (b, b_orig);
305:
306: mpz_init_set (b_period, a);
307: if (mpz_mod4 (a) == 3 && mpz_even_p (b))
308: mpz_mul_ui (b_period, b_period, 8);
309: else if (mpz_mod4 (a) >= 2)
310: mpz_mul_ui (b_period, b_period, 4);
311:
312: /* don't bother with these tests if they're only going to produce
313: even/even */
314: if (mpz_even_p (a) && mpz_even_p (b) && mpz_even_p (b_period))
315: goto done;
316:
317: for (i = 0; i < 10; i++)
318: {
319: mpz_add (b, b, b_period);
320: try_pn (a, b, answer);
321: }
322:
323: mpz_set (b, b_orig);
324: for (i = 0; i < 10; i++)
325: {
326: mpz_sub (b, b, b_period);
327: try_pn (a, b, answer);
328: }
329:
330: done:
331: mpz_clear (b);
332: mpz_clear (b_period);
333: }
334:
335:
336: /* Try (a/b*2^k) for various k. If it happens mpz_ui_kronecker() gets (a/2)
337: wrong it will show up as wrong answers demanded. */
338: void
339: try_2den (mpz_srcptr a, mpz_srcptr b_orig, int answer)
340: {
341: mpz_t b;
342: int i;
343: int answer_two;
344:
345: /* don't bother when b==0 */
346: if (mpz_sgn (b_orig) == 0)
347: return;
348:
349: mpz_init_set (b, b_orig);
350: answer_two = mpz_kronecker_ui (a, 2);
351:
352: for (i = 0; i < 3 * BITS_PER_MP_LIMB; i++)
353: {
354: answer *= answer_two;
355: mpz_mul_2exp (b, b, 1);
356: try_pn (a, b, answer);
357: }
358:
359: mpz_clear (b);
360: }
361:
362:
363: /* Try (a*2^k/b) for various k. If it happens mpz_ui_kronecker() gets (2/b)
364: wrong it will show up as wrong answers demanded. */
365: void
366: try_2num (mpz_srcptr a_orig, mpz_srcptr b, int answer)
367: {
368: mpz_t a;
369: int i;
370: int answer_twos;
371:
372: /* don't bother when a==0 */
373: if (mpz_sgn (a_orig) == 0)
374: return;
375:
376: mpz_init_set (a, a_orig);
377: answer_twos = mpz_ui_kronecker (2, b);
378:
379: for (i = 0; i < 3 * BITS_PER_MP_LIMB; i++)
380: {
381: answer *= answer_twos;
382: mpz_mul_2exp (a, a, 1);
383: try_pn (a, b, answer);
384: }
385:
386: mpz_clear (a);
387: }
388:
389:
390: /* The try_2num() and try_2den() routines don't in turn call
391: try_periodic_num() and try_periodic_den() because it hugely increases the
392: number of tests performed, without obviously increasing coverage.
393:
394: Useful extra derived cases can be added here. */
395:
396: void
397: try_all (mpz_t a, mpz_t b, int answer)
398: {
399: try_pn (a, b, answer);
400: try_periodic_num (a, b, answer);
401: try_periodic_den (a, b, answer);
402: try_2num (a, b, answer);
403: try_2den (a, b, answer);
404: }
405:
406:
407: void
408: check_data (void)
409: {
410: static const struct {
411: const char *a;
412: const char *b;
413: int answer;
414:
415: } data[] = {
416:
417: /* Note that the various derived checks in try_all() reduce the cases
418: that need to be given here. */
419:
420: /* some zeros */
421: { "0", "0", 0 },
422: { "0", "2", 0 },
423: { "0", "6", 0 },
424: { "5", "0", 0 },
425: { "24", "60", 0 },
426:
427: /* (a/1) = 1, any a
428: In particular note (0/1)=1 so that (a/b)=(a mod b/b). */
429: { "0", "1", 1 },
430: { "1", "1", 1 },
431: { "2", "1", 1 },
432: { "3", "1", 1 },
433: { "4", "1", 1 },
434: { "5", "1", 1 },
435:
436: /* (0/b) = 0, b != 1 */
437: { "0", "3", 0 },
438: { "0", "5", 0 },
439: { "0", "7", 0 },
440: { "0", "9", 0 },
441: { "0", "11", 0 },
442: { "0", "13", 0 },
443: { "0", "15", 0 },
444:
445: /* (1/b) = 1 */
446: { "1", "1", 1 },
447: { "1", "3", 1 },
448: { "1", "5", 1 },
449: { "1", "7", 1 },
450: { "1", "9", 1 },
451: { "1", "11", 1 },
452:
453: /* (-1/b) = (-1)^((b-1)/2) which is -1 for b==3 mod 4 */
454: { "-1", "1", 1 },
455: { "-1", "3", -1 },
456: { "-1", "5", 1 },
457: { "-1", "7", -1 },
458: { "-1", "9", 1 },
459: { "-1", "11", -1 },
460: { "-1", "13", 1 },
461: { "-1", "15", -1 },
462: { "-1", "17", 1 },
463: { "-1", "19", -1 },
464:
465: /* (2/b) = (-1)^((b^2-1)/8) which is -1 for b==3,5 mod 8.
466: try_2num() will exercise multiple powers of 2 in the numerator. */
467: { "2", "1", 1 },
468: { "2", "3", -1 },
469: { "2", "5", -1 },
470: { "2", "7", 1 },
471: { "2", "9", 1 },
472: { "2", "11", -1 },
473: { "2", "13", -1 },
474: { "2", "15", 1 },
475: { "2", "17", 1 },
476:
477: /* (-2/b) = (-1)^((b^2-1)/8)*(-1)^((b-1)/2) which is -1 for b==5,7mod8.
478: try_2num() will exercise multiple powers of 2 in the numerator, which
479: will test that the shift in mpz_si_kronecker() uses unsigned not
480: signed. */
481: { "-2", "1", 1 },
482: { "-2", "3", 1 },
483: { "-2", "5", -1 },
484: { "-2", "7", -1 },
485: { "-2", "9", 1 },
486: { "-2", "11", 1 },
487: { "-2", "13", -1 },
488: { "-2", "15", -1 },
489: { "-2", "17", 1 },
490:
491: /* (a/2)=(2/a).
492: try_2den() will exercise multiple powers of 2 in the denominator. */
493: { "3", "2", -1 },
494: { "5", "2", -1 },
495: { "7", "2", 1 },
496: { "9", "2", 1 },
497: { "11", "2", -1 },
498:
499: /* Harriet Griffin, "Elementary Theory of Numbers", page 155, various
500: examples. */
501: { "2", "135", 1 },
502: { "135", "19", -1 },
503: { "2", "19", -1 },
504: { "19", "135", 1 },
505: { "173", "135", 1 },
506: { "38", "135", 1 },
507: { "135", "173", 1 },
508: { "173", "5", -1 },
509: { "3", "5", -1 },
510: { "5", "173", -1 },
511: { "173", "3", -1 },
512: { "2", "3", -1 },
513: { "3", "173", -1 },
514: { "253", "21", 1 },
515: { "1", "21", 1 },
516: { "21", "253", 1 },
517: { "21", "11", -1 },
518: { "-1", "11", -1 },
519:
520: /* Griffin page 147 */
521: { "-1", "17", 1 },
522: { "2", "17", 1 },
523: { "-2", "17", 1 },
524: { "-1", "89", 1 },
525: { "2", "89", 1 },
526:
527: /* Griffin page 148 */
528: { "89", "11", 1 },
529: { "1", "11", 1 },
530: { "89", "3", -1 },
531: { "2", "3", -1 },
532: { "3", "89", -1 },
533: { "11", "89", 1 },
534: { "33", "89", -1 },
535:
536: /* H. Davenport, "The Higher Arithmetic", page 65, the quadratic
537: residues and non-residues mod 19. */
538: { "1", "19", 1 },
539: { "4", "19", 1 },
540: { "5", "19", 1 },
541: { "6", "19", 1 },
542: { "7", "19", 1 },
543: { "9", "19", 1 },
544: { "11", "19", 1 },
545: { "16", "19", 1 },
546: { "17", "19", 1 },
547: { "2", "19", -1 },
548: { "3", "19", -1 },
549: { "8", "19", -1 },
550: { "10", "19", -1 },
551: { "12", "19", -1 },
552: { "13", "19", -1 },
553: { "14", "19", -1 },
554: { "15", "19", -1 },
555: { "18", "19", -1 },
556:
557: /* Residues and non-residues mod 13 */
558: { "0", "13", 0 },
559: { "1", "13", 1 },
560: { "2", "13", -1 },
561: { "3", "13", 1 },
562: { "4", "13", 1 },
563: { "5", "13", -1 },
564: { "6", "13", -1 },
565: { "7", "13", -1 },
566: { "8", "13", -1 },
567: { "9", "13", 1 },
568: { "10", "13", 1 },
569: { "11", "13", -1 },
570: { "12", "13", 1 },
571:
572: /* various */
573: { "5", "7", -1 },
574: { "15", "17", 1 },
575: { "67", "89", 1 },
576:
577: };
578:
579: int i, answer;
580: mpz_t a, b;
581:
582: mpz_init (a);
583: mpz_init (b);
584:
585: for (i = 0; i < numberof (data); i++)
586: {
587: mpz_set_str (a, data[i].a, 0);
588: mpz_set_str (b, data[i].b, 0);
589:
590: answer = data[i].answer;
591: try_all (a, b, data[i].answer);
592: }
593:
594: mpz_clear (a);
595: mpz_clear (b);
596: }
597:
598:
599: /* (a^2/b)=1 if gcd(a,b)=1, or (a^2/b)=0 if gcd(a,b)!=1. */
600: void
601: check_squares_zi (void)
602: {
603: mpz_t a, b, g;
604: int i, answer;
605:
606: mpz_init (a);
607: mpz_init (b);
608: mpz_init (g);
609:
610: for (i = 0; i < 200; i++)
611: {
612: mpz_random2 (a, rand() % 50);
613: mpz_random2 (b, rand() % 50);
614: mpz_mul (a, a, b);
615:
616: mpz_gcd (g, a, b);
617: if (mpz_cmp_ui (g, 1) == 0)
618: answer = 1;
619: else
620: answer = 0;
621:
622: try_all (a, b, answer);
623: }
624:
625: mpz_clear (a);
626: mpz_clear (b);
627: mpz_clear (g);
628: }
629:
630:
631: int
632: main (int argc, char *argv[])
633: {
634: if (argc >= 2 && strcmp (argv[1], "-p") == 0)
635: {
636: option_pari = 1;
637:
638: printf ("\
639: try(a,b,answer) =\n\
640: {\n\
641: if (kronecker(a,b) != answer,\n\
642: print(\"wrong at \", a, \",\", b,\n\
643: \" expected \", answer,\n\
644: \" pari says \", kronecker(a,b)))\n\
645: }\n");
646: }
647:
648: check_data ();
649: check_squares_zi ();
650:
651: exit (0);
652: }
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 / mpz / tests