Annotation of OpenXM_contrib/gmp/mpz/tests/t-jac.c, Revision 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