1 /* Exercise mpz_*_kronecker_*() and mpz_jacobi() functions.
3 Copyright 1999, 2000, 2001, 2002, 2003, 2004 Free Software Foundation, Inc.
5 This file is part of the GNU MP Library.
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.
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.
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/. */
21 /* With no arguments the various Kronecker/Jacobi symbol routines are
22 checked against some test data and a lot of derived data.
24 To check the test data against PARI-GP, run
28 It takes a while because the output from "t-jac -p" is big.
33 More big test cases than those given by check_squares_zi would be good. */
45 #ifdef _LONG_LONG_LIMB
56 mpz_mod4 (mpz_srcptr z)
62 mpz_fdiv_r_2exp (m, z, 2);
69 mpz_fits_ulimb_p (mpz_srcptr z)
71 return (SIZ(z) == 1 || SIZ(z) == 0);
75 mpz_get_ulimb (mpz_srcptr z)
85 try_base (mp_limb_t a, mp_limb_t b, int answer)
89 if ((b & 1) == 0 || b == 1 || a > b)
92 got = mpn_jacobi_base (a, b, 0);
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"),
104 try_zi_ui (mpz_srcptr a, unsigned long b, int answer)
108 got = mpz_kronecker_ui (a, b);
111 printf ("mpz_kronecker_ui (");
112 mpz_out_str (stdout, 10, a);
113 printf (", %lu) is %d should be %d\n", b, got, answer);
120 try_zi_si (mpz_srcptr a, long b, int answer)
124 got = mpz_kronecker_si (a, b);
127 printf ("mpz_kronecker_si (");
128 mpz_out_str (stdout, 10, a);
129 printf (", %ld) is %d should be %d\n", b, got, answer);
136 try_ui_zi (unsigned long a, mpz_srcptr b, int answer)
140 got = mpz_ui_kronecker (a, b);
143 printf ("mpz_ui_kronecker (%lu, ", a);
144 mpz_out_str (stdout, 10, b);
145 printf (") is %d should be %d\n", got, answer);
152 try_si_zi (long a, mpz_srcptr b, int answer)
156 got = mpz_si_kronecker (a, b);
159 printf ("mpz_si_kronecker (%ld, ", a);
160 mpz_out_str (stdout, 10, b);
161 printf (") is %d should be %d\n", got, answer);
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. */
171 try_zi_zi (mpz_srcptr a, mpz_srcptr b, int answer)
175 got = mpz_kronecker (a, b);
178 printf ("mpz_kronecker (");
179 mpz_out_str (stdout, 10, a);
181 mpz_out_str (stdout, 10, b);
182 printf (") is %d should be %d\n", got, answer);
189 try_pari (mpz_srcptr a, mpz_srcptr b, int answer)
192 mpz_out_str (stdout, 10, a);
194 mpz_out_str (stdout, 10, b);
195 printf (",%d)\n", answer);
200 try_each (mpz_srcptr a, mpz_srcptr b, int answer)
204 try_pari (a, b, answer);
208 if (mpz_fits_ulimb_p (a) && mpz_fits_ulimb_p (b))
209 try_base (mpz_get_ulimb (a), mpz_get_ulimb (b), answer);
211 if (mpz_fits_ulong_p (b))
212 try_zi_ui (a, mpz_get_ui (b), answer);
214 if (mpz_fits_slong_p (b))
215 try_zi_si (a, mpz_get_si (b), answer);
217 if (mpz_fits_ulong_p (a))
218 try_ui_zi (mpz_get_ui (a), b, answer);
220 if (mpz_fits_sint_p (a))
221 try_si_zi (mpz_get_si (a), b, answer);
223 try_zi_zi (a, b, answer);
227 /* Try (a/b) and (a/-b). */
229 try_pn (mpz_srcptr a, mpz_srcptr b_orig, int answer)
233 mpz_init_set (b, b_orig);
234 try_each (a, b, answer);
240 try_each (a, b, answer);
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. */
250 try_periodic_num (mpz_srcptr a_orig, mpz_srcptr b, int answer)
255 if (mpz_sgn (b) <= 0)
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);
263 /* don't bother with these tests if they're only going to produce
265 if (mpz_even_p (a) && mpz_even_p (b) && mpz_even_p (a_period))
268 for (i = 0; i < 6; i++)
270 mpz_add (a, a, a_period);
271 try_pn (a, b, answer);
275 for (i = 0; i < 6; i++)
277 mpz_sub (a, a, a_period);
278 try_pn (a, b, answer);
283 mpz_clear (a_period);
287 /* Try (a/b+k*p) for various k, using the fact (a/b) is periodic in b of
293 a==3mod4 and b odd 4*a
294 a==3mod4 and b even 8*a
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. */
303 try_periodic_den (mpz_srcptr a, mpz_srcptr b_orig, int answer)
308 if (mpz_sgn (a) == 0 || mpz_sgn (b_orig) == 0)
311 mpz_init_set (b, b_orig);
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);
319 /* don't bother with these tests if they're only going to produce
321 if (mpz_even_p (a) && mpz_even_p (b) && mpz_even_p (b_period))
324 for (i = 0; i < 6; i++)
326 mpz_add (b, b, b_period);
327 try_pn (a, b, answer);
331 for (i = 0; i < 6; i++)
333 mpz_sub (b, b, b_period);
334 try_pn (a, b, answer);
339 mpz_clear (b_period);
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
351 /* Try (a/b*2^k) for various k. */
353 try_2den (mpz_srcptr a, mpz_srcptr b_orig, int answer)
357 int answer_a2, answer_k;
360 /* don't bother when b==0 */
361 if (mpz_sgn (b_orig) == 0)
364 mpz_init_set (b, b_orig);
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
371 for (kindex = 0; kindex < numberof (ktable); kindex++)
375 /* answer_k = answer*(answer_a2^k) */
376 answer_k = (answer_a2 == 0 && k != 0 ? 0
377 : (k & 1) == 1 && answer_a2 == -1 ? -answer
380 mpz_mul_2exp (b, b_orig, k);
381 try_pn (a, b, answer_k);
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. */
391 try_2num (mpz_srcptr a_orig, mpz_srcptr b, int answer)
395 int answer_2b, answer_k;
398 /* don't bother when a==0 */
399 if (mpz_sgn (a_orig) == 0)
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
409 for (kindex = 0; kindex < numberof (ktable); kindex++)
413 /* answer_k = answer*(answer_2b^k) */
414 answer_k = (answer_2b == 0 && k != 0 ? 0
415 : (k & 1) == 1 && answer_2b == -1 ? -answer
418 mpz_mul_2exp (a, a_orig, k);
419 try_pn (a, b, answer_k);
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.
430 Useful extra derived cases can be added here. */
433 try_all (mpz_t a, mpz_t b, int answer)
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);
446 static const struct {
453 /* Note that the various derived checks in try_all() reduce the cases
454 that need to be given here. */
464 In particular note (0/1)=1 so that (a/b)=(a mod b/b). */
472 /* (0/b) = 0, b != 1 */
489 /* (-1/b) = (-1)^((b-1)/2) which is -1 for b==3 mod 4 */
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. */
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
528 try_2den() will exercise multiple powers of 2 in the denominator. */
535 /* Harriet Griffin, "Elementary Theory of Numbers", page 155, various
556 /* Griffin page 147 */
563 /* Griffin page 148 */
572 /* H. Davenport, "The Higher Arithmetic", page 65, the quadratic
573 residues and non-residues mod 19. */
593 /* Residues and non-residues mod 13 */
613 /* special values inducing a==b==1 at the end of jac_or_kron() */
614 { "0x10000000000000000000000000000000000000000000000001",
615 "0x10000000000000000000000000000000000000000000000003", 1 },
624 for (i = 0; i < numberof (data); i++)
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);
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. */
639 check_squares_zi (void)
641 gmp_randstate_ptr rands = RANDS;
644 mp_size_t size_range, an, bn;
652 for (i = 0; i < 50; i++)
654 mpz_urandomb (bs, rands, 32);
655 size_range = mpz_get_ui (bs) % 10 + 2;
657 mpz_urandomb (bs, rands, size_range);
658 an = mpz_get_ui (bs);
659 mpz_rrandomb (a, rands, an);
661 mpz_urandomb (bs, rands, size_range);
662 bn = mpz_get_ui (bs);
663 mpz_rrandomb (b, rands, bn);
666 if (mpz_cmp_ui (g, 1L) == 0)
673 try_all (a, b, answer);
683 /* Check the handling of asize==0, make sure it isn't affected by the low
690 mpz_init_set_ui (a, 0);
695 try_all (a, b, 1); /* (0/1)=1 */
697 try_all (a, b, 1); /* (0/1)=1 */
701 try_all (a, b, 1); /* (0/-1)=1 */
703 try_all (a, b, 1); /* (0/-1)=1 */
707 try_all (a, b, 0); /* (0/0)=0 */
709 try_all (a, b, 0); /* (0/0)=0 */
713 try_all (a, b, 0); /* (0/2)=0 */
715 try_all (a, b, 0); /* (0/2)=0 */
723 main (int argc, char *argv[])
727 if (argc >= 2 && strcmp (argv[1], "-p") == 0)
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\