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 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
47 int option_pari = 0;
48
49
50 unsigned long
51 mpz_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
63 int
64 mpz_fits_ulimb_p (mpz_srcptr z)
65 {
66 return (SIZ(z) == 1 || SIZ(z) == 0);
67 }
68
69 mp_limb_t
70 mpz_get_ulimb (mpz_srcptr z)
71 {
72 if (SIZ(z) == 0)
73 return 0;
74 else
75 return PTR(z)[0];
76 }
77
78
79 void
80 try_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
98 void
99 try_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
114 void
115 try_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
130 void
131 try_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
146 void
147 try_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. */
165 void
166 try_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
183 void
184 try_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
194 void
195 try_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). */
227 void
228 try_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
248 void
249 try_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
301 void
302 try_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
342 static 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. */
351 void
352 try_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. */
389 void
390 try_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
431 void
432 try_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
442 void
443 check_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. */
664 void
665 check_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. */
711 void
712 check_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 */
749 int
750 ref_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
768 void
769 check_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
872 void
873 check_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
982 int
983 main (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 ("\
992 try(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 }