contrib/gmp/mpz/kronuz.c

   1 /* mpz_ui_kronecker -- ulong+mpz Kronecker/Jacobi symbol.
   2
   3 Copyright 1999, 2000, 2001, 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 the GNU Lesser General Public License as published by
   9 the Free Software Foundation; either version 3 of the License, or (at your
  10 option) any later version.
  11
  12 The GNU MP Library is distributed in the hope that it will be useful, but
  13 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  14 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  15 License for more details.
  16
  17 You should have received a copy of the GNU Lesser General Public License
  18 along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  19
  20 #include "gmp.h"
  21 #include "gmp-impl.h"
  22 #include "longlong.h"
  23
  24
  25 int
  26 mpz_ui_kronecker (unsigned long a, mpz_srcptr b)
  27 {
  28   mp_srcptr  b_ptr;
  29   mp_limb_t  b_low;
  30   int        b_abs_size;
  31   mp_limb_t  b_rem;
  32   int        twos;
  33   int        result_bit1;
  34
  35   /* (a/-1)=1 when a>=0, so the sign of b is ignored */
  36   b_abs_size = ABSIZ (b);
  37
  38   if (b_abs_size == 0)
  39     return JACOBI_U0 (a);  /* (a/0) */
  40
  41   if (a > GMP_NUMB_MAX)
  42     {
  43       mp_limb_t  alimbs[2];
  44       mpz_t      az;
  45       ALLOC(az) = numberof (alimbs);
  46       PTR(az) = alimbs;
  47       mpz_set_ui (az, a);
  48       return mpz_kronecker (az, b);
  49     }
  50
  51   b_ptr = PTR(b);
  52   b_low = b_ptr[0];
  53   result_bit1 = 0;
  54
  55   if (! (b_low & 1))
  56     {
  57       /* (0/b)=0 for b!=+/-1; and (even/even)=0 */
  58       if (! (a & 1))
  59         return 0;
  60
  61       /* a odd, b even
  62
  63          Establish shifted b_low with valid bit1 for the RECIP below.  Zero
  64          limbs stripped are accounted for, but zero bits on b_low are not
  65          because they remain in {b_ptr,b_abs_size} for
  66          JACOBI_MOD_OR_MODEXACT_1_ODD. */
  67
  68       JACOBI_STRIP_LOW_ZEROS (result_bit1, a, b_ptr, b_abs_size, b_low);
  69       if (! (b_low & 1))
  70         {
  71           if (UNLIKELY (b_low == GMP_NUMB_HIGHBIT))
  72             {
  73               /* need b_ptr[1] to get bit1 in b_low */
  74               if (b_abs_size == 1)
  75                 {
  76                   /* (a/0x80...00) == (a/2)^(NUMB-1) */
  77                   if ((GMP_NUMB_BITS % 2) == 0)
  78                     {
  79                       /* JACOBI_STRIP_LOW_ZEROS does nothing to result_bit1
  80                          when GMP_NUMB_BITS is even, so it's still 0. */
  81                       ASSERT (result_bit1 == 0);
  82                       result_bit1 = JACOBI_TWO_U_BIT1 (a);
  83                     }
  84                   return JACOBI_BIT1_TO_PN (result_bit1);
  85                 }
  86
  87               /* b_abs_size > 1 */
  88               b_low = b_ptr[1] << 1;
  89             }
  90           else
  91             {
  92               count_trailing_zeros (twos, b_low);
  93               b_low >>= twos;
  94             }
  95         }
  96     }
  97   else
  98     {
  99       if (a == 0)        /* (0/b)=1 for b=+/-1, 0 otherwise */
 100         return (b_abs_size == 1 && b_low == 1);
 101
 102       if (! (a & 1))
 103         {
 104           /* a even, b odd */
 105           count_trailing_zeros (twos, a);
 106           a >>= twos;
 107           /* (a*2^n/b) = (a/b) * (2/a)^n */
 108           result_bit1 = JACOBI_TWOS_U_BIT1 (twos, b_low);
 109         }
 110     }
 111
 112   if (a == 1)
 113     return JACOBI_BIT1_TO_PN (result_bit1);  /* (1/b)=1 */
 114
 115   /* (a/b*2^n) = (b*2^n mod a / a) * RECIP(a,b) */
 116   JACOBI_MOD_OR_MODEXACT_1_ODD (result_bit1, b_rem, b_ptr, b_abs_size, a);
 117   result_bit1 ^= JACOBI_RECIP_UU_BIT1 (a, b_low);
 118   return mpn_jacobi_base (b_rem, (mp_limb_t) a, result_bit1);
 119 }