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