diff options
author | Kevin Ryde <user42@zip.com.au> | 2001-02-02 00:55:24 +0100 |
---|---|---|
committer | Kevin Ryde <user42@zip.com.au> | 2001-02-02 00:55:24 +0100 |
commit | 9c532b656dd173c189512ec2c2375b535542b933 (patch) | |
tree | a97bc3768c11152b166692f18fa17e7b1c33ab18 /tests/refmpz.c | |
parent | bd03a38085e0efc5466a8b1d98b835b3070ef3ae (diff) | |
download | gmp-9c532b656dd173c189512ec2c2375b535542b933.tar.gz |
* tests/refmpz.c: New file.
Diffstat (limited to 'tests/refmpz.c')
-rw-r--r-- | tests/refmpz.c | 199 |
1 files changed, 199 insertions, 0 deletions
diff --git a/tests/refmpz.c b/tests/refmpz.c new file mode 100644 index 000000000..5b9fc1367 --- /dev/null +++ b/tests/refmpz.c @@ -0,0 +1,199 @@ +/* Reference mpz functions */ + +/* +Copyright 2000, 2001 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 of the License, or (at your +option) any later version. + +The GNU MP Library 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 Lesser General Public +License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. +*/ + +/* always do assertion checking */ +#define WANT_ASSERT 1 + +#include <stdio.h> +#include "gmp.h" +#include "gmp-impl.h" +#include "longlong.h" +#include "tests.h" + + +/* Change this to "#define TRACE(x) x" for some traces. */ +#define TRACE(x) + + +/* (0/b), with mpz b; is 1 if b=+/-1, 0 otherwise */ +#define JACOBI_0Z(b) JACOBI_0LS (PTR(b)[0], SIZ(b)) + +/* (a/b) effect due to sign of b: mpz/mpz */ +#define JACOBI_BSGN_ZZ_BIT1(a, b) JACOBI_BSGN_SS_BIT1 (SIZ(a), SIZ(b)) + +/* (a/b) effect due to sign of a: mpz/unsigned-mpz, b odd; + is (-1/b) if a<0, or +1 if a>=0 */ +#define JACOBI_ASGN_ZZU_BIT1(a, b) JACOBI_ASGN_SU_BIT1 (SIZ(a), PTR(b)[0]) + + +int +refmpz_kronecker (mpz_srcptr a_orig, mpz_srcptr b_orig) +{ + unsigned long twos; + mpz_t a, b; + int result_bit1 = 0; + + if (mpz_sgn (b_orig) == 0) + return JACOBI_Z0 (a_orig); /* (a/0) */ + + if (mpz_sgn (a_orig) == 0) + return JACOBI_0Z (b_orig); /* (0/b) */ + + if (mpz_even_p (a_orig) && mpz_even_p (b_orig)) + return 0; + + if (mpz_cmp_ui (b_orig, 1) == 0) + return 1; + + mpz_init_set (a, a_orig); + mpz_init_set (b, b_orig); + + if (mpz_sgn (b) < 0) + { + result_bit1 ^= JACOBI_BSGN_ZZ_BIT1 (a, b); + mpz_neg (b, b); + } + if (mpz_even_p (b)) + { + twos = mpz_scan1 (b, 0L); + mpz_tdiv_q_2exp (b, b, twos); + result_bit1 ^= JACOBI_TWOS_U_BIT1 (twos, PTR(a)[0]); + } + + if (mpz_sgn (a) < 0) + { + result_bit1 ^= JACOBI_N1B_BIT1 (PTR(b)[0]); + mpz_neg (a, a); + } + if (mpz_even_p (a)) + { + twos = mpz_scan1 (a, 0L); + mpz_tdiv_q_2exp (a, a, twos); + result_bit1 ^= JACOBI_TWOS_U_BIT1 (twos, PTR(b)[0]); + } + + for (;;) + { + ASSERT (mpz_odd_p (a)); + ASSERT (mpz_odd_p (b)); + ASSERT (mpz_sgn (a) > 0); + ASSERT (mpz_sgn (b) > 0); + + TRACE (printf ("top\n"); + mpz_trace (" a", a); + mpz_trace (" b", b)); + + if (mpz_cmp (a, b) < 0) + { + TRACE (printf ("swap\n")); + mpz_swap (a, b); + result_bit1 ^= JACOBI_RECIP_UU_BIT1 (PTR(a)[0], PTR(b)[0]); + } + + if (mpz_cmp_ui (b, 1) == 0) + break; + + mpz_sub (a, a, b); + TRACE (printf ("sub\n"); + mpz_trace (" a", a)); + if (mpz_sgn (a) == 0) + goto zero; + + twos = mpz_scan1 (a, 0L); + mpz_fdiv_q_2exp (a, a, twos); + TRACE (printf ("twos %lu\n", twos); + mpz_trace (" a", a)); + result_bit1 ^= JACOBI_TWOS_U_BIT1 (twos, PTR(b)[0]); + } + + mpz_clear (a); + mpz_clear (b); + return JACOBI_BIT1_TO_PN (result_bit1); + + zero: + mpz_clear (a); + mpz_clear (b); + return 0; +} + +/* Same as mpz_kronecker, but ignoring factors of 2 on b */ +int +refmpz_jacobi (mpz_srcptr a, mpz_srcptr b) +{ + mpz_t b_odd; + mpz_init_set (b_odd, b); + if (mpz_sgn (b_odd) != 0) + mpz_fdiv_q_2exp (b_odd, b_odd, mpz_scan1 (b_odd, 0L)); + return refmpz_kronecker (a, b_odd); +} + +int +refmpz_legendre (mpz_srcptr a, mpz_srcptr b) +{ + return refmpz_jacobi (a, b); +} + + +int +refmpz_kronecker_ui (mpz_srcptr a, unsigned long b) +{ + mpz_t bz; + int ret; + mpz_init_set_ui (bz, b); + ret = refmpz_kronecker (a, bz); + mpz_clear (bz); + return ret; +} + +int +refmpz_kronecker_si (mpz_srcptr a, long b) +{ + mpz_t bz; + int ret; + mpz_init_set_si (bz, b); + ret = refmpz_kronecker (a, bz); + mpz_clear (bz); + return ret; +} + +int +refmpz_ui_kronecker (unsigned long a, mpz_srcptr b) +{ + mpz_t az; + int ret; + mpz_init_set_ui (az, a); + ret = refmpz_kronecker (az, b); + mpz_clear (az); + return ret; +} + +int +refmpz_si_kronecker (long a, mpz_srcptr b) +{ + mpz_t az; + int ret; + mpz_init_set_si (az, a); + ret = refmpz_kronecker (az, b); + mpz_clear (az); + return ret; +} |