source/libs/gmp/gmp-src/mpz/jacobi.c

   1 /* mpz_jacobi, mpz_legendre, mpz_kronecker -- mpz/mpz Jacobi symbols.
   2
   3 Copyright 2000-2002, 2005, 2010-2012 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 <stdio.h>
  32 #include "gmp.h"
  33 #include "gmp-impl.h"
  34 #include "longlong.h"
  35
  36
  37 /* This code does triple duty as mpz_jacobi, mpz_legendre and
  38    mpz_kronecker. For ABI compatibility, the link symbol is
  39    __gmpz_jacobi, not __gmpz_kronecker, even though the latter would
  40    be more logical.
  41
  42    mpz_jacobi could assume b is odd, but the improvements from that seem
  43    small compared to other operations, and anything significant should be
  44    checked at run-time since we'd like odd b to go fast in mpz_kronecker
  45    too.
  46
  47    mpz_legendre could assume b is an odd prime, but knowing this doesn't
  48    present any obvious benefits.  Result 0 wouldn't arise (unless "a" is a
  49    multiple of b), but the checking for that takes little time compared to
  50    other operations.
  51
  52    Enhancements:
  53
  54    mpn_bdiv_qr should be used instead of mpn_tdiv_qr.
  55
  56 */
  57
  58 int
  59 mpz_jacobi (mpz_srcptr a, mpz_srcptr b)
  60 {
  61   mp_srcptr  asrcp, bsrcp;
  62   mp_size_t  asize, bsize;
  63   mp_limb_t  alow, blow;
  64   mp_ptr     ap, bp;
  65   unsigned   btwos;
  66   int        result_bit1;
  67   int        res;
  68   TMP_DECL;
  69
  70   asize = SIZ(a);
  71   asrcp = PTR(a);
  72   alow = asrcp[0];
  73
  74   bsize = SIZ(b);
  75   bsrcp = PTR(b);
  76   blow = bsrcp[0];
  77
  78   /* The MPN jacobi functions require positive a and b, and b odd. So
  79      we must to handle the cases of a or b zero, then signs, and then
  80      the case of even b.
  81   */
  82
  83   if (bsize == 0)
  84     /* (a/0) = [ a = 1 or a = -1 ] */
  85     return JACOBI_LS0 (alow, asize);
  86
  87   if (asize == 0)
  88     /* (0/b) = [ b = 1 or b = - 1 ] */
  89     return JACOBI_0LS (blow, bsize);
  90
  91   if ( (((alow | blow) & 1) == 0))
  92     /* Common factor of 2 ==> (a/b) = 0 */
  93     return 0;
  94
  95   if (bsize < 0)
  96     {
  97       /* (a/-1) = -1 if a < 0, +1 if a >= 0 */
  98       result_bit1 = (asize < 0) << 1;
  99       bsize = -bsize;
 100     }
 101   else
 102     result_bit1 = 0;
 103
 104   JACOBI_STRIP_LOW_ZEROS (result_bit1, alow, bsrcp, bsize, blow);
 105
 106   count_trailing_zeros (btwos, blow);
 107   blow >>= btwos;
 108
 109   if (bsize > 1 && btwos > 0)
 110     {
 111       mp_limb_t b1 = bsrcp[1];
 112       blow |= b1 << (GMP_NUMB_BITS - btwos);
 113       if (bsize == 2 && (b1 >> btwos) == 0)
 114         bsize = 1;
 115     }
 116
 117   if (asize < 0)
 118     {
 119       /* (-1/b) = -1 iff b = 3 (mod 4) */
 120       result_bit1 ^= JACOBI_N1B_BIT1(blow);
 121       asize = -asize;
 122     }
 123
 124   JACOBI_STRIP_LOW_ZEROS (result_bit1, blow, asrcp, asize, alow);
 125
 126   /* Ensure asize >= bsize. Take advantage of the generalized
 127      reciprocity law (a/b*2^n) = (b*2^n / a) * RECIP(a,b) */
 128
 129   if (asize < bsize)
 130     {
 131       MPN_SRCPTR_SWAP (asrcp, asize, bsrcp, bsize);
 132       MP_LIMB_T_SWAP (alow, blow);
 133
 134       /* NOTE: The value of alow (old blow) is a bit subtle. For this code
 135          path, we get alow as the low, always odd, limb of shifted A. Which is
 136          what we need for the reciprocity update below.
 137
 138          However, all other uses of alow assumes that it is *not*
 139          shifted. Luckily, alow matters only when either
 140
 141          + btwos > 0, in which case A is always odd
 142
 143          + asize == bsize == 1, in which case this code path is never
 144            taken. */
 145
 146       count_trailing_zeros (btwos, blow);
 147       blow >>= btwos;
 148
 149       if (bsize > 1 && btwos > 0)
 150         {
 151           mp_limb_t b1 = bsrcp[1];
 152           blow |= b1 << (GMP_NUMB_BITS - btwos);
 153           if (bsize == 2 && (b1 >> btwos) == 0)
 154             bsize = 1;
 155         }
 156
 157       result_bit1 ^= JACOBI_RECIP_UU_BIT1 (alow, blow);
 158     }
 159
 160   if (bsize == 1)
 161     {
 162       result_bit1 ^= JACOBI_TWOS_U_BIT1(btwos, alow);
 163
 164       if (blow == 1)
 165         return JACOBI_BIT1_TO_PN (result_bit1);
 166
 167       if (asize > 1)
 168         JACOBI_MOD_OR_MODEXACT_1_ODD (result_bit1, alow, asrcp, asize, blow);
 169
 170       return mpn_jacobi_base (alow, blow, result_bit1);
 171     }
 172
 173   /* Allocation strategy: For A, we allocate a working copy only for A % B, but
 174      when A is much larger than B, we have to allocate space for the large
 175      quotient. We use the same area, pointed to by bp, for both the quotient
 176      A/B and the working copy of B. */
 177
 178   TMP_MARK;
 179
 180   if (asize >= 2*bsize)
 181     TMP_ALLOC_LIMBS_2 (ap, bsize, bp, asize - bsize + 1);
 182   else
 183     TMP_ALLOC_LIMBS_2 (ap, bsize, bp, bsize);
 184
 185   /* In the case of even B, we conceptually shift out the powers of two first,
 186      and then divide A mod B. Hence, when taking those powers of two into
 187      account, we must use alow *before* the division. Doing the actual division
 188      first is ok, because the point is to remove multiples of B from A, and
 189      multiples of 2^k B are good enough. */
 190   if (asize > bsize)
 191     mpn_tdiv_qr (bp, ap, 0, asrcp, asize, bsrcp, bsize);
 192   else
 193     MPN_COPY (ap, asrcp, bsize);
 194
 195   if (btwos > 0)
 196     {
 197       result_bit1 ^= JACOBI_TWOS_U_BIT1(btwos, alow);
 198
 199       ASSERT_NOCARRY (mpn_rshift (bp, bsrcp, bsize, btwos));
 200       bsize -= (ap[bsize-1] | bp[bsize-1]) == 0;
 201     }
 202   else
 203     MPN_COPY (bp, bsrcp, bsize);
 204
 205   ASSERT (blow == bp[0]);
 206   res = mpn_jacobi_n (ap, bp, bsize,
 207                       mpn_jacobi_init (ap[0], blow, (result_bit1>>1) & 1));
 208
 209   TMP_FREE;
 210   return res;
 211 }