rts/gmp/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 /*
   7 Copyright (C) 1999, 2000 Free Software Foundation, Inc.
   8
   9 This file is part of the GNU MP Library.
  10
  11 The GNU MP Library is free software; you can redistribute it and/or modify
  12 it under the terms of the GNU Lesser General Public License as published by
  13 the Free Software Foundation; either version 2.1 of the License, or (at your
  14 option) any later version.
  15
  16 The GNU MP Library is distributed in the hope that it will be useful, but
  17 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  18 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  19 License for more details.
  20
  21 You should have received a copy of the GNU Lesser General Public License
  22 along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
  23 the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
  24 MA 02111-1307, USA.  */
  25
  26 #include "gmp.h"
  27 #include "gmp-impl.h"
  28 #include "longlong.h"
  29
  30
  31 #if COUNT_TRAILING_ZEROS_TIME <= 7
  32 /* If count_trailing_zeros is fast, use it.
  33    K7 at 7 cycles and P6 at 2 are good here.  K6 at 12-27 and P5 at 18-42
  34    are not.  The default 15 in longlong.h is meant to mean not good here.  */
  35
  36 #define PROCESS_TWOS_ANY                                \
  37   {                                                     \
  38     mp_limb_t  twos;                                    \
  39     count_trailing_zeros (twos, a);                     \
  40     result_bit1 ^= JACOBI_TWOS_U_BIT1 (twos, b);        \
  41     a >>= twos;                                         \
  42   }
  43
  44 #define PROCESS_TWOS_EVEN  PROCESS_TWOS_ANY
  45
  46 #else
  47 /* Use a loop instead.  With "a" uniformly distributed there will usually be
  48    only a few trailing zeros.
  49
  50    Unfortunately the branch for the while loop here will be on a 50/50
  51    chance of a 1 or 0, which is bad for branch prediction.  */
  52
  53 #define PROCESS_TWOS_EVEN               \
  54   {                                     \
  55     int  two;                           \
  56     two = JACOBI_TWO_U_BIT1 (b);        \
  57     do                                  \
  58       {                                 \
  59         a >>= 1;                        \
  60         result_bit1 ^= two;             \
  61         ASSERT (a != 0);                \
  62       }                                 \
  63     while ((a & 1) == 0);               \
  64   }
  65
  66 #define PROCESS_TWOS_ANY        \
  67   if ((a & 1) == 0)             \
  68     PROCESS_TWOS_EVEN;
  69
  70 #endif
  71
  72
  73 /* Calculate the value of the Jacobi symbol (a/b) of two mp_limb_t's, but
  74    with a restricted range of inputs accepted, namely b>1, b odd, and a<=b.
  75
  76    The initial result_bit1 is taken as a parameter for the convenience of
  77    mpz_kronecker_zi_ui() et al.  The sign changes both here and in those
  78    routines accumulate nicely in bit 1, see the JACOBI macros.
  79
  80    The return value here is the normal +1, 0, or -1.  Note that +1 and -1
  81    have bit 1 in the "BIT1" sense, which could be useful if the caller is
  82    accumulating it into some extended calculation.
  83
  84    Duplicating the loop body to avoid the MP_LIMB_T_SWAP(a,b) would be
  85    possible, but a couple of tests suggest it's not a significant speedup,
  86    and may even be a slowdown, so what's here is good enough for now.
  87
  88    Future: The code doesn't demand a<=b actually, so maybe this could be
  89    relaxed.  All the places this is used currently call with a<=b though.  */
  90
  91 int
  92 #if __STDC__
  93 mpn_jacobi_base (mp_limb_t a, mp_limb_t b, int result_bit1)
  94 #else
  95 mpn_jacobi_base (a, b, result_bit1)
  96      mp_limb_t a;
  97      mp_limb_t b;
  98      int       result_bit1;
  99 #endif
 100 {
 101   ASSERT (b & 1);  /* b odd */
 102   ASSERT (b != 1);
 103   ASSERT (a <= b);
 104
 105   if (a == 0)
 106     return 0;
 107
 108   PROCESS_TWOS_ANY;
 109   if (a == 1)
 110     goto done;
 111
 112   for (;;)
 113     {
 114       result_bit1 ^= JACOBI_RECIP_UU_BIT1 (a, b);
 115       MP_LIMB_T_SWAP (a, b);
 116
 117       do
 118         {
 119           /* working on (a/b), a,b odd, a>=b */
 120           ASSERT (a & 1);
 121           ASSERT (b & 1);
 122           ASSERT (a >= b);
 123
 124           if ((a -= b) == 0)
 125             return 0;
 126
 127           PROCESS_TWOS_EVEN;
 128           if (a == 1)
 129             goto done;
 130         }
 131       while (a >= b);
 132     }
 133
 134  done:
 135   return JACOBI_BIT1_TO_PN (result_bit1);
 136 }