tests/mpz/t-jac.c

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