186d7f5d3SJohn Marino /* mpz_kronecker_si -- mpz+long Kronecker/Jacobi symbol.
286d7f5d3SJohn Marino
386d7f5d3SJohn Marino Copyright 1999, 2000, 2001, 2002 Free Software Foundation, Inc.
486d7f5d3SJohn Marino
586d7f5d3SJohn Marino This file is part of the GNU MP Library.
686d7f5d3SJohn Marino
786d7f5d3SJohn Marino The GNU MP Library is free software; you can redistribute it and/or modify
886d7f5d3SJohn Marino it under the terms of the GNU Lesser General Public License as published by
986d7f5d3SJohn Marino the Free Software Foundation; either version 3 of the License, or (at your
1086d7f5d3SJohn Marino option) any later version.
1186d7f5d3SJohn Marino
1286d7f5d3SJohn Marino The GNU MP Library is distributed in the hope that it will be useful, but
1386d7f5d3SJohn Marino WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
1486d7f5d3SJohn Marino or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
1586d7f5d3SJohn Marino License for more details.
1686d7f5d3SJohn Marino
1786d7f5d3SJohn Marino You should have received a copy of the GNU Lesser General Public License
1886d7f5d3SJohn Marino along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
1986d7f5d3SJohn Marino
2086d7f5d3SJohn Marino #include "gmp.h"
2186d7f5d3SJohn Marino #include "gmp-impl.h"
2286d7f5d3SJohn Marino #include "longlong.h"
2386d7f5d3SJohn Marino
2486d7f5d3SJohn Marino
2586d7f5d3SJohn Marino /* After the absolute value of b is established it's treated as an unsigned
2686d7f5d3SJohn Marino long, because 0x80..00 doesn't fit in a signed long. */
2786d7f5d3SJohn Marino
2886d7f5d3SJohn Marino int
mpz_kronecker_si(mpz_srcptr a,long b)2986d7f5d3SJohn Marino mpz_kronecker_si (mpz_srcptr a, long b)
3086d7f5d3SJohn Marino {
3186d7f5d3SJohn Marino mp_srcptr a_ptr;
3286d7f5d3SJohn Marino mp_size_t a_size;
3386d7f5d3SJohn Marino mp_limb_t a_rem, b_limb;
3486d7f5d3SJohn Marino int result_bit1;
3586d7f5d3SJohn Marino
3686d7f5d3SJohn Marino a_size = SIZ(a);
3786d7f5d3SJohn Marino if (a_size == 0)
3886d7f5d3SJohn Marino return JACOBI_0S (b);
3986d7f5d3SJohn Marino
4086d7f5d3SJohn Marino #if GMP_NUMB_BITS < BITS_PER_ULONG
4186d7f5d3SJohn Marino if (b > GMP_NUMB_MAX || b < -GMP_NUMB_MAX)
4286d7f5d3SJohn Marino {
4386d7f5d3SJohn Marino mp_limb_t blimbs[2];
4486d7f5d3SJohn Marino mpz_t bz;
4586d7f5d3SJohn Marino ALLOC(bz) = numberof (blimbs);
4686d7f5d3SJohn Marino PTR(bz) = blimbs;
4786d7f5d3SJohn Marino mpz_set_si (bz, b);
4886d7f5d3SJohn Marino return mpz_kronecker (a, bz);
4986d7f5d3SJohn Marino }
5086d7f5d3SJohn Marino #endif
5186d7f5d3SJohn Marino
5286d7f5d3SJohn Marino result_bit1 = JACOBI_BSGN_SS_BIT1 (a_size, b);
5386d7f5d3SJohn Marino b_limb = (unsigned long) ABS (b);
5486d7f5d3SJohn Marino a_ptr = PTR(a);
5586d7f5d3SJohn Marino
5686d7f5d3SJohn Marino if ((b_limb & 1) == 0)
5786d7f5d3SJohn Marino {
5886d7f5d3SJohn Marino mp_limb_t a_low = a_ptr[0];
5986d7f5d3SJohn Marino int twos;
6086d7f5d3SJohn Marino
6186d7f5d3SJohn Marino if (b_limb == 0)
6286d7f5d3SJohn Marino return JACOBI_LS0 (a_low, a_size); /* (a/0) */
6386d7f5d3SJohn Marino
6486d7f5d3SJohn Marino if (! (a_low & 1))
6586d7f5d3SJohn Marino return 0; /* (even/even)=0 */
6686d7f5d3SJohn Marino
6786d7f5d3SJohn Marino /* (a/2)=(2/a) for a odd */
6886d7f5d3SJohn Marino count_trailing_zeros (twos, b_limb);
6986d7f5d3SJohn Marino b_limb >>= twos;
7086d7f5d3SJohn Marino result_bit1 ^= JACOBI_TWOS_U_BIT1 (twos, a_low);
7186d7f5d3SJohn Marino }
7286d7f5d3SJohn Marino
7386d7f5d3SJohn Marino if (b_limb == 1)
7486d7f5d3SJohn Marino return JACOBI_BIT1_TO_PN (result_bit1); /* (a/1)=1 for any a */
7586d7f5d3SJohn Marino
7686d7f5d3SJohn Marino result_bit1 ^= JACOBI_ASGN_SU_BIT1 (a_size, b_limb);
7786d7f5d3SJohn Marino a_size = ABS(a_size);
7886d7f5d3SJohn Marino
7986d7f5d3SJohn Marino /* (a/b) = (a mod b / b) */
8086d7f5d3SJohn Marino JACOBI_MOD_OR_MODEXACT_1_ODD (result_bit1, a_rem, a_ptr, a_size, b_limb);
8186d7f5d3SJohn Marino return mpn_jacobi_base (a_rem, b_limb, result_bit1);
8286d7f5d3SJohn Marino }
83