1 /* Exercise mpz_*_kronecker_*() and mpz_jacobi() functions. 2 3 Copyright 1999-2004, 2013 Free Software Foundation, Inc. 4 5 This file is part of the GNU MP Library test suite. 6 7 The GNU MP Library test suite is free software; you can redistribute it 8 and/or modify it under the terms of the GNU General Public License as 9 published by the Free Software Foundation; either version 3 of the License, 10 or (at your option) any later version. 11 12 The GNU MP Library test suite is distributed in the hope that it will be 13 useful, but WITHOUT ANY WARRANTY; without even the implied warranty of 14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General 15 Public License for more details. 16 17 You should have received a copy of the GNU General Public License along with 18 the GNU MP Library test suite. If not, see https://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 Enhancements: 29 30 More big test cases than those given by check_squares_zi would be good. */ 31 32 33 #include <stdio.h> 34 #include <stdlib.h> 35 #include <string.h> 36 37 #include "gmp-impl.h" 38 #include "tests.h" 39 40 #ifdef _LONG_LONG_LIMB 41 #define LL(l,ll) ll 42 #else 43 #define LL(l,ll) l 44 #endif 45 46 47 int option_pari = 0; 48 49 50 unsigned long 51 mpz_mod4 (mpz_srcptr z) 52 { 53 mpz_t m; 54 unsigned long ret; 55 56 mpz_init (m); 57 mpz_fdiv_r_2exp (m, z, 2); 58 ret = mpz_get_ui (m); 59 mpz_clear (m); 60 return ret; 61 } 62 63 int 64 mpz_fits_ulimb_p (mpz_srcptr z) 65 { 66 return (SIZ(z) == 1 || SIZ(z) == 0); 67 } 68 69 mp_limb_t 70 mpz_get_ulimb (mpz_srcptr z) 71 { 72 if (SIZ(z) == 0) 73 return 0; 74 else 75 return PTR(z)[0]; 76 } 77 78 79 void 80 try_base (mp_limb_t a, mp_limb_t b, int answer) 81 { 82 int got; 83 84 if ((b & 1) == 0 || b == 1 || a > b) 85 return; 86 87 got = mpn_jacobi_base (a, b, 0); 88 if (got != answer) 89 { 90 printf (LL("mpn_jacobi_base (%lu, %lu) is %d should be %d\n", 91 "mpn_jacobi_base (%llu, %llu) is %d should be %d\n"), 92 a, b, got, answer); 93 abort (); 94 } 95 } 96 97 98 void 99 try_zi_ui (mpz_srcptr a, unsigned long b, int answer) 100 { 101 int got; 102 103 got = mpz_kronecker_ui (a, b); 104 if (got != answer) 105 { 106 printf ("mpz_kronecker_ui ("); 107 mpz_out_str (stdout, 10, a); 108 printf (", %lu) is %d should be %d\n", b, got, answer); 109 abort (); 110 } 111 } 112 113 114 void 115 try_zi_si (mpz_srcptr a, long b, int answer) 116 { 117 int got; 118 119 got = mpz_kronecker_si (a, b); 120 if (got != answer) 121 { 122 printf ("mpz_kronecker_si ("); 123 mpz_out_str (stdout, 10, a); 124 printf (", %ld) is %d should be %d\n", b, got, answer); 125 abort (); 126 } 127 } 128 129 130 void 131 try_ui_zi (unsigned long a, mpz_srcptr b, int answer) 132 { 133 int got; 134 135 got = mpz_ui_kronecker (a, b); 136 if (got != answer) 137 { 138 printf ("mpz_ui_kronecker (%lu, ", a); 139 mpz_out_str (stdout, 10, b); 140 printf (") is %d should be %d\n", got, answer); 141 abort (); 142 } 143 } 144 145 146 void 147 try_si_zi (long a, mpz_srcptr b, int answer) 148 { 149 int got; 150 151 got = mpz_si_kronecker (a, b); 152 if (got != answer) 153 { 154 printf ("mpz_si_kronecker (%ld, ", a); 155 mpz_out_str (stdout, 10, b); 156 printf (") is %d should be %d\n", got, answer); 157 abort (); 158 } 159 } 160 161 162 /* Don't bother checking mpz_jacobi, since it only differs for b even, and 163 we don't have an actual expected answer for it. tests/devel/try.c does 164 some checks though. */ 165 void 166 try_zi_zi (mpz_srcptr a, mpz_srcptr b, int answer) 167 { 168 int got; 169 170 got = mpz_kronecker (a, b); 171 if (got != answer) 172 { 173 printf ("mpz_kronecker ("); 174 mpz_out_str (stdout, 10, a); 175 printf (", "); 176 mpz_out_str (stdout, 10, b); 177 printf (") is %d should be %d\n", got, answer); 178 abort (); 179 } 180 } 181 182 183 void 184 try_pari (mpz_srcptr a, mpz_srcptr b, int answer) 185 { 186 printf ("try("); 187 mpz_out_str (stdout, 10, a); 188 printf (","); 189 mpz_out_str (stdout, 10, b); 190 printf (",%d)\n", answer); 191 } 192 193 194 void 195 try_each (mpz_srcptr a, mpz_srcptr b, int answer) 196 { 197 #if 0 198 fprintf(stderr, "asize = %d, bsize = %d\n", 199 mpz_sizeinbase (a, 2), mpz_sizeinbase (b, 2)); 200 #endif 201 if (option_pari) 202 { 203 try_pari (a, b, answer); 204 return; 205 } 206 207 if (mpz_fits_ulimb_p (a) && mpz_fits_ulimb_p (b)) 208 try_base (mpz_get_ulimb (a), mpz_get_ulimb (b), answer); 209 210 if (mpz_fits_ulong_p (b)) 211 try_zi_ui (a, mpz_get_ui (b), answer); 212 213 if (mpz_fits_slong_p (b)) 214 try_zi_si (a, mpz_get_si (b), answer); 215 216 if (mpz_fits_ulong_p (a)) 217 try_ui_zi (mpz_get_ui (a), b, answer); 218 219 if (mpz_fits_sint_p (a)) 220 try_si_zi (mpz_get_si (a), b, answer); 221 222 try_zi_zi (a, b, answer); 223 } 224 225 226 /* Try (a/b) and (a/-b). */ 227 void 228 try_pn (mpz_srcptr a, mpz_srcptr b_orig, int answer) 229 { 230 mpz_t b; 231 232 mpz_init_set (b, b_orig); 233 try_each (a, b, answer); 234 235 mpz_neg (b, b); 236 if (mpz_sgn (a) < 0) 237 answer = -answer; 238 239 try_each (a, b, answer); 240 241 mpz_clear (b); 242 } 243 244 245 /* Try (a+k*p/b) for various k, using the fact (a/b) is periodic in a with 246 period p. For b>0, p=b if b!=2mod4 or p=4*b if b==2mod4. */ 247 248 void 249 try_periodic_num (mpz_srcptr a_orig, mpz_srcptr b, int answer) 250 { 251 mpz_t a, a_period; 252 int i; 253 254 if (mpz_sgn (b) <= 0) 255 return; 256 257 mpz_init_set (a, a_orig); 258 mpz_init_set (a_period, b); 259 if (mpz_mod4 (b) == 2) 260 mpz_mul_ui (a_period, a_period, 4); 261 262 /* don't bother with these tests if they're only going to produce 263 even/even */ 264 if (mpz_even_p (a) && mpz_even_p (b) && mpz_even_p (a_period)) 265 goto done; 266 267 for (i = 0; i < 6; i++) 268 { 269 mpz_add (a, a, a_period); 270 try_pn (a, b, answer); 271 } 272 273 mpz_set (a, a_orig); 274 for (i = 0; i < 6; i++) 275 { 276 mpz_sub (a, a, a_period); 277 try_pn (a, b, answer); 278 } 279 280 done: 281 mpz_clear (a); 282 mpz_clear (a_period); 283 } 284 285 286 /* Try (a/b+k*p) for various k, using the fact (a/b) is periodic in b of 287 period p. 288 289 period p 290 a==0,1mod4 a 291 a==2mod4 4*a 292 a==3mod4 and b odd 4*a 293 a==3mod4 and b even 8*a 294 295 In Henri Cohen's book the period is given as 4*a for all a==2,3mod4, but 296 a counterexample would seem to be (3/2)=-1 which with (3/14)=+1 doesn't 297 have period 4*a (but rather 8*a with (3/26)=-1). Maybe the plain 4*a is 298 to be read as applying to a plain Jacobi symbol with b odd, rather than 299 the Kronecker extension to b even. */ 300 301 void 302 try_periodic_den (mpz_srcptr a, mpz_srcptr b_orig, int answer) 303 { 304 mpz_t b, b_period; 305 int i; 306 307 if (mpz_sgn (a) == 0 || mpz_sgn (b_orig) == 0) 308 return; 309 310 mpz_init_set (b, b_orig); 311 312 mpz_init_set (b_period, a); 313 if (mpz_mod4 (a) == 3 && mpz_even_p (b)) 314 mpz_mul_ui (b_period, b_period, 8L); 315 else if (mpz_mod4 (a) >= 2) 316 mpz_mul_ui (b_period, b_period, 4L); 317 318 /* don't bother with these tests if they're only going to produce 319 even/even */ 320 if (mpz_even_p (a) && mpz_even_p (b) && mpz_even_p (b_period)) 321 goto done; 322 323 for (i = 0; i < 6; i++) 324 { 325 mpz_add (b, b, b_period); 326 try_pn (a, b, answer); 327 } 328 329 mpz_set (b, b_orig); 330 for (i = 0; i < 6; i++) 331 { 332 mpz_sub (b, b, b_period); 333 try_pn (a, b, answer); 334 } 335 336 done: 337 mpz_clear (b); 338 mpz_clear (b_period); 339 } 340 341 342 static const unsigned long ktable[] = { 343 0, 1, 2, 3, 4, 5, 6, 7, 344 GMP_NUMB_BITS-1, GMP_NUMB_BITS, GMP_NUMB_BITS+1, 345 2*GMP_NUMB_BITS-1, 2*GMP_NUMB_BITS, 2*GMP_NUMB_BITS+1, 346 3*GMP_NUMB_BITS-1, 3*GMP_NUMB_BITS, 3*GMP_NUMB_BITS+1 347 }; 348 349 350 /* Try (a/b*2^k) for various k. */ 351 void 352 try_2den (mpz_srcptr a, mpz_srcptr b_orig, int answer) 353 { 354 mpz_t b; 355 int kindex; 356 int answer_a2, answer_k; 357 unsigned long k; 358 359 /* don't bother when b==0 */ 360 if (mpz_sgn (b_orig) == 0) 361 return; 362 363 mpz_init_set (b, b_orig); 364 365 /* (a/2) is 0 if a even, 1 if a==1 or 7 mod 8, -1 if a==3 or 5 mod 8 */ 366 answer_a2 = (mpz_even_p (a) ? 0 367 : (((SIZ(a) >= 0 ? PTR(a)[0] : -PTR(a)[0]) + 2) & 7) < 4 ? 1 368 : -1); 369 370 for (kindex = 0; kindex < numberof (ktable); kindex++) 371 { 372 k = ktable[kindex]; 373 374 /* answer_k = answer*(answer_a2^k) */ 375 answer_k = (answer_a2 == 0 && k != 0 ? 0 376 : (k & 1) == 1 && answer_a2 == -1 ? -answer 377 : answer); 378 379 mpz_mul_2exp (b, b_orig, k); 380 try_pn (a, b, answer_k); 381 } 382 383 mpz_clear (b); 384 } 385 386 387 /* Try (a*2^k/b) for various k. If it happens mpz_ui_kronecker() gets (2/b) 388 wrong it will show up as wrong answers demanded. */ 389 void 390 try_2num (mpz_srcptr a_orig, mpz_srcptr b, int answer) 391 { 392 mpz_t a; 393 int kindex; 394 int answer_2b, answer_k; 395 unsigned long k; 396 397 /* don't bother when a==0 */ 398 if (mpz_sgn (a_orig) == 0) 399 return; 400 401 mpz_init (a); 402 403 /* (2/b) is 0 if b even, 1 if b==1 or 7 mod 8, -1 if b==3 or 5 mod 8 */ 404 answer_2b = (mpz_even_p (b) ? 0 405 : (((SIZ(b) >= 0 ? PTR(b)[0] : -PTR(b)[0]) + 2) & 7) < 4 ? 1 406 : -1); 407 408 for (kindex = 0; kindex < numberof (ktable); kindex++) 409 { 410 k = ktable[kindex]; 411 412 /* answer_k = answer*(answer_2b^k) */ 413 answer_k = (answer_2b == 0 && k != 0 ? 0 414 : (k & 1) == 1 && answer_2b == -1 ? -answer 415 : answer); 416 417 mpz_mul_2exp (a, a_orig, k); 418 try_pn (a, b, answer_k); 419 } 420 421 mpz_clear (a); 422 } 423 424 425 /* The try_2num() and try_2den() routines don't in turn call 426 try_periodic_num() and try_periodic_den() because it hugely increases the 427 number of tests performed, without obviously increasing coverage. 428 429 Useful extra derived cases can be added here. */ 430 431 void 432 try_all (mpz_t a, mpz_t b, int answer) 433 { 434 try_pn (a, b, answer); 435 try_periodic_num (a, b, answer); 436 try_periodic_den (a, b, answer); 437 try_2num (a, b, answer); 438 try_2den (a, b, answer); 439 } 440 441 442 void 443 check_data (void) 444 { 445 static const struct { 446 const char *a; 447 const char *b; 448 int answer; 449 450 } data[] = { 451 452 /* Note that the various derived checks in try_all() reduce the cases 453 that need to be given here. */ 454 455 /* some zeros */ 456 { "0", "0", 0 }, 457 { "0", "2", 0 }, 458 { "0", "6", 0 }, 459 { "5", "0", 0 }, 460 { "24", "60", 0 }, 461 462 /* (a/1) = 1, any a 463 In particular note (0/1)=1 so that (a/b)=(a mod b/b). */ 464 { "0", "1", 1 }, 465 { "1", "1", 1 }, 466 { "2", "1", 1 }, 467 { "3", "1", 1 }, 468 { "4", "1", 1 }, 469 { "5", "1", 1 }, 470 471 /* (0/b) = 0, b != 1 */ 472 { "0", "3", 0 }, 473 { "0", "5", 0 }, 474 { "0", "7", 0 }, 475 { "0", "9", 0 }, 476 { "0", "11", 0 }, 477 { "0", "13", 0 }, 478 { "0", "15", 0 }, 479 480 /* (1/b) = 1 */ 481 { "1", "1", 1 }, 482 { "1", "3", 1 }, 483 { "1", "5", 1 }, 484 { "1", "7", 1 }, 485 { "1", "9", 1 }, 486 { "1", "11", 1 }, 487 488 /* (-1/b) = (-1)^((b-1)/2) which is -1 for b==3 mod 4 */ 489 { "-1", "1", 1 }, 490 { "-1", "3", -1 }, 491 { "-1", "5", 1 }, 492 { "-1", "7", -1 }, 493 { "-1", "9", 1 }, 494 { "-1", "11", -1 }, 495 { "-1", "13", 1 }, 496 { "-1", "15", -1 }, 497 { "-1", "17", 1 }, 498 { "-1", "19", -1 }, 499 500 /* (2/b) = (-1)^((b^2-1)/8) which is -1 for b==3,5 mod 8. 501 try_2num() will exercise multiple powers of 2 in the numerator. */ 502 { "2", "1", 1 }, 503 { "2", "3", -1 }, 504 { "2", "5", -1 }, 505 { "2", "7", 1 }, 506 { "2", "9", 1 }, 507 { "2", "11", -1 }, 508 { "2", "13", -1 }, 509 { "2", "15", 1 }, 510 { "2", "17", 1 }, 511 512 /* (-2/b) = (-1)^((b^2-1)/8)*(-1)^((b-1)/2) which is -1 for b==5,7mod8. 513 try_2num() will exercise multiple powers of 2 in the numerator, which 514 will test that the shift in mpz_si_kronecker() uses unsigned not 515 signed. */ 516 { "-2", "1", 1 }, 517 { "-2", "3", 1 }, 518 { "-2", "5", -1 }, 519 { "-2", "7", -1 }, 520 { "-2", "9", 1 }, 521 { "-2", "11", 1 }, 522 { "-2", "13", -1 }, 523 { "-2", "15", -1 }, 524 { "-2", "17", 1 }, 525 526 /* (a/2)=(2/a). 527 try_2den() will exercise multiple powers of 2 in the denominator. */ 528 { "3", "2", -1 }, 529 { "5", "2", -1 }, 530 { "7", "2", 1 }, 531 { "9", "2", 1 }, 532 { "11", "2", -1 }, 533 534 /* Harriet Griffin, "Elementary Theory of Numbers", page 155, various 535 examples. */ 536 { "2", "135", 1 }, 537 { "135", "19", -1 }, 538 { "2", "19", -1 }, 539 { "19", "135", 1 }, 540 { "173", "135", 1 }, 541 { "38", "135", 1 }, 542 { "135", "173", 1 }, 543 { "173", "5", -1 }, 544 { "3", "5", -1 }, 545 { "5", "173", -1 }, 546 { "173", "3", -1 }, 547 { "2", "3", -1 }, 548 { "3", "173", -1 }, 549 { "253", "21", 1 }, 550 { "1", "21", 1 }, 551 { "21", "253", 1 }, 552 { "21", "11", -1 }, 553 { "-1", "11", -1 }, 554 555 /* Griffin page 147 */ 556 { "-1", "17", 1 }, 557 { "2", "17", 1 }, 558 { "-2", "17", 1 }, 559 { "-1", "89", 1 }, 560 { "2", "89", 1 }, 561 562 /* Griffin page 148 */ 563 { "89", "11", 1 }, 564 { "1", "11", 1 }, 565 { "89", "3", -1 }, 566 { "2", "3", -1 }, 567 { "3", "89", -1 }, 568 { "11", "89", 1 }, 569 { "33", "89", -1 }, 570 571 /* H. Davenport, "The Higher Arithmetic", page 65, the quadratic 572 residues and non-residues mod 19. */ 573 { "1", "19", 1 }, 574 { "4", "19", 1 }, 575 { "5", "19", 1 }, 576 { "6", "19", 1 }, 577 { "7", "19", 1 }, 578 { "9", "19", 1 }, 579 { "11", "19", 1 }, 580 { "16", "19", 1 }, 581 { "17", "19", 1 }, 582 { "2", "19", -1 }, 583 { "3", "19", -1 }, 584 { "8", "19", -1 }, 585 { "10", "19", -1 }, 586 { "12", "19", -1 }, 587 { "13", "19", -1 }, 588 { "14", "19", -1 }, 589 { "15", "19", -1 }, 590 { "18", "19", -1 }, 591 592 /* Residues and non-residues mod 13 */ 593 { "0", "13", 0 }, 594 { "1", "13", 1 }, 595 { "2", "13", -1 }, 596 { "3", "13", 1 }, 597 { "4", "13", 1 }, 598 { "5", "13", -1 }, 599 { "6", "13", -1 }, 600 { "7", "13", -1 }, 601 { "8", "13", -1 }, 602 { "9", "13", 1 }, 603 { "10", "13", 1 }, 604 { "11", "13", -1 }, 605 { "12", "13", 1 }, 606 607 /* various */ 608 { "5", "7", -1 }, 609 { "15", "17", 1 }, 610 { "67", "89", 1 }, 611 612 /* special values inducing a==b==1 at the end of jac_or_kron() */ 613 { "0x10000000000000000000000000000000000000000000000001", 614 "0x10000000000000000000000000000000000000000000000003", 1 }, 615 616 /* Test for previous bugs in jacobi_2. */ 617 { "0x43900000000", "0x42400000439", -1 }, /* 32-bit limbs */ 618 { "0x4390000000000000000", "0x4240000000000000439", -1 }, /* 64-bit limbs */ 619 620 { "198158408161039063", "198158360916398807", -1 }, 621 622 /* Some tests involving large quotients in the continued fraction 623 expansion. */ 624 { "37200210845139167613356125645445281805", 625 "451716845976689892447895811408978421929", -1 }, 626 { "67674091930576781943923596701346271058970643542491743605048620644676477275152701774960868941561652032482173612421015", 627 "4902678867794567120224500687210807069172039735", 0 }, 628 { "2666617146103764067061017961903284334497474492754652499788571378062969111250584288683585223600172138551198546085281683283672592", "2666617146103764067061017961903284334497474492754652499788571378062969111250584288683585223600172138551198546085281683290481773", 1 }, 629 630 /* Exercises the case asize == 1, btwos > 0 in mpz_jacobi. */ 631 { "804609", "421248363205206617296534688032638102314410556521742428832362659824", 1 } , 632 { "4190209", "2239744742177804210557442048984321017460028974602978995388383905961079286530650825925074203175536427000", 1 }, 633 634 /* Exercises the case asize == 1, btwos = 63 in mpz_jacobi 635 (relevant when GMP_LIMB_BITS == 64). */ 636 { "17311973299000934401", "1675975991242824637446753124775689449936871337036614677577044717424700351103148799107651171694863695242089956242888229458836426332300124417011114380886016", 1 }, 637 { "3220569220116583677", "41859917623035396746", -1 }, 638 639 /* Other test cases that triggered bugs during development. */ 640 { "37200210845139167613356125645445281805", "340116213441272389607827434472642576514", -1 }, 641 { "74400421690278335226712251290890563610", "451716845976689892447895811408978421929", -1 }, 642 }; 643 644 int i; 645 mpz_t a, b; 646 647 mpz_init (a); 648 mpz_init (b); 649 650 for (i = 0; i < numberof (data); i++) 651 { 652 mpz_set_str_or_abort (a, data[i].a, 0); 653 mpz_set_str_or_abort (b, data[i].b, 0); 654 try_all (a, b, data[i].answer); 655 } 656 657 mpz_clear (a); 658 mpz_clear (b); 659 } 660 661 662 /* (a^2/b)=1 if gcd(a,b)=1, or (a^2/b)=0 if gcd(a,b)!=1. 663 This includes when a=0 or b=0. */ 664 void 665 check_squares_zi (void) 666 { 667 gmp_randstate_ptr rands = RANDS; 668 mpz_t a, b, g; 669 int i, answer; 670 mp_size_t size_range, an, bn; 671 mpz_t bs; 672 673 mpz_init (bs); 674 mpz_init (a); 675 mpz_init (b); 676 mpz_init (g); 677 678 for (i = 0; i < 50; i++) 679 { 680 mpz_urandomb (bs, rands, 32); 681 size_range = mpz_get_ui (bs) % 10 + i/8 + 2; 682 683 mpz_urandomb (bs, rands, size_range); 684 an = mpz_get_ui (bs); 685 mpz_rrandomb (a, rands, an); 686 687 mpz_urandomb (bs, rands, size_range); 688 bn = mpz_get_ui (bs); 689 mpz_rrandomb (b, rands, bn); 690 691 mpz_gcd (g, a, b); 692 if (mpz_cmp_ui (g, 1L) == 0) 693 answer = 1; 694 else 695 answer = 0; 696 697 mpz_mul (a, a, a); 698 699 try_all (a, b, answer); 700 } 701 702 mpz_clear (bs); 703 mpz_clear (a); 704 mpz_clear (b); 705 mpz_clear (g); 706 } 707 708 709 /* Check the handling of asize==0, make sure it isn't affected by the low 710 limb. */ 711 void 712 check_a_zero (void) 713 { 714 mpz_t a, b; 715 716 mpz_init_set_ui (a, 0); 717 mpz_init (b); 718 719 mpz_set_ui (b, 1L); 720 PTR(a)[0] = 0; 721 try_all (a, b, 1); /* (0/1)=1 */ 722 PTR(a)[0] = 1; 723 try_all (a, b, 1); /* (0/1)=1 */ 724 725 mpz_set_si (b, -1L); 726 PTR(a)[0] = 0; 727 try_all (a, b, 1); /* (0/-1)=1 */ 728 PTR(a)[0] = 1; 729 try_all (a, b, 1); /* (0/-1)=1 */ 730 731 mpz_set_ui (b, 0); 732 PTR(a)[0] = 0; 733 try_all (a, b, 0); /* (0/0)=0 */ 734 PTR(a)[0] = 1; 735 try_all (a, b, 0); /* (0/0)=0 */ 736 737 mpz_set_ui (b, 2); 738 PTR(a)[0] = 0; 739 try_all (a, b, 0); /* (0/2)=0 */ 740 PTR(a)[0] = 1; 741 try_all (a, b, 0); /* (0/2)=0 */ 742 743 mpz_clear (a); 744 mpz_clear (b); 745 } 746 747 748 /* Assumes that b = prod p_k^e_k */ 749 int 750 ref_jacobi (mpz_srcptr a, mpz_srcptr b, unsigned nprime, 751 mpz_t prime[], unsigned *exp) 752 { 753 unsigned i; 754 int res; 755 756 for (i = 0, res = 1; i < nprime; i++) 757 if (exp[i]) 758 { 759 int legendre = refmpz_legendre (a, prime[i]); 760 if (!legendre) 761 return 0; 762 if (exp[i] & 1) 763 res *= legendre; 764 } 765 return res; 766 } 767 768 void 769 check_jacobi_factored (void) 770 { 771 #define PRIME_N 10 772 #define PRIME_MAX_SIZE 50 773 #define PRIME_MAX_EXP 4 774 #define PRIME_A_COUNT 10 775 #define PRIME_B_COUNT 5 776 #define PRIME_MAX_B_SIZE 2000 777 778 gmp_randstate_ptr rands = RANDS; 779 mpz_t prime[PRIME_N]; 780 unsigned exp[PRIME_N]; 781 mpz_t a, b, t, bs; 782 unsigned i; 783 784 mpz_init (a); 785 mpz_init (b); 786 mpz_init (t); 787 mpz_init (bs); 788 789 /* Generate primes */ 790 for (i = 0; i < PRIME_N; i++) 791 { 792 mp_size_t size; 793 mpz_init (prime[i]); 794 mpz_urandomb (bs, rands, 32); 795 size = mpz_get_ui (bs) % PRIME_MAX_SIZE + 2; 796 mpz_rrandomb (prime[i], rands, size); 797 if (mpz_cmp_ui (prime[i], 3) <= 0) 798 mpz_set_ui (prime[i], 3); 799 else 800 mpz_nextprime (prime[i], prime[i]); 801 } 802 803 for (i = 0; i < PRIME_B_COUNT; i++) 804 { 805 unsigned j, k; 806 mp_bitcnt_t bsize; 807 808 mpz_set_ui (b, 1); 809 bsize = 1; 810 811 for (j = 0; j < PRIME_N && bsize < PRIME_MAX_B_SIZE; j++) 812 { 813 mpz_urandomb (bs, rands, 32); 814 exp[j] = mpz_get_ui (bs) % PRIME_MAX_EXP; 815 mpz_pow_ui (t, prime[j], exp[j]); 816 mpz_mul (b, b, t); 817 bsize = mpz_sizeinbase (b, 2); 818 } 819 for (k = 0; k < PRIME_A_COUNT; k++) 820 { 821 int answer; 822 mpz_rrandomb (a, rands, bsize + 2); 823 answer = ref_jacobi (a, b, j, prime, exp); 824 try_all (a, b, answer); 825 } 826 } 827 for (i = 0; i < PRIME_N; i++) 828 mpz_clear (prime[i]); 829 830 mpz_clear (a); 831 mpz_clear (b); 832 mpz_clear (t); 833 mpz_clear (bs); 834 835 #undef PRIME_N 836 #undef PRIME_MAX_SIZE 837 #undef PRIME_MAX_EXP 838 #undef PRIME_A_COUNT 839 #undef PRIME_B_COUNT 840 #undef PRIME_MAX_B_SIZE 841 } 842 843 /* These tests compute (a|n), where the quotient sequence includes 844 large quotients, and n has a known factorization. Such inputs are 845 generated as follows. First, construct a large n, as a power of a 846 prime p of moderate size. 847 848 Next, compute a matrix from factors (q,1;1,0), with q chosen with 849 uniformly distributed size. We must stop with matrix elements of 850 roughly half the size of n. Denote elements of M as M = (m00, m01; 851 m10, m11). 852 853 We now look for solutions to 854 855 n = m00 x + m01 y 856 a = m10 x + m11 y 857 858 with x,y > 0. Since n >= m00 * m01, there exists a positive 859 solution to the first equation. Find those x, y, and substitute in 860 the second equation to get a. Then the quotient sequence for (a|n) 861 is precisely the quotients used when constructing M, followed by 862 the quotient sequence for (x|y). 863 864 Numbers should also be large enough that we exercise hgcd_jacobi, 865 which means that they should be larger than 866 867 max (GCD_DC_THRESHOLD, 3 * HGCD_THRESHOLD) 868 869 With an n of roughly 40000 bits, this should hold on most machines. 870 */ 871 872 void 873 check_large_quotients (void) 874 { 875 #define COUNT 50 876 #define PBITS 200 877 #define PPOWER 201 878 #define MAX_QBITS 500 879 880 gmp_randstate_ptr rands = RANDS; 881 882 mpz_t p, n, q, g, s, t, x, y, bs; 883 mpz_t M[2][2]; 884 mp_bitcnt_t nsize; 885 unsigned i; 886 887 mpz_init (p); 888 mpz_init (n); 889 mpz_init (q); 890 mpz_init (g); 891 mpz_init (s); 892 mpz_init (t); 893 mpz_init (x); 894 mpz_init (y); 895 mpz_init (bs); 896 mpz_init (M[0][0]); 897 mpz_init (M[0][1]); 898 mpz_init (M[1][0]); 899 mpz_init (M[1][1]); 900 901 /* First generate a number with known factorization, as a random 902 smallish prime raised to an odd power. Then (a|n) = (a|p). */ 903 mpz_rrandomb (p, rands, PBITS); 904 mpz_nextprime (p, p); 905 mpz_pow_ui (n, p, PPOWER); 906 907 nsize = mpz_sizeinbase (n, 2); 908 909 for (i = 0; i < COUNT; i++) 910 { 911 int answer; 912 mp_bitcnt_t msize; 913 914 mpz_set_ui (M[0][0], 1); 915 mpz_set_ui (M[0][1], 0); 916 mpz_set_ui (M[1][0], 0); 917 mpz_set_ui (M[1][1], 1); 918 919 for (msize = 1; 2*(msize + MAX_QBITS) + 1 < nsize ;) 920 { 921 unsigned i; 922 mpz_rrandomb (bs, rands, 32); 923 mpz_rrandomb (q, rands, 1 + mpz_get_ui (bs) % MAX_QBITS); 924 925 /* Multiply by (q, 1; 1,0) from the right */ 926 for (i = 0; i < 2; i++) 927 { 928 mp_bitcnt_t size; 929 mpz_swap (M[i][0], M[i][1]); 930 mpz_addmul (M[i][0], M[i][1], q); 931 size = mpz_sizeinbase (M[i][0], 2); 932 if (size > msize) 933 msize = size; 934 } 935 } 936 mpz_gcdext (g, s, t, M[0][0], M[0][1]); 937 ASSERT_ALWAYS (mpz_cmp_ui (g, 1) == 0); 938 939 /* Solve n = M[0][0] * x + M[0][1] * y */ 940 if (mpz_sgn (s) > 0) 941 { 942 mpz_mul (x, n, s); 943 mpz_fdiv_qr (q, x, x, M[0][1]); 944 mpz_mul (y, q, M[0][0]); 945 mpz_addmul (y, t, n); 946 ASSERT_ALWAYS (mpz_sgn (y) > 0); 947 } 948 else 949 { 950 mpz_mul (y, n, t); 951 mpz_fdiv_qr (q, y, y, M[0][0]); 952 mpz_mul (x, q, M[0][1]); 953 mpz_addmul (x, s, n); 954 ASSERT_ALWAYS (mpz_sgn (x) > 0); 955 } 956 mpz_mul (x, x, M[1][0]); 957 mpz_addmul (x, y, M[1][1]); 958 959 /* Now (x|n) has the selected large quotients */ 960 answer = refmpz_legendre (x, p); 961 try_zi_zi (x, n, answer); 962 } 963 mpz_clear (p); 964 mpz_clear (n); 965 mpz_clear (q); 966 mpz_clear (g); 967 mpz_clear (s); 968 mpz_clear (t); 969 mpz_clear (x); 970 mpz_clear (y); 971 mpz_clear (bs); 972 mpz_clear (M[0][0]); 973 mpz_clear (M[0][1]); 974 mpz_clear (M[1][0]); 975 mpz_clear (M[1][1]); 976 #undef COUNT 977 #undef PBITS 978 #undef PPOWER 979 #undef MAX_QBITS 980 } 981 982 int 983 main (int argc, char *argv[]) 984 { 985 tests_start (); 986 987 if (argc >= 2 && strcmp (argv[1], "-p") == 0) 988 { 989 option_pari = 1; 990 991 printf ("\ 992 try(a,b,answer) =\n\ 993 {\n\ 994 if (kronecker(a,b) != answer,\n\ 995 print(\"wrong at \", a, \",\", b,\n\ 996 \" expected \", answer,\n\ 997 \" pari says \", kronecker(a,b)))\n\ 998 }\n"); 999 } 1000 1001 check_data (); 1002 check_squares_zi (); 1003 check_a_zero (); 1004 check_jacobi_factored (); 1005 check_large_quotients (); 1006 tests_end (); 1007 exit (0); 1008 } 1009