source/libs/gmp/gmp-src/mpn/generic/jacbase.c

   1 /* mpn_jacobi_base -- limb/limb Jacobi symbol with restricted arguments.
   2
   3    THIS INTERFACE IS PRELIMINARY AND MIGHT DISAPPEAR OR BE SUBJECT TO
   4    INCOMPATIBLE CHANGES IN A FUTURE RELEASE OF GMP.
   5
   6 Copyright 1999-2002, 2010 Free Software Foundation, Inc.
   7
   8 This file is part of the GNU MP Library.
   9
  10 The GNU MP Library is free software; you can redistribute it and/or modify
  11 it under the terms of either:
  12
  13   * the GNU Lesser General Public License as published by the Free
  14     Software Foundation; either version 3 of the License, or (at your
  15     option) any later version.
  16
  17 or
  18
  19   * the GNU General Public License as published by the Free Software
  20     Foundation; either version 2 of the License, or (at your option) any
  21     later version.
  22
  23 or both in parallel, as here.
  24
  25 The GNU MP Library is distributed in the hope that it will be useful, but
  26 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  27 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  28 for more details.
  29
  30 You should have received copies of the GNU General Public License and the
  31 GNU Lesser General Public License along with the GNU MP Library.  If not,
  32 see https://www.gnu.org/licenses/.  */
  33
  34 #include "gmp.h"
  35 #include "gmp-impl.h"
  36 #include "longlong.h"
  37
  38
  39 /* Use the simple loop by default.  The generic count_trailing_zeros is not
  40    very fast, and the extra trickery of method 3 has proven to be less use
  41    than might have been though.  */
  42 #ifndef JACOBI_BASE_METHOD
  43 #define JACOBI_BASE_METHOD  2
  44 #endif
  45
  46
  47 /* Use count_trailing_zeros.  */
  48 #if JACOBI_BASE_METHOD == 1
  49 #define PROCESS_TWOS_ANY                                \
  50   {                                                     \
  51     mp_limb_t  twos;                                    \
  52     count_trailing_zeros (twos, a);                     \
  53     result_bit1 ^= JACOBI_TWOS_U_BIT1 (twos, b);        \
  54     a >>= twos;                                         \
  55   }
  56 #define PROCESS_TWOS_EVEN  PROCESS_TWOS_ANY
  57 #endif
  58
  59 /* Use a simple loop.  A disadvantage of this is that there's a branch on a
  60    50/50 chance of a 0 or 1 low bit.  */
  61 #if JACOBI_BASE_METHOD == 2
  62 #define PROCESS_TWOS_EVEN               \
  63   {                                     \
  64     int  two;                           \
  65     two = JACOBI_TWO_U_BIT1 (b);        \
  66     do                                  \
  67       {                                 \
  68         a >>= 1;                        \
  69         result_bit1 ^= two;             \
  70         ASSERT (a != 0);                \
  71       }                                 \
  72     while ((a & 1) == 0);               \
  73   }
  74 #define PROCESS_TWOS_ANY        \
  75   if ((a & 1) == 0)             \
  76     PROCESS_TWOS_EVEN;
  77 #endif
  78
  79 /* Process one bit arithmetically, then a simple loop.  This cuts the loop
  80    condition down to a 25/75 chance, which should branch predict better.
  81    The CPU will need a reasonable variable left shift.  */
  82 #if JACOBI_BASE_METHOD == 3
  83 #define PROCESS_TWOS_EVEN               \
  84   {                                     \
  85     int  two, mask, shift;              \
  86                                         \
  87     two = JACOBI_TWO_U_BIT1 (b);        \
  88     mask = (~a & 2);                    \
  89     a >>= 1;                            \
  90                                         \
  91     shift = (~a & 1);                   \
  92     a >>= shift;                        \
  93     result_bit1 ^= two ^ (two & mask);  \
  94                                         \
  95     while ((a & 1) == 0)                \
  96       {                                 \
  97         a >>= 1;                        \
  98         result_bit1 ^= two;             \
  99         ASSERT (a != 0);                \
 100       }                                 \
 101   }
 102 #define PROCESS_TWOS_ANY                \
 103   {                                     \
 104     int  two, mask, shift;              \
 105                                         \
 106     two = JACOBI_TWO_U_BIT1 (b);        \
 107     shift = (~a & 1);                   \
 108     a >>= shift;                        \
 109                                         \
 110     mask = shift << 1;                  \
 111     result_bit1 ^= (two & mask);        \
 112                                         \
 113     while ((a & 1) == 0)                \
 114       {                                 \
 115         a >>= 1;                        \
 116         result_bit1 ^= two;             \
 117         ASSERT (a != 0);                \
 118       }                                 \
 119   }
 120 #endif
 121
 122 #if JACOBI_BASE_METHOD < 4
 123 /* Calculate the value of the Jacobi symbol (a/b) of two mp_limb_t's, but
 124    with a restricted range of inputs accepted, namely b>1, b odd.
 125
 126    The initial result_bit1 is taken as a parameter for the convenience of
 127    mpz_kronecker_ui() et al.  The sign changes both here and in those
 128    routines accumulate nicely in bit 1, see the JACOBI macros.
 129
 130    The return value here is the normal +1, 0, or -1.  Note that +1 and -1
 131    have bit 1 in the "BIT1" sense, which could be useful if the caller is
 132    accumulating it into some extended calculation.
 133
 134    Duplicating the loop body to avoid the MP_LIMB_T_SWAP(a,b) would be
 135    possible, but a couple of tests suggest it's not a significant speedup,
 136    and may even be a slowdown, so what's here is good enough for now. */
 137
 138 int
 139 mpn_jacobi_base (mp_limb_t a, mp_limb_t b, int result_bit1)
 140 {
 141   ASSERT (b & 1);  /* b odd */
 142   ASSERT (b != 1);
 143
 144   if (a == 0)
 145     return 0;
 146
 147   PROCESS_TWOS_ANY;
 148   if (a == 1)
 149     goto done;
 150
 151   if (a >= b)
 152     goto a_gt_b;
 153
 154   for (;;)
 155     {
 156       result_bit1 ^= JACOBI_RECIP_UU_BIT1 (a, b);
 157       MP_LIMB_T_SWAP (a, b);
 158
 159     a_gt_b:
 160       do
 161         {
 162           /* working on (a/b), a,b odd, a>=b */
 163           ASSERT (a & 1);
 164           ASSERT (b & 1);
 165           ASSERT (a >= b);
 166
 167           if ((a -= b) == 0)
 168             return 0;
 169
 170           PROCESS_TWOS_EVEN;
 171           if (a == 1)
 172             goto done;
 173         }
 174       while (a >= b);
 175     }
 176
 177  done:
 178   return JACOBI_BIT1_TO_PN (result_bit1);
 179 }
 180 #endif
 181
 182 #if JACOBI_BASE_METHOD == 4
 183 /* Computes (a/b) for odd b > 1 and any a. The initial bit is taken as a
 184  * parameter. We have no need for the convention that the sign is in
 185  * bit 1, internally we use bit 0. */
 186
 187 /* FIXME: Could try table-based count_trailing_zeros. */
 188 int
 189 mpn_jacobi_base (mp_limb_t a, mp_limb_t b, int bit)
 190 {
 191   int c;
 192
 193   ASSERT (b & 1);
 194   ASSERT (b > 1);
 195
 196   if (a == 0)
 197     /* This is the only line which depends on b > 1 */
 198     return 0;
 199
 200   bit >>= 1;
 201
 202   /* Below, we represent a and b shifted right so that the least
 203      significant one bit is implicit. */
 204
 205   b >>= 1;
 206
 207   count_trailing_zeros (c, a);
 208   bit ^= c & (b ^ (b >> 1));
 209
 210   /* We may have c==GMP_LIMB_BITS-1, so we can't use a>>c+1. */
 211   a >>= c;
 212   a >>= 1;
 213
 214   do
 215     {
 216       mp_limb_t t = a - b;
 217       mp_limb_t bgta = LIMB_HIGHBIT_TO_MASK (t);
 218
 219       if (t == 0)
 220         return 0;
 221
 222       /* If b > a, invoke reciprocity */
 223       bit ^= (bgta & a & b);
 224
 225       /* b <-- min (a, b) */
 226       b += (bgta & t);
 227
 228       /* a <-- |a - b| */
 229       a = (t ^ bgta) - bgta;
 230
 231       /* Number of trailing zeros is the same no matter if we look at
 232        * t or a, but using t gives more parallelism. */
 233       count_trailing_zeros (c, t);
 234       c ++;
 235       /* (2/b) = -1 if b = 3 or 5 mod 8 */
 236       bit ^= c & (b ^ (b >> 1));
 237       a >>= c;
 238     }
 239   while (b > 0);
 240
 241   return 1-2*(bit & 1);
 242 }
 243 #endif /* JACOBI_BASE_METHOD == 4 */