xref: /dflybsd-src/contrib/gmp/mpz/kronzu.c (revision 86d7f5d305c6adaa56ff4582ece9859d73106103)
186d7f5d3SJohn Marino /* mpz_kronecker_ui -- mpz+ulong 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 int
mpz_kronecker_ui(mpz_srcptr a,unsigned long b)2686d7f5d3SJohn Marino mpz_kronecker_ui (mpz_srcptr a, unsigned long b)
2786d7f5d3SJohn Marino {
2886d7f5d3SJohn Marino   mp_srcptr  a_ptr;
2986d7f5d3SJohn Marino   mp_size_t  a_size;
3086d7f5d3SJohn Marino   mp_limb_t  a_rem;
3186d7f5d3SJohn Marino   int        result_bit1;
3286d7f5d3SJohn Marino 
3386d7f5d3SJohn Marino   a_size = SIZ(a);
3486d7f5d3SJohn Marino   if (a_size == 0)
3586d7f5d3SJohn Marino     return JACOBI_0U (b);
3686d7f5d3SJohn Marino 
3786d7f5d3SJohn Marino   if (b > GMP_NUMB_MAX)
3886d7f5d3SJohn Marino     {
3986d7f5d3SJohn Marino       mp_limb_t  blimbs[2];
4086d7f5d3SJohn Marino       mpz_t      bz;
4186d7f5d3SJohn Marino       ALLOC(bz) = numberof (blimbs);
4286d7f5d3SJohn Marino       PTR(bz) = blimbs;
4386d7f5d3SJohn Marino       mpz_set_ui (bz, b);
4486d7f5d3SJohn Marino       return mpz_kronecker (a, bz);
4586d7f5d3SJohn Marino     }
4686d7f5d3SJohn Marino 
4786d7f5d3SJohn Marino   a_ptr = PTR(a);
4886d7f5d3SJohn Marino   if ((b & 1) != 0)
4986d7f5d3SJohn Marino     {
5086d7f5d3SJohn Marino       result_bit1 = JACOBI_ASGN_SU_BIT1 (a_size, b);
5186d7f5d3SJohn Marino     }
5286d7f5d3SJohn Marino   else
5386d7f5d3SJohn Marino     {
5486d7f5d3SJohn Marino       mp_limb_t  a_low = a_ptr[0];
5586d7f5d3SJohn Marino       int        twos;
5686d7f5d3SJohn Marino 
5786d7f5d3SJohn Marino       if (b == 0)
5886d7f5d3SJohn Marino         return JACOBI_LS0 (a_low, a_size);   /* (a/0) */
5986d7f5d3SJohn Marino 
6086d7f5d3SJohn Marino       if (! (a_low & 1))
6186d7f5d3SJohn Marino         return 0;  /* (even/even)=0 */
6286d7f5d3SJohn Marino 
6386d7f5d3SJohn Marino       /* (a/2)=(2/a) for a odd */
6486d7f5d3SJohn Marino       count_trailing_zeros (twos, b);
6586d7f5d3SJohn Marino       b >>= twos;
6686d7f5d3SJohn Marino       result_bit1 = (JACOBI_TWOS_U_BIT1 (twos, a_low)
6786d7f5d3SJohn Marino                      ^ JACOBI_ASGN_SU_BIT1 (a_size, b));
6886d7f5d3SJohn Marino     }
6986d7f5d3SJohn Marino 
7086d7f5d3SJohn Marino   if (b == 1)
7186d7f5d3SJohn Marino     return JACOBI_BIT1_TO_PN (result_bit1);  /* (a/1)=1 for any a */
7286d7f5d3SJohn Marino 
7386d7f5d3SJohn Marino   a_size = ABS(a_size);
7486d7f5d3SJohn Marino 
7586d7f5d3SJohn Marino   /* (a/b) = (a mod b / b) */
7686d7f5d3SJohn Marino   JACOBI_MOD_OR_MODEXACT_1_ODD (result_bit1, a_rem, a_ptr, a_size, b);
7786d7f5d3SJohn Marino   return mpn_jacobi_base (a_rem, (mp_limb_t) b, result_bit1);
7886d7f5d3SJohn Marino }
79