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
mpz_mod4(mpz_srcptr z)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
mpz_fits_ulimb_p(mpz_srcptr z)68 mpz_fits_ulimb_p (mpz_srcptr z)
69 {
70 return (SIZ(z) == 1 || SIZ(z) == 0);
71 }
72
73 mp_limb_t
mpz_get_ulimb(mpz_srcptr z)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
try_base(mp_limb_t a,mp_limb_t b,int answer)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
try_zi_ui(mpz_srcptr a,unsigned long b,int answer)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
try_zi_si(mpz_srcptr a,long b,int answer)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
try_ui_zi(unsigned long a,mpz_srcptr b,int answer)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
try_si_zi(long a,mpz_srcptr b,int answer)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
try_zi_zi(mpz_srcptr a,mpz_srcptr b,int answer)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
try_pari(mpz_srcptr a,mpz_srcptr b,int answer)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
try_each(mpz_srcptr a,mpz_srcptr b,int answer)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
try_pn(mpz_srcptr a,mpz_srcptr b_orig,int answer)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
try_periodic_num(mpz_srcptr a_orig,mpz_srcptr b,int answer)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
try_periodic_den(mpz_srcptr a,mpz_srcptr b_orig,int answer)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
try_2den(mpz_srcptr a,mpz_srcptr b_orig,int answer)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
try_2num(mpz_srcptr a_orig,mpz_srcptr b,int answer)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
try_all(mpz_t a,mpz_t b,int answer)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
check_data(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
check_squares_zi(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
check_a_zero(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
ref_jacobi(mpz_srcptr a,mpz_srcptr b,unsigned nprime,mpz_t prime[],unsigned * exp)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
check_jacobi_factored(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
check_large_quotients(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
main(int argc,char * argv[])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 }
1013