1/* Exercise mpz_*_kronecker_*() and mpz_jacobi() functions. 2 3Copyright 1999-2004, 2013 Free Software Foundation, Inc. 4 5This file is part of the GNU MP Library test suite. 6 7The GNU MP Library test suite is free software; you can redistribute it 8and/or modify it under the terms of the GNU General Public License as 9published by the Free Software Foundation; either version 3 of the License, 10or (at your option) any later version. 11 12The GNU MP Library test suite is distributed in the hope that it will be 13useful, but WITHOUT ANY WARRANTY; without even the implied warranty of 14MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General 15Public License for more details. 16 17You should have received a copy of the GNU General Public License along with 18the 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 47int option_pari = 0; 48 49 50unsigned long 51mpz_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 63int 64mpz_fits_ulimb_p (mpz_srcptr z) 65{ 66 return (SIZ(z) == 1 || SIZ(z) == 0); 67} 68 69mp_limb_t 70mpz_get_ulimb (mpz_srcptr z) 71{ 72 if (SIZ(z) == 0) 73 return 0; 74 else 75 return PTR(z)[0]; 76} 77 78 79void 80try_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 98void 99try_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 114void 115try_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 130void 131try_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 146void 147try_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. */ 165void 166try_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 183void 184try_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 194void 195try_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). */ 227void 228try_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 248void 249try_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 301void 302try_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 342static 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. */ 351void 352try_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. */ 389void 390try_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 431void 432try_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 442void 443check_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. */ 664void 665check_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. */ 711void 712check_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 */ 749int 750ref_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 768void 769check_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 872void 873check_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 982int 983main (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 ("\ 992try(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