/* Exercise mpz_*_kronecker_*() and mpz_jacobi() functions. Copyright 1999, 2000, 2001, 2002, 2003, 2004 Free Software Foundation, Inc. This file is part of the GNU MP Library test suite. The GNU MP Library test suite is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 3 of the License, or (at your option) any later version. The GNU MP Library test suite is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with the GNU MP Library test suite. If not, see http://www.gnu.org/licenses/. */ /* With no arguments the various Kronecker/Jacobi symbol routines are checked against some test data and a lot of derived data. To check the test data against PARI-GP, run t-jac -p | gp -q It takes a while because the output from "t-jac -p" is big. Enhancements: More big test cases than those given by check_squares_zi would be good. */ #include #include #include #include "gmp.h" #include "gmp-impl.h" #include "tests.h" #ifdef _LONG_LONG_LIMB #define LL(l,ll) ll #else #define LL(l,ll) l #endif int option_pari = 0; unsigned long mpz_mod4 (mpz_srcptr z) { mpz_t m; unsigned long ret; mpz_init (m); mpz_fdiv_r_2exp (m, z, 2); ret = mpz_get_ui (m); mpz_clear (m); return ret; } int mpz_fits_ulimb_p (mpz_srcptr z) { return (SIZ(z) == 1 || SIZ(z) == 0); } mp_limb_t mpz_get_ulimb (mpz_srcptr z) { if (SIZ(z) == 0) return 0; else return PTR(z)[0]; } void try_base (mp_limb_t a, mp_limb_t b, int answer) { int got; if ((b & 1) == 0 || b == 1 || a > b) return; got = mpn_jacobi_base (a, b, 0); if (got != answer) { printf (LL("mpn_jacobi_base (%lu, %lu) is %d should be %d\n", "mpn_jacobi_base (%llu, %llu) is %d should be %d\n"), a, b, got, answer); abort (); } } void try_zi_ui (mpz_srcptr a, unsigned long b, int answer) { int got; got = mpz_kronecker_ui (a, b); if (got != answer) { printf ("mpz_kronecker_ui ("); mpz_out_str (stdout, 10, a); printf (", %lu) is %d should be %d\n", b, got, answer); abort (); } } void try_zi_si (mpz_srcptr a, long b, int answer) { int got; got = mpz_kronecker_si (a, b); if (got != answer) { printf ("mpz_kronecker_si ("); mpz_out_str (stdout, 10, a); printf (", %ld) is %d should be %d\n", b, got, answer); abort (); } } void try_ui_zi (unsigned long a, mpz_srcptr b, int answer) { int got; got = mpz_ui_kronecker (a, b); if (got != answer) { printf ("mpz_ui_kronecker (%lu, ", a); mpz_out_str (stdout, 10, b); printf (") is %d should be %d\n", got, answer); abort (); } } void try_si_zi (long a, mpz_srcptr b, int answer) { int got; got = mpz_si_kronecker (a, b); if (got != answer) { printf ("mpz_si_kronecker (%ld, ", a); mpz_out_str (stdout, 10, b); printf (") is %d should be %d\n", got, answer); abort (); } } /* Don't bother checking mpz_jacobi, since it only differs for b even, and we don't have an actual expected answer for it. tests/devel/try.c does some checks though. */ void try_zi_zi (mpz_srcptr a, mpz_srcptr b, int answer) { int got; got = mpz_kronecker (a, b); if (got != answer) { printf ("mpz_kronecker ("); mpz_out_str (stdout, 10, a); printf (", "); mpz_out_str (stdout, 10, b); printf (") is %d should be %d\n", got, answer); abort (); } } void try_pari (mpz_srcptr a, mpz_srcptr b, int answer) { printf ("try("); mpz_out_str (stdout, 10, a); printf (","); mpz_out_str (stdout, 10, b); printf (",%d)\n", answer); } void try_each (mpz_srcptr a, mpz_srcptr b, int answer) { #if 0 fprintf(stderr, "asize = %d, bsize = %d\n", mpz_sizeinbase (a, 2), mpz_sizeinbase (b, 2)); #endif if (option_pari) { try_pari (a, b, answer); return; } if (mpz_fits_ulimb_p (a) && mpz_fits_ulimb_p (b)) try_base (mpz_get_ulimb (a), mpz_get_ulimb (b), answer); if (mpz_fits_ulong_p (b)) try_zi_ui (a, mpz_get_ui (b), answer); if (mpz_fits_slong_p (b)) try_zi_si (a, mpz_get_si (b), answer); if (mpz_fits_ulong_p (a)) try_ui_zi (mpz_get_ui (a), b, answer); if (mpz_fits_sint_p (a)) try_si_zi (mpz_get_si (a), b, answer); try_zi_zi (a, b, answer); } /* Try (a/b) and (a/-b). */ void try_pn (mpz_srcptr a, mpz_srcptr b_orig, int answer) { mpz_t b; mpz_init_set (b, b_orig); try_each (a, b, answer); mpz_neg (b, b); if (mpz_sgn (a) < 0) answer = -answer; try_each (a, b, answer); mpz_clear (b); } /* Try (a+k*p/b) for various k, using the fact (a/b) is periodic in a with period p. For b>0, p=b if b!=2mod4 or p=4*b if b==2mod4. */ void try_periodic_num (mpz_srcptr a_orig, mpz_srcptr b, int answer) { mpz_t a, a_period; int i; if (mpz_sgn (b) <= 0) return; mpz_init_set (a, a_orig); mpz_init_set (a_period, b); if (mpz_mod4 (b) == 2) mpz_mul_ui (a_period, a_period, 4); /* don't bother with these tests if they're only going to produce even/even */ if (mpz_even_p (a) && mpz_even_p (b) && mpz_even_p (a_period)) goto done; for (i = 0; i < 6; i++) { mpz_add (a, a, a_period); try_pn (a, b, answer); } mpz_set (a, a_orig); for (i = 0; i < 6; i++) { mpz_sub (a, a, a_period); try_pn (a, b, answer); } done: mpz_clear (a); mpz_clear (a_period); } /* Try (a/b+k*p) for various k, using the fact (a/b) is periodic in b of period p. period p a==0,1mod4 a a==2mod4 4*a a==3mod4 and b odd 4*a a==3mod4 and b even 8*a In Henri Cohen's book the period is given as 4*a for all a==2,3mod4, but a counterexample would seem to be (3/2)=-1 which with (3/14)=+1 doesn't have period 4*a (but rather 8*a with (3/26)=-1). Maybe the plain 4*a is to be read as applying to a plain Jacobi symbol with b odd, rather than the Kronecker extension to b even. */ void try_periodic_den (mpz_srcptr a, mpz_srcptr b_orig, int answer) { mpz_t b, b_period; int i; if (mpz_sgn (a) == 0 || mpz_sgn (b_orig) == 0) return; mpz_init_set (b, b_orig); mpz_init_set (b_period, a); if (mpz_mod4 (a) == 3 && mpz_even_p (b)) mpz_mul_ui (b_period, b_period, 8L); else if (mpz_mod4 (a) >= 2) mpz_mul_ui (b_period, b_period, 4L); /* don't bother with these tests if they're only going to produce even/even */ if (mpz_even_p (a) && mpz_even_p (b) && mpz_even_p (b_period)) goto done; for (i = 0; i < 6; i++) { mpz_add (b, b, b_period); try_pn (a, b, answer); } mpz_set (b, b_orig); for (i = 0; i < 6; i++) { mpz_sub (b, b, b_period); try_pn (a, b, answer); } done: mpz_clear (b); mpz_clear (b_period); } static const unsigned long ktable[] = { 0, 1, 2, 3, 4, 5, 6, 7, GMP_NUMB_BITS-1, GMP_NUMB_BITS, GMP_NUMB_BITS+1, 2*GMP_NUMB_BITS-1, 2*GMP_NUMB_BITS, 2*GMP_NUMB_BITS+1, 3*GMP_NUMB_BITS-1, 3*GMP_NUMB_BITS, 3*GMP_NUMB_BITS+1 }; /* Try (a/b*2^k) for various k. */ void try_2den (mpz_srcptr a, mpz_srcptr b_orig, int answer) { mpz_t b; int kindex; int answer_a2, answer_k; unsigned long k; /* don't bother when b==0 */ if (mpz_sgn (b_orig) == 0) return; mpz_init_set (b, b_orig); /* (a/2) is 0 if a even, 1 if a==1 or 7 mod 8, -1 if a==3 or 5 mod 8 */ answer_a2 = (mpz_even_p (a) ? 0 : (((SIZ(a) >= 0 ? PTR(a)[0] : -PTR(a)[0]) + 2) & 7) < 4 ? 1 : -1); for (kindex = 0; kindex < numberof (ktable); kindex++) { k = ktable[kindex]; /* answer_k = answer*(answer_a2^k) */ answer_k = (answer_a2 == 0 && k != 0 ? 0 : (k & 1) == 1 && answer_a2 == -1 ? -answer : answer); mpz_mul_2exp (b, b_orig, k); try_pn (a, b, answer_k); } mpz_clear (b); } /* Try (a*2^k/b) for various k. If it happens mpz_ui_kronecker() gets (2/b) wrong it will show up as wrong answers demanded. */ void try_2num (mpz_srcptr a_orig, mpz_srcptr b, int answer) { mpz_t a; int kindex; int answer_2b, answer_k; unsigned long k; /* don't bother when a==0 */ if (mpz_sgn (a_orig) == 0) return; mpz_init (a); /* (2/b) is 0 if b even, 1 if b==1 or 7 mod 8, -1 if b==3 or 5 mod 8 */ answer_2b = (mpz_even_p (b) ? 0 : (((SIZ(b) >= 0 ? PTR(b)[0] : -PTR(b)[0]) + 2) & 7) < 4 ? 1 : -1); for (kindex = 0; kindex < numberof (ktable); kindex++) { k = ktable[kindex]; /* answer_k = answer*(answer_2b^k) */ answer_k = (answer_2b == 0 && k != 0 ? 0 : (k & 1) == 1 && answer_2b == -1 ? -answer : answer); mpz_mul_2exp (a, a_orig, k); try_pn (a, b, answer_k); } mpz_clear (a); } /* The try_2num() and try_2den() routines don't in turn call try_periodic_num() and try_periodic_den() because it hugely increases the number of tests performed, without obviously increasing coverage. Useful extra derived cases can be added here. */ void try_all (mpz_t a, mpz_t b, int answer) { try_pn (a, b, answer); try_periodic_num (a, b, answer); try_periodic_den (a, b, answer); try_2num (a, b, answer); try_2den (a, b, answer); } void check_data (void) { static const struct { const char *a; const char *b; int answer; } data[] = { /* Note that the various derived checks in try_all() reduce the cases that need to be given here. */ /* some zeros */ { "0", "0", 0 }, { "0", "2", 0 }, { "0", "6", 0 }, { "5", "0", 0 }, { "24", "60", 0 }, /* (a/1) = 1, any a In particular note (0/1)=1 so that (a/b)=(a mod b/b). */ { "0", "1", 1 }, { "1", "1", 1 }, { "2", "1", 1 }, { "3", "1", 1 }, { "4", "1", 1 }, { "5", "1", 1 }, /* (0/b) = 0, b != 1 */ { "0", "3", 0 }, { "0", "5", 0 }, { "0", "7", 0 }, { "0", "9", 0 }, { "0", "11", 0 }, { "0", "13", 0 }, { "0", "15", 0 }, /* (1/b) = 1 */ { "1", "1", 1 }, { "1", "3", 1 }, { "1", "5", 1 }, { "1", "7", 1 }, { "1", "9", 1 }, { "1", "11", 1 }, /* (-1/b) = (-1)^((b-1)/2) which is -1 for b==3 mod 4 */ { "-1", "1", 1 }, { "-1", "3", -1 }, { "-1", "5", 1 }, { "-1", "7", -1 }, { "-1", "9", 1 }, { "-1", "11", -1 }, { "-1", "13", 1 }, { "-1", "15", -1 }, { "-1", "17", 1 }, { "-1", "19", -1 }, /* (2/b) = (-1)^((b^2-1)/8) which is -1 for b==3,5 mod 8. try_2num() will exercise multiple powers of 2 in the numerator. */ { "2", "1", 1 }, { "2", "3", -1 }, { "2", "5", -1 }, { "2", "7", 1 }, { "2", "9", 1 }, { "2", "11", -1 }, { "2", "13", -1 }, { "2", "15", 1 }, { "2", "17", 1 }, /* (-2/b) = (-1)^((b^2-1)/8)*(-1)^((b-1)/2) which is -1 for b==5,7mod8. try_2num() will exercise multiple powers of 2 in the numerator, which will test that the shift in mpz_si_kronecker() uses unsigned not signed. */ { "-2", "1", 1 }, { "-2", "3", 1 }, { "-2", "5", -1 }, { "-2", "7", -1 }, { "-2", "9", 1 }, { "-2", "11", 1 }, { "-2", "13", -1 }, { "-2", "15", -1 }, { "-2", "17", 1 }, /* (a/2)=(2/a). try_2den() will exercise multiple powers of 2 in the denominator. */ { "3", "2", -1 }, { "5", "2", -1 }, { "7", "2", 1 }, { "9", "2", 1 }, { "11", "2", -1 }, /* Harriet Griffin, "Elementary Theory of Numbers", page 155, various examples. */ { "2", "135", 1 }, { "135", "19", -1 }, { "2", "19", -1 }, { "19", "135", 1 }, { "173", "135", 1 }, { "38", "135", 1 }, { "135", "173", 1 }, { "173", "5", -1 }, { "3", "5", -1 }, { "5", "173", -1 }, { "173", "3", -1 }, { "2", "3", -1 }, { "3", "173", -1 }, { "253", "21", 1 }, { "1", "21", 1 }, { "21", "253", 1 }, { "21", "11", -1 }, { "-1", "11", -1 }, /* Griffin page 147 */ { "-1", "17", 1 }, { "2", "17", 1 }, { "-2", "17", 1 }, { "-1", "89", 1 }, { "2", "89", 1 }, /* Griffin page 148 */ { "89", "11", 1 }, { "1", "11", 1 }, { "89", "3", -1 }, { "2", "3", -1 }, { "3", "89", -1 }, { "11", "89", 1 }, { "33", "89", -1 }, /* H. Davenport, "The Higher Arithmetic", page 65, the quadratic residues and non-residues mod 19. */ { "1", "19", 1 }, { "4", "19", 1 }, { "5", "19", 1 }, { "6", "19", 1 }, { "7", "19", 1 }, { "9", "19", 1 }, { "11", "19", 1 }, { "16", "19", 1 }, { "17", "19", 1 }, { "2", "19", -1 }, { "3", "19", -1 }, { "8", "19", -1 }, { "10", "19", -1 }, { "12", "19", -1 }, { "13", "19", -1 }, { "14", "19", -1 }, { "15", "19", -1 }, { "18", "19", -1 }, /* Residues and non-residues mod 13 */ { "0", "13", 0 }, { "1", "13", 1 }, { "2", "13", -1 }, { "3", "13", 1 }, { "4", "13", 1 }, { "5", "13", -1 }, { "6", "13", -1 }, { "7", "13", -1 }, { "8", "13", -1 }, { "9", "13", 1 }, { "10", "13", 1 }, { "11", "13", -1 }, { "12", "13", 1 }, /* various */ { "5", "7", -1 }, { "15", "17", 1 }, { "67", "89", 1 }, /* special values inducing a==b==1 at the end of jac_or_kron() */ { "0x10000000000000000000000000000000000000000000000001", "0x10000000000000000000000000000000000000000000000003", 1 }, /* Test for previous bugs in jacobi_2. */ { "0x43900000000", "0x42400000439", -1 }, /* 32-bit limbs */ { "0x4390000000000000000", "0x4240000000000000439", -1 }, /* 64-bit limbs */ { "198158408161039063", "198158360916398807", -1 }, /* Some tests involving large quotients in the continued fraction expansion. */ { "37200210845139167613356125645445281805", "451716845976689892447895811408978421929", -1 }, { "67674091930576781943923596701346271058970643542491743605048620644676477275152701774960868941561652032482173612421015", "4902678867794567120224500687210807069172039735", 0 }, { "2666617146103764067061017961903284334497474492754652499788571378062969111250584288683585223600172138551198546085281683283672592", "2666617146103764067061017961903284334497474492754652499788571378062969111250584288683585223600172138551198546085281683290481773", 1 }, /* Exersizes the case asize == 1, btwos > 0 in mpz_jacobi. */ { "804609", "421248363205206617296534688032638102314410556521742428832362659824", 1 } , { "4190209", "2239744742177804210557442048984321017460028974602978995388383905961079286530650825925074203175536427000", 1 }, /* Exersizes the case asize == 1, btwos = 63 in mpz_jacobi (relevant when GMP_LIMB_BITS == 64). */ { "17311973299000934401", "1675975991242824637446753124775689449936871337036614677577044717424700351103148799107651171694863695242089956242888229458836426332300124417011114380886016", 1 }, { "3220569220116583677", "41859917623035396746", -1 }, /* Other test cases that triggered bugs during development. */ { "37200210845139167613356125645445281805", "340116213441272389607827434472642576514", -1 }, { "74400421690278335226712251290890563610", "451716845976689892447895811408978421929", -1 }, }; int i; mpz_t a, b; mpz_init (a); mpz_init (b); for (i = 0; i < numberof (data); i++) { mpz_set_str_or_abort (a, data[i].a, 0); mpz_set_str_or_abort (b, data[i].b, 0); try_all (a, b, data[i].answer); } mpz_clear (a); mpz_clear (b); } /* (a^2/b)=1 if gcd(a,b)=1, or (a^2/b)=0 if gcd(a,b)!=1. This includes when a=0 or b=0. */ void check_squares_zi (void) { gmp_randstate_ptr rands = RANDS; mpz_t a, b, g; int i, answer; mp_size_t size_range, an, bn; mpz_t bs; mpz_init (bs); mpz_init (a); mpz_init (b); mpz_init (g); for (i = 0; i < 50; i++) { mpz_urandomb (bs, rands, 32); size_range = mpz_get_ui (bs) % 10 + i/8 + 2; mpz_urandomb (bs, rands, size_range); an = mpz_get_ui (bs); mpz_rrandomb (a, rands, an); mpz_urandomb (bs, rands, size_range); bn = mpz_get_ui (bs); mpz_rrandomb (b, rands, bn); mpz_gcd (g, a, b); if (mpz_cmp_ui (g, 1L) == 0) answer = 1; else answer = 0; mpz_mul (a, a, a); try_all (a, b, answer); } mpz_clear (bs); mpz_clear (a); mpz_clear (b); mpz_clear (g); } /* Check the handling of asize==0, make sure it isn't affected by the low limb. */ void check_a_zero (void) { mpz_t a, b; mpz_init_set_ui (a, 0); mpz_init (b); mpz_set_ui (b, 1L); PTR(a)[0] = 0; try_all (a, b, 1); /* (0/1)=1 */ PTR(a)[0] = 1; try_all (a, b, 1); /* (0/1)=1 */ mpz_set_si (b, -1L); PTR(a)[0] = 0; try_all (a, b, 1); /* (0/-1)=1 */ PTR(a)[0] = 1; try_all (a, b, 1); /* (0/-1)=1 */ mpz_set_ui (b, 0); PTR(a)[0] = 0; try_all (a, b, 0); /* (0/0)=0 */ PTR(a)[0] = 1; try_all (a, b, 0); /* (0/0)=0 */ mpz_set_ui (b, 2); PTR(a)[0] = 0; try_all (a, b, 0); /* (0/2)=0 */ PTR(a)[0] = 1; try_all (a, b, 0); /* (0/2)=0 */ mpz_clear (a); mpz_clear (b); } /* Assumes that b = prod p_k^e_k */ int ref_jacobi (mpz_srcptr a, mpz_srcptr b, unsigned nprime, mpz_t prime[], unsigned *exp) { unsigned i; int res; for (i = 0, res = 1; i < nprime; i++) if (exp[i]) { int legendre = refmpz_legendre (a, prime[i]); if (!legendre) return 0; if (exp[i] & 1) res *= legendre; } return res; } void check_jacobi_factored (void) { #define PRIME_N 10 #define PRIME_MAX_SIZE 50 #define PRIME_MAX_EXP 4 #define PRIME_A_COUNT 10 #define PRIME_B_COUNT 5 #define PRIME_MAX_B_SIZE 2000 gmp_randstate_ptr rands = RANDS; mpz_t prime[PRIME_N]; unsigned exp[PRIME_N]; mpz_t a, b, t, bs; unsigned i; mpz_init (a); mpz_init (b); mpz_init (t); mpz_init (bs); /* Generate primes */ for (i = 0; i < PRIME_N; i++) { mp_size_t size; mpz_init (prime[i]); mpz_urandomb (bs, rands, 32); size = mpz_get_ui (bs) % PRIME_MAX_SIZE + 2; mpz_rrandomb (prime[i], rands, size); if (mpz_cmp_ui (prime[i], 3) <= 0) mpz_set_ui (prime[i], 3); else mpz_nextprime (prime[i], prime[i]); } for (i = 0; i < PRIME_B_COUNT; i++) { unsigned j, k; mp_bitcnt_t bsize; mpz_set_ui (b, 1); bsize = 1; for (j = 0; j < PRIME_N && bsize < PRIME_MAX_B_SIZE; j++) { mpz_urandomb (bs, rands, 32); exp[j] = mpz_get_ui (bs) % PRIME_MAX_EXP; mpz_pow_ui (t, prime[j], exp[j]); mpz_mul (b, b, t); bsize = mpz_sizeinbase (b, 2); } for (k = 0; k < PRIME_A_COUNT; k++) { int answer; mpz_rrandomb (a, rands, bsize + 2); answer = ref_jacobi (a, b, j, prime, exp); try_all (a, b, answer); } } for (i = 0; i < PRIME_N; i++) mpz_clear (prime[i]); mpz_clear (a); mpz_clear (b); mpz_clear (t); mpz_clear (bs); #undef PRIME_N #undef PRIME_MAX_SIZE #undef PRIME_MAX_EXP #undef PRIME_A_COUNT #undef PRIME_B_COUNT #undef PRIME_MAX_B_SIZE } /* These tests compute (a|n), where the quotient sequence includes large quotients, and n has a known factorization. Such inputs are generated as follows. First, construct a large n, as a power of a prime p of moderate size. Next, compute a matrix from factors (q,1;1,0), with q chosen with uniformly distributed size. We must stop with matrix elements of roughly half the size of n. Denote elements of M as M = (m00, m01; m10, m11). We now look for solutions to n = m00 x + m01 y a = m10 x + m11 y with x,y > 0. Since n >= m00 * m01, there exists a positive solution to the first equation. Find those x, y, and substitute in the second equation to get a. Then the quotient sequence for (a|n) is precisely the quotients used when constructing M, followed by the quotient sequence for (x|y). Numbers should also be large enough that we exercise hgcd_jacobi, which means that they should be larger than max (GCD_DC_THRESHOLD, 3 * HGCD_THRESHOLD) With an n of roughly 40000 bits, this should hold on most machines. */ void check_large_quotients (void) { #define COUNT 50 #define PBITS 200 #define PPOWER 201 #define MAX_QBITS 500 gmp_randstate_ptr rands = RANDS; mpz_t p, n, q, g, s, t, x, y, bs; mpz_t M[2][2]; mp_bitcnt_t nsize; unsigned i; mpz_init (p); mpz_init (n); mpz_init (q); mpz_init (g); mpz_init (s); mpz_init (t); mpz_init (x); mpz_init (y); mpz_init (bs); mpz_init (M[0][0]); mpz_init (M[0][1]); mpz_init (M[1][0]); mpz_init (M[1][1]); /* First generate a number with known factorization, as a random smallish prime raised to an odd power. Then (a|n) = (a|p). */ mpz_rrandomb (p, rands, PBITS); mpz_nextprime (p, p); mpz_pow_ui (n, p, PPOWER); nsize = mpz_sizeinbase (n, 2); for (i = 0; i < COUNT; i++) { unsigned j; unsigned chain_len; int answer; mp_bitcnt_t msize; mpz_set_ui (M[0][0], 1); mpz_set_ui (M[0][1], 0); mpz_set_ui (M[1][0], 0); mpz_set_ui (M[1][1], 1); for (msize = 1; 2*(msize + MAX_QBITS) + 1 < nsize ;) { unsigned i; mpz_rrandomb (bs, rands, 32); mpz_rrandomb (q, rands, 1 + mpz_get_ui (bs) % MAX_QBITS); /* Multiply by (q, 1; 1,0) from the right */ for (i = 0; i < 2; i++) { mp_bitcnt_t size; mpz_swap (M[i][0], M[i][1]); mpz_addmul (M[i][0], M[i][1], q); size = mpz_sizeinbase (M[i][0], 2); if (size > msize) msize = size; } } mpz_gcdext (g, s, t, M[0][0], M[0][1]); ASSERT_ALWAYS (mpz_cmp_ui (g, 1) == 0); /* Solve n = M[0][0] * x + M[0][1] * y */ if (mpz_sgn (s) > 0) { mpz_mul (x, n, s); mpz_fdiv_qr (q, x, x, M[0][1]); mpz_mul (y, q, M[0][0]); mpz_addmul (y, t, n); ASSERT_ALWAYS (mpz_sgn (y) > 0); } else { mpz_mul (y, n, t); mpz_fdiv_qr (q, y, y, M[0][0]); mpz_mul (x, q, M[0][1]); mpz_addmul (x, s, n); ASSERT_ALWAYS (mpz_sgn (x) > 0); } mpz_mul (x, x, M[1][0]); mpz_addmul (x, y, M[1][1]); /* Now (x|n) has the selected large quotients */ answer = refmpz_legendre (x, p); try_zi_zi (x, n, answer); } mpz_clear (p); mpz_clear (n); mpz_clear (q); mpz_clear (g); mpz_clear (s); mpz_clear (t); mpz_clear (x); mpz_clear (y); mpz_clear (bs); mpz_clear (M[0][0]); mpz_clear (M[0][1]); mpz_clear (M[1][0]); mpz_clear (M[1][1]); #undef COUNT #undef PBITS #undef PPOWER #undef MAX_QBITS } int main (int argc, char *argv[]) { tests_start (); if (argc >= 2 && strcmp (argv[1], "-p") == 0) { option_pari = 1; printf ("\ try(a,b,answer) =\n\ {\n\ if (kronecker(a,b) != answer,\n\ print(\"wrong at \", a, \",\", b,\n\ \" expected \", answer,\n\ \" pari says \", kronecker(a,b)))\n\ }\n"); } check_data (); check_squares_zi (); check_a_zero (); check_jacobi_factored (); check_large_quotients (); tests_end (); exit (0); }