Home | History | Annotate | Line # | Download | only in mpz 
      1 /* mpz_si_kronecker -- long+mpz Kronecker/Jacobi symbol.
      2 
      3 Copyright 1999-2002 Free Software Foundation, Inc.
      4 
      5 This file is part of the GNU MP Library.
      6 
      7 The GNU MP Library is free software; you can redistribute it and/or modify
      8 it under the terms of either:
      9 
     10   * the GNU Lesser General Public License as published by the Free
     11     Software Foundation; either version 3 of the License, or (at your
     12     option) any later version.
     13 
     14 or
     15 
     16   * the GNU General Public License as published by the Free Software
     17     Foundation; either version 2 of the License, or (at your option) any
     18     later version.
     19 
     20 or both in parallel, as here.
     21 
     22 The GNU MP Library is distributed in the hope that it will be useful, but
     23 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
     24 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
     25 for more details.
     26 
     27 You should have received copies of the GNU General Public License and the
     28 GNU Lesser General Public License along with the GNU MP Library.  If not,
     29 see https://www.gnu.org/licenses/.  */
     30 
     31 #include "gmp-impl.h"
     32 #include "longlong.h"
     33 
     34 
     35 int
     36 mpz_si_kronecker (long a, mpz_srcptr b)
     37 {
     38   mp_srcptr  b_ptr;
     39   mp_limb_t  b_low;
     40   mp_size_t  b_size;
     41   mp_size_t  b_abs_size;
     42   mp_limb_t  a_limb, b_rem;
     43   unsigned   twos;
     44   int        result_bit1;
     45 
     46 #if GMP_NUMB_BITS < BITS_PER_ULONG
     47   if (a > GMP_NUMB_MAX || a < -GMP_NUMB_MAX)
     48     {
     49       mp_limb_t  alimbs[2];
     50       mpz_t      az;
     51       ALLOC(az) = numberof (alimbs);
     52       PTR(az) = alimbs;
     53       mpz_set_si (az, a);
     54       return mpz_kronecker (az, b);
     55     }
     56 #endif
     57 
     58   b_size = SIZ (b);
     59   if (b_size == 0)
     60     return JACOBI_S0 (a);  /* (a/0) */
     61 
     62   /* account for the effect of the sign of b, then ignore it */
     63   result_bit1 = JACOBI_BSGN_SS_BIT1 (a, b_size);
     64 
     65   b_ptr = PTR(b);
     66   b_low = b_ptr[0];
     67   b_abs_size = ABS (b_size);
     68 
     69   if ((b_low & 1) != 0)
     70     {
     71       /* b odd */
     72 
     73       result_bit1 ^= JACOBI_ASGN_SU_BIT1 (a, b_low);
     74       a_limb = ABS_CAST(mp_limb_t, a);
     75 
     76       if ((a_limb & 1) == 0)
     77 	{
     78 	  /* (0/b)=1 for b=+/-1, 0 otherwise */
     79 	  if (a_limb == 0)
     80 	    return (b_abs_size == 1 && b_low == 1);
     81 
     82 	  /* a even, b odd */
     83 	  count_trailing_zeros (twos, a_limb);
     84 	  a_limb >>= twos;
     85 	  /* (a*2^n/b) = (a/b) * twos(n,a) */
     86 	  result_bit1 ^= JACOBI_TWOS_U_BIT1 (twos, b_low);
     87 	}
     88     }
     89   else
     90     {
     91       /* (even/even)=0, and (0/b)=0 for b!=+/-1 */
     92       if ((a & 1) == 0)
     93 	return 0;
     94 
     95       /* a odd, b even
     96 
     97 	 Establish shifted b_low with valid bit1 for ASGN and RECIP below.
     98 	 Zero limbs stripped are accounted for, but zero bits on b_low are
     99 	 not because they remain in {b_ptr,b_abs_size} for the
    100 	 JACOBI_MOD_OR_MODEXACT_1_ODD. */
    101 
    102       JACOBI_STRIP_LOW_ZEROS (result_bit1, a, b_ptr, b_abs_size, b_low);
    103       if ((b_low & 1) == 0)
    104 	{
    105 	  if (UNLIKELY (b_low == GMP_NUMB_HIGHBIT))
    106 	    {
    107 	      /* need b_ptr[1] to get bit1 in b_low */
    108 	      if (b_abs_size == 1)
    109 		{
    110 		  /* (a/0x80000000) = (a/2)^(BPML-1) */
    111 		  if ((GMP_NUMB_BITS % 2) == 0)
    112 		    result_bit1 ^= JACOBI_TWO_U_BIT1 (a);
    113 		  return JACOBI_BIT1_TO_PN (result_bit1);
    114 		}
    115 
    116 	      /* b_abs_size > 1 */
    117 	      b_low = b_ptr[1] << 1;
    118 	    }
    119 	  else
    120 	    {
    121 	      count_trailing_zeros (twos, b_low);
    122 	      b_low >>= twos;
    123 	    }
    124 	}
    125 
    126       result_bit1 ^= JACOBI_ASGN_SU_BIT1 (a, b_low);
    127       a_limb = (unsigned long) ABS(a);
    128     }
    129 
    130   if (a_limb == 1)
    131     return JACOBI_BIT1_TO_PN (result_bit1);  /* (1/b)=1 */
    132 
    133   /* (a/b*2^n) = (b*2^n mod a / a) * recip(a,b) */
    134   JACOBI_MOD_OR_MODEXACT_1_ODD (result_bit1, b_rem, b_ptr, b_abs_size, a_limb);
    135   result_bit1 ^= JACOBI_RECIP_UU_BIT1 (a_limb, b_low);
    136   return mpn_jacobi_base (b_rem, a_limb, result_bit1);
    137 }
    138