From 11da511c784eca003deb90c23570f0873954e0de Mon Sep 17 00:00:00 2001 From: Duncan Wilkie Date: Sat, 18 Nov 2023 06:11:09 -0600 Subject: Initial commit. --- gmp-6.3.0/tests/mpz/t-jac.c | 1008 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 1008 insertions(+) create mode 100644 gmp-6.3.0/tests/mpz/t-jac.c (limited to 'gmp-6.3.0/tests/mpz/t-jac.c') diff --git a/gmp-6.3.0/tests/mpz/t-jac.c b/gmp-6.3.0/tests/mpz/t-jac.c new file mode 100644 index 0000000..ed8aff8 --- /dev/null +++ b/gmp-6.3.0/tests/mpz/t-jac.c @@ -0,0 +1,1008 @@ +/* Exercise mpz_*_kronecker_*() and mpz_jacobi() functions. + +Copyright 1999-2004, 2013 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 https://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 + + Enhancements: + + More big test cases than those given by check_squares_zi would be good. */ + + +#include +#include +#include + +#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 }, + + /* Exercises the case asize == 1, btwos > 0 in mpz_jacobi. */ + { "804609", "421248363205206617296534688032638102314410556521742428832362659824", 1 } , + { "4190209", "2239744742177804210557442048984321017460028974602978995388383905961079286530650825925074203175536427000", 1 }, + + /* Exercises 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++) + { + 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); +} -- cgit v1.2.3