=================================================================== RCS file: /home/cvs/OpenXM_contrib/gmp/mpz/Attic/kronsz.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/kronsz.c 2000/09/09 14:12:54 1.1.1.1 +++ OpenXM_contrib/gmp/mpz/Attic/kronsz.c 2003/08/25 16:06:33 1.1.1.2 @@ -1,7 +1,6 @@ -/* mpz_si_kronecker -- Kronecker/Jacobi symbol. */ +/* mpz_si_kronecker -- long+mpz 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,109 +17,110 @@ 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 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 + + int -#if __STDC__ mpz_si_kronecker (long a, mpz_srcptr b) -#else -mpz_si_kronecker (a, b) - long a; - mpz_srcptr b; -#endif { - int b_abs_size; - mp_srcptr b_ptr; - mp_limb_t b_low; - int twos; + mp_srcptr b_ptr = PTR(b); + mp_limb_t b_low = b_ptr[0]; + mp_size_t b_size = SIZ (b); + mp_size_t b_abs_size = ABS (b_size); + mp_limb_t a_limb, b_rem; + unsigned twos; int result_bit1; - b_abs_size = ABSIZ (b); +#if GMP_NUMB_BITS < BITS_PER_ULONG + if (a > GMP_NUMB_MAX || a < -GMP_NUMB_MAX) + { + mp_limb_t alimbs[2]; + mpz_t az; + ALLOC(az) = numberof (alimbs); + PTR(az) = alimbs; + mpz_set_si (az, a); + return mpz_kronecker (az, b); + } +#endif + if (b_abs_size == 0) return JACOBI_S0 (a); /* (a/0) */ - b_ptr = PTR(b); - b_low = b_ptr[0]; + /* account for the effect of the sign of b, then ignore it */ + result_bit1 = JACOBI_BSGN_SS_BIT1 (a, b_size); - /* (0/b) = 1 if b=+/-1, 0 otherwise */ - if (a == 0) - return (b_abs_size == 1) & (b_low == 1); + if ((b_low & 1) != 0) + { + /* b odd */ - /* account for the effect of the sign of b, so can then ignore it */ - result_bit1 = JACOBI_BSGN_SZ_BIT1 (a, b); + result_bit1 ^= JACOBI_ASGN_SU_BIT1 (a, b_low); + a_limb = (unsigned long) ABS(a); - if ((b_low & 1) == 0) - { - /* b even */ + if ((a_limb & 1) == 0) + { + /* (0/b)=1 for b=+/-1, 0 otherwise */ + if (a_limb == 0) + return (b_abs_size == 1 && b_low == 1); + /* a even, b odd */ + count_trailing_zeros (twos, a_limb); + a_limb >>= twos; + /* (a*2^n/b) = (a/b) * twos(n,a) */ + result_bit1 ^= JACOBI_TWOS_U_BIT1 (twos, b_low); + } + } + else + { + /* (even/even)=0, and (0/b)=0 for b!=+/-1 */ if ((a & 1) == 0) - return 0; /* (a/b)=0 if both a,b even */ + return 0; - /* Require MP_BITS_PER_LIMB even, so that (a/2)^MP_BITS_PER_LIMB = 1, - and so that therefore there's no need to account for how many zero - limbs are stripped. */ - ASSERT ((BITS_PER_MP_LIMB & 1) == 0); + /* a odd, b even */ - MPN_STRIP_LOW_ZEROS_NOT_ZERO (b_ptr, b_abs_size); - b_low = b_ptr[0]; - + /* establish shifted b_low with valid bit1 for use with ASGN and RECIP */ + MPN_STRIP_LOW_ZEROS_NOT_ZERO (b_ptr, b_abs_size, b_low); if ((b_low & 1) == 0) { - /* odd a, even b */ + if (b_low == GMP_NUMB_HIGHBIT) + { + if (b_abs_size == 1) + { + /* (a/0x80000000) = (a/2)^(BPML-1) */ + result_bit1 ^= JACOBI_TWOS_U_BIT1 (GMP_NUMB_BITS - 1, a); + return JACOBI_BIT1_TO_PN (result_bit1); + } - mp_limb_t b_shl_bit1; - - count_trailing_zeros (twos, b_low); - - /* b_shl_bit1 is b>>twos, but with only bit 1 guaranteed */ - if (twos == BITS_PER_MP_LIMB-1) - b_shl_bit1 = (b_abs_size == 1) ? 0 : (b_ptr[1] << 1); + /* b_abs_size > 1 */ + b_low = b_ptr[1] << 1; + } else - b_shl_bit1 = (b_low >> twos); - - result_bit1 ^= JACOBI_ASGN_SU_BIT1 (a, b_shl_bit1); - a = ABS(a); - - if (a == 1) - return JACOBI_BIT1_TO_PN (result_bit1); /* (1/b)=1 */ - - /* twos (a/2), reciprocity to (b/a), and (b/a) = (b mod a / b) */ - return mpn_jacobi_base (mpn_mod_1_rshift (b_ptr, b_abs_size, - twos, a), - a, - result_bit1 - ^ JACOBI_TWOS_U_BIT1 (twos, a) - ^ JACOBI_RECIP_UU_BIT1 (a, b_shl_bit1)); + { + count_trailing_zeros (twos, b_low); + b_low >>= twos; + } } - } - /* b odd */ - - result_bit1 ^= JACOBI_ASGN_SU_BIT1 (a, b_low); - a = ABS(a); - - /* (a/1) = 1 for any a */ - if (b_abs_size == 1 && b_low == 1) - return JACOBI_BIT1_TO_PN (result_bit1); - - /* Note a is cast to unsigned because 0x80..00 doesn't fit in a signed. */ - if ((a & 1) == 0) - { - count_trailing_zeros (twos, a); - a = ((unsigned long) a) >> twos; - result_bit1 ^= JACOBI_TWOS_U_BIT1 (twos, b_low); + result_bit1 ^= JACOBI_ASGN_SU_BIT1 (a, b_low); + a_limb = (unsigned long) ABS(a); } - if (a == 1) + if (a_limb == 1) return JACOBI_BIT1_TO_PN (result_bit1); /* (1/b)=1 */ - /* reciprocity to (b/a), and (b/a) == (b mod a / a) */ - return mpn_jacobi_base (mpn_mod_1 (b_ptr, b_abs_size, a), a, - result_bit1 ^ JACOBI_RECIP_UU_BIT1 (a, b_low)); + /* (a/b*2^n) = (b*2^n mod a / a) * recip(a,b) */ + JACOBI_MOD_OR_MODEXACT_1_ODD (result_bit1, b_rem, b_ptr, b_abs_size, a_limb); + result_bit1 ^= JACOBI_RECIP_UU_BIT1 (a_limb, b_low); + return mpn_jacobi_base (b_rem, a_limb, result_bit1); }