1
/* Exercise mpz_*_kronecker_*() and mpz_jacobi() functions. */
4
Copyright 1999, 2000, 2001 Free Software Foundation, Inc.
6
This file is part of the GNU MP Library.
8
The GNU MP Library is free software; you can redistribute it and/or modify
9
it under the terms of the GNU Lesser General Public License as published by
10
the Free Software Foundation; either version 2.1 of the License, or (at your
11
option) any later version.
13
The GNU MP Library is distributed in the hope that it will be useful, but
14
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16
License for more details.
18
You should have received a copy of the GNU Lesser General Public License
19
along with the GNU MP Library; see the file COPYING.LIB. If not, write to
20
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
25
/* With no arguments the various Kronecker/Jacobi symbol routines are
26
checked against some test data and a lot of derived data.
28
To check the test data against PARI-GP, run
32
It takes a while because the output from "t-jac -p" is big.
37
More big test cases than those given by check_squares_zi would be good. */
49
#ifdef _LONG_LONG_LIMB
60
mpz_mod4 (mpz_srcptr z)
66
mpz_fdiv_r_2exp (m, z, 2);
73
mpz_fits_ulimb_p (mpz_srcptr z)
75
return (SIZ(z) == 1 || SIZ(z) == 0);
79
mpz_get_ulimb (mpz_srcptr z)
89
try_base (mp_limb_t a, mp_limb_t b, int answer)
93
if ((b & 1) == 0 || b == 1 || a > b)
96
got = mpn_jacobi_base (a, b, 0);
99
printf (LL("mpn_jacobi_base (%lu, %lu) is %d should be %d\n",
100
"mpn_jacobi_base (%llu, %llu) is %d should be %d\n"),
108
try_zi_ui (mpz_srcptr a, unsigned long b, int answer)
112
got = mpz_kronecker_ui (a, b);
115
printf ("mpz_kronecker_ui (");
116
mpz_out_str (stdout, 10, a);
117
printf (", %lu) is %d should be %d\n", b, got, answer);
124
try_zi_si (mpz_srcptr a, long b, int answer)
128
got = mpz_kronecker_si (a, b);
131
printf ("mpz_kronecker_si (");
132
mpz_out_str (stdout, 10, a);
133
printf (", %ld) is %d should be %d\n", b, got, answer);
140
try_ui_zi (unsigned long a, mpz_srcptr b, int answer)
144
got = mpz_ui_kronecker (a, b);
147
printf ("mpz_ui_kronecker (%lu, ", a);
148
mpz_out_str (stdout, 10, b);
149
printf (") is %d should be %d\n", got, answer);
156
try_si_zi (int a, mpz_srcptr b, int answer)
160
got = mpz_si_kronecker (a, b);
163
printf ("mpz_si_kronecker (%d, ", a);
164
mpz_out_str (stdout, 10, b);
165
printf (") is %d should be %d\n", got, answer);
171
/* Don't bother checking mpz_jacobi, since it only differs for b even, and
172
we don't have an actual expected answer for it. tests/devel/try.c does
173
some checks though. */
175
try_zi_zi (mpz_srcptr a, mpz_srcptr b, int answer)
179
got = mpz_kronecker (a, b);
182
printf ("mpz_kronecker (");
183
mpz_out_str (stdout, 10, a);
185
mpz_out_str (stdout, 10, b);
186
printf (") is %d should be %d\n", got, answer);
193
try_pari (mpz_srcptr a, mpz_srcptr b, int answer)
196
mpz_out_str (stdout, 10, a);
198
mpz_out_str (stdout, 10, b);
199
printf (",%d)\n", answer);
204
try_each (mpz_srcptr a, mpz_srcptr b, int answer)
208
try_pari (a, b, answer);
212
if (mpz_fits_ulimb_p (a) && mpz_fits_ulimb_p (b))
213
try_base (mpz_get_ulimb (a), mpz_get_ulimb (b), answer);
215
if (mpz_fits_ulong_p (b))
216
try_zi_ui (a, mpz_get_ui (b), answer);
218
if (mpz_fits_slong_p (b))
219
try_zi_si (a, mpz_get_si (b), answer);
221
if (mpz_fits_ulong_p (a))
222
try_ui_zi (mpz_get_ui (a), b, answer);
224
if (mpz_fits_sint_p (a))
225
try_si_zi (mpz_get_si (a), b, answer);
227
try_zi_zi (a, b, answer);
231
/* Try (a/b) and (a/-b). */
233
try_pn (mpz_srcptr a, mpz_srcptr b_orig, int answer)
237
mpz_init_set (b, b_orig);
238
try_each (a, b, answer);
244
try_each (a, b, answer);
250
/* Try (a+k*p/b) for various k, using the fact (a/b) is periodic in a with
251
period p. For b>0, p=b if b!=2mod4 or p=4*b if b==2mod4. */
254
try_periodic_num (mpz_srcptr a_orig, mpz_srcptr b, int answer)
259
if (mpz_sgn (b) <= 0)
262
mpz_init_set (a, a_orig);
263
mpz_init_set (a_period, b);
264
if (mpz_mod4 (b) == 2)
265
mpz_mul_ui (a_period, a_period, 4);
267
/* don't bother with these tests if they're only going to produce
269
if (mpz_even_p (a) && mpz_even_p (b) && mpz_even_p (a_period))
272
for (i = 0; i < 10; i++)
274
mpz_add (a, a, a_period);
275
try_pn (a, b, answer);
279
for (i = 0; i < 10; i++)
281
mpz_sub (a, a, a_period);
282
try_pn (a, b, answer);
287
mpz_clear (a_period);
291
/* Try (a/b+k*p) for various k, using the fact (a/b) is periodic in b of
297
a==3mod4 and b odd 4*a
298
a==3mod4 and b even 8*a
300
In Henri Cohen's book the period is given as 4*a for all a==2,3mod4, but
301
a counterexample would seem to be (3/2)=-1 which with (3/14)=+1 doesn't
302
have period 4*a (but rather 8*a with (3/26)=-1). Maybe the plain 4*a is
303
to be read as applying to a plain Jacobi symbol with b odd, rather than
304
the Kronecker extension to b even. */
307
try_periodic_den (mpz_srcptr a, mpz_srcptr b_orig, int answer)
312
if (mpz_sgn (a) == 0 || mpz_sgn (b_orig) == 0)
315
mpz_init_set (b, b_orig);
317
mpz_init_set (b_period, a);
318
if (mpz_mod4 (a) == 3 && mpz_even_p (b))
319
mpz_mul_ui (b_period, b_period, 8);
320
else if (mpz_mod4 (a) >= 2)
321
mpz_mul_ui (b_period, b_period, 4);
323
/* don't bother with these tests if they're only going to produce
325
if (mpz_even_p (a) && mpz_even_p (b) && mpz_even_p (b_period))
328
for (i = 0; i < 10; i++)
330
mpz_add (b, b, b_period);
331
try_pn (a, b, answer);
335
for (i = 0; i < 10; i++)
337
mpz_sub (b, b, b_period);
338
try_pn (a, b, answer);
343
mpz_clear (b_period);
347
/* Try (a/b*2^k) for various k. */
349
try_2den (mpz_srcptr a, mpz_srcptr b_orig, int answer)
355
/* don't bother when b==0 */
356
if (mpz_sgn (b_orig) == 0)
359
mpz_init_set (b, b_orig);
361
/* answer_two = (a/2) */
362
switch (mpz_fdiv_ui (a, 8)) {
377
for (i = 0; i < 3 * BITS_PER_MP_LIMB; i++)
379
answer *= answer_two;
380
mpz_mul_2exp (b, b, 1);
381
try_pn (a, b, answer);
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)
397
/* don't bother when a==0 */
398
if (mpz_sgn (a_orig) == 0)
401
mpz_init_set (a, a_orig);
402
answer_twos = mpz_ui_kronecker (2, b);
404
for (i = 0; i < 3 * BITS_PER_MP_LIMB; i++)
406
answer *= answer_twos;
407
mpz_mul_2exp (a, a, 1);
408
try_pn (a, b, answer);
415
/* The try_2num() and try_2den() routines don't in turn call
416
try_periodic_num() and try_periodic_den() because it hugely increases the
417
number of tests performed, without obviously increasing coverage.
419
Useful extra derived cases can be added here. */
422
try_all (mpz_t a, mpz_t b, int answer)
424
try_pn (a, b, answer);
425
try_periodic_num (a, b, answer);
426
try_periodic_den (a, b, answer);
427
try_2num (a, b, answer);
428
try_2den (a, b, answer);
435
static const struct {
442
/* Note that the various derived checks in try_all() reduce the cases
443
that need to be given here. */
453
In particular note (0/1)=1 so that (a/b)=(a mod b/b). */
461
/* (0/b) = 0, b != 1 */
478
/* (-1/b) = (-1)^((b-1)/2) which is -1 for b==3 mod 4 */
490
/* (2/b) = (-1)^((b^2-1)/8) which is -1 for b==3,5 mod 8.
491
try_2num() will exercise multiple powers of 2 in the numerator. */
502
/* (-2/b) = (-1)^((b^2-1)/8)*(-1)^((b-1)/2) which is -1 for b==5,7mod8.
503
try_2num() will exercise multiple powers of 2 in the numerator, which
504
will test that the shift in mpz_si_kronecker() uses unsigned not
517
try_2den() will exercise multiple powers of 2 in the denominator. */
524
/* Harriet Griffin, "Elementary Theory of Numbers", page 155, various
545
/* Griffin page 147 */
552
/* Griffin page 148 */
561
/* H. Davenport, "The Higher Arithmetic", page 65, the quadratic
562
residues and non-residues mod 19. */
582
/* Residues and non-residues mod 13 */
602
/* special values inducing a==b==1 at the end of jac_or_kron() */
603
{ "0x10000000000000000000000000000000000000000000000001",
604
"0x10000000000000000000000000000000000000000000000003", 1 },
613
for (i = 0; i < numberof (data); i++)
615
mpz_set_str_or_abort (a, data[i].a, 0);
616
mpz_set_str_or_abort (b, data[i].b, 0);
618
answer = data[i].answer;
619
try_all (a, b, data[i].answer);
627
/* (a^2/b)=1 if gcd(a,b)=1, or (a^2/b)=0 if gcd(a,b)!=1.
628
This includes when a=0 or b=0. */
630
check_squares_zi (void)
632
gmp_randstate_ptr rands = RANDS;
635
mp_size_t size_range, an, bn;
643
for (i = 0; i < 200; i++)
645
mpz_urandomb (bs, rands, 32);
646
size_range = mpz_get_ui (bs) % 10 + 2;
648
mpz_urandomb (bs, rands, size_range);
649
an = mpz_get_ui (bs);
650
mpz_rrandomb (a, rands, an);
652
mpz_urandomb (bs, rands, size_range);
653
bn = mpz_get_ui (bs);
654
mpz_rrandomb (b, rands, bn);
657
if (mpz_cmp_ui (g, 1L) == 0)
664
try_all (a, b, answer);
674
/* Check the handling of asize==0, make sure it isn't affected by the low
681
mpz_init_set_ui (a, 0);
686
try_all (a, b, 1); /* (0/1)=1 */
688
try_all (a, b, 1); /* (0/1)=1 */
692
try_all (a, b, 1); /* (0/-1)=1 */
694
try_all (a, b, 1); /* (0/-1)=1 */
698
try_all (a, b, 0); /* (0/0)=0 */
700
try_all (a, b, 0); /* (0/0)=0 */
704
try_all (a, b, 0); /* (0/2)=0 */
706
try_all (a, b, 0); /* (0/2)=0 */
714
main (int argc, char *argv[])
718
if (argc >= 2 && strcmp (argv[1], "-p") == 0)
725
if (kronecker(a,b) != answer,\n\
726
print(\"wrong at \", a, \",\", b,\n\
727
\" expected \", answer,\n\
728
\" pari says \", kronecker(a,b)))\n\
added
removed
gmp3/tests/mpz/t-jac.c
