=================================================================== RCS file: /home/cvs/OpenXM_contrib/gmp/mpz/Attic/kronzs.c,v retrieving revision 1.1.1.1 retrieving revision 1.1.1.2 diff -u -p -r1.1.1.1 -r1.1.1.2 --- OpenXM_contrib/gmp/mpz/Attic/kronzs.c 2000/09/09 14:12:54 1.1.1.1 +++ OpenXM_contrib/gmp/mpz/Attic/kronzs.c 2003/08/25 16:06:33 1.1.1.2 @@ -1,7 +1,6 @@ -/* mpz_kronecker_si -- Kronecker/Jacobi symbol. */ +/* mpz_kronecker_si -- mpz+long Kronecker/Jacobi symbol. -/* -Copyright (C) 1999, 2000 Free Software Foundation, Inc. +Copyright 1999, 2000, 2001, 2002 Free Software Foundation, Inc. This file is part of the GNU MP Library. @@ -18,57 +17,77 @@ 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. -*/ +MA 02111-1307, USA. */ #include "gmp.h" #include "gmp-impl.h" #include "longlong.h" -/* This function is expected to be often used with b odd, so there's a test - for this before invoking count_trailing_zeros(). +/* This implementation depends on BITS_PER_MP_LIMB being even, so that + (a/2)^BITS_PER_MP_LIMB = 1 and so there's no need to pay attention to how + many low zero limbs are stripped. */ +#if BITS_PER_MP_LIMB % 2 != 0 +Error, error, unsupported BITS_PER_MP_LIMB +#endif - After the absolute value of b is established it's treated as an unsigned + +/* After the absolute value of b is established it's treated as an unsigned long, because 0x80..00 doesn't fit in a signed long. */ int -#if __STDC__ mpz_kronecker_si (mpz_srcptr a, long b) -#else -mpz_kronecker_si (a, b) - mpz_srcptr a; - long b; -#endif { - int result_bit1; - int twos; + mp_srcptr a_ptr = PTR(a); + mp_size_t a_size; + mp_limb_t a_rem, b_limb; + int result_bit1; - if (b == 0) - return JACOBI_Z0 (a); + a_size = SIZ(a); + if (a_size == 0) + return JACOBI_0S (b); - result_bit1 = JACOBI_BSGN_ZS_BIT1(a, b); - b = ABS (b); +#if GMP_NUMB_BITS < BITS_PER_ULONG + if (b > GMP_NUMB_MAX || b < -GMP_NUMB_MAX) + { + mp_limb_t blimbs[2]; + mpz_t bz; + ALLOC(bz) = numberof (blimbs); + PTR(bz) = blimbs; + mpz_set_si (bz, b); + return mpz_kronecker (a, bz); + } +#endif - if (b == 1) - return JACOBI_BIT1_TO_PN (result_bit1); /* (a/1) = 1 for any a */ + result_bit1 = JACOBI_BSGN_SS_BIT1 (a_size, b); + b_limb = (unsigned long) ABS (b); - if (b & 1) - return mpn_jacobi_base (mpz_fdiv_ui (a, b), b, result_bit1); - - /* result 0 if both a,b even */ - if (mpz_even_p (a)) - return 0; + if ((b_limb & 1) == 0) + { + mp_limb_t a_low = a_ptr[0]; + int twos; - /* (a/2)=(2/a) when a odd */ - count_trailing_zeros (twos, b); - result_bit1 ^= JACOBI_TWOS_U_BIT1 (twos, PTR(a)[0]); + if (b_limb == 0) + return JACOBI_LS0 (a_low, a_size); /* (a/0) */ - b = ((unsigned long) b) >> twos; - if (b == 1) - return JACOBI_BIT1_TO_PN (result_bit1); - else - return mpn_jacobi_base (mpz_fdiv_ui (a, b), b, result_bit1); + if (! (a_low & 1)) + return 0; /* (even/even)=0 */ + + /* (a/2)=(2/a) for a odd */ + count_trailing_zeros (twos, b_limb); + b_limb >>= twos; + result_bit1 ^= JACOBI_TWOS_U_BIT1 (twos, a_low); + } + + if (b_limb == 1) + return JACOBI_BIT1_TO_PN (result_bit1); /* (a/1)=1 for any a */ + + result_bit1 ^= JACOBI_ASGN_SU_BIT1 (a_size, b_limb); + a_size = ABS(a_size); + + /* (a/b) = (a mod b / b) */ + JACOBI_MOD_OR_MODEXACT_1_ODD (result_bit1, a_rem, a_ptr, a_size, b_limb); + return mpn_jacobi_base (a_rem, b_limb, result_bit1); }