contrib/gmp/mpn/generic/mod_1.c

   1 /* mpn_mod_1(dividend_ptr, dividend_size, divisor_limb) --
   2    Divide (DIVIDEND_PTR,,DIVIDEND_SIZE) by DIVISOR_LIMB.
   3    Return the single-limb remainder.
   4    There are no constraints on the value of the divisor.
   5
   6 Copyright 1991, 1993, 1994, 1999, 2000, 2002, 2007, 2008, 2009 Free
   7 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 3 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.  If not, see http://www.gnu.org/licenses/.  */
  23
  24 #include "gmp.h"
  25 #include "gmp-impl.h"
  26 #include "longlong.h"
  27
  28
  29 /* The size where udiv_qrnnd_preinv should be used rather than udiv_qrnnd,
  30    meaning the quotient size where that should happen, the quotient size
  31    being how many udiv divisions will be done.
  32
  33    The default is to use preinv always, CPUs where this doesn't suit have
  34    tuned thresholds.  Note in particular that preinv should certainly be
  35    used if that's the only division available (USE_PREINV_ALWAYS).  */
  36
  37 #ifndef MOD_1_NORM_THRESHOLD
  38 #define MOD_1_NORM_THRESHOLD  0
  39 #endif
  40
  41 #ifndef MOD_1_UNNORM_THRESHOLD
  42 #define MOD_1_UNNORM_THRESHOLD  0
  43 #endif
  44
  45 #ifndef MOD_1_1_THRESHOLD
  46 #define MOD_1_1_THRESHOLD  MP_SIZE_T_MAX /* default is not to use mpn_mod_1s */
  47 #endif
  48
  49 #ifndef MOD_1_2_THRESHOLD
  50 #define MOD_1_2_THRESHOLD  10
  51 #endif
  52
  53 #ifndef MOD_1_4_THRESHOLD
  54 #define MOD_1_4_THRESHOLD  120
  55 #endif
  56
  57
  58 /* The comments in mpn/generic/divrem_1.c apply here too.
  59
  60    As noted in the algorithms section of the manual, the shifts in the loop
  61    for the unnorm case can be avoided by calculating r = a%(d*2^n), followed
  62    by a final (r*2^n)%(d*2^n).  In fact if it happens that a%(d*2^n) can
  63    skip a division where (a*2^n)%(d*2^n) can't then there's the same number
  64    of divide steps, though how often that happens depends on the assumed
  65    distributions of dividend and divisor.  In any case this idea is left to
  66    CPU specific implementations to consider.  */
  67
  68 static mp_limb_t
  69 mpn_mod_1_unnorm (mp_srcptr up, mp_size_t un, mp_limb_t d)
  70 {
  71   mp_size_t  i;
  72   mp_limb_t  n1, n0, r;
  73   mp_limb_t  dummy;
  74   int cnt;
  75
  76   ASSERT (un > 0);
  77   ASSERT (d != 0);
  78
  79   d <<= GMP_NAIL_BITS;
  80
  81   /* Skip a division if high < divisor.  Having the test here before
  82      normalizing will still skip as often as possible.  */
  83   r = up[un - 1] << GMP_NAIL_BITS;
  84   if (r < d)
  85     {
  86       r >>= GMP_NAIL_BITS;
  87       un--;
  88       if (un == 0)
  89         return r;
  90     }
  91   else
  92     r = 0;
  93
  94   /* If udiv_qrnnd doesn't need a normalized divisor, can use the simple
  95      code above. */
  96   if (! UDIV_NEEDS_NORMALIZATION
  97       && BELOW_THRESHOLD (un, MOD_1_UNNORM_THRESHOLD))
  98     {
  99       for (i = un - 1; i >= 0; i--)
 100         {
 101           n0 = up[i] << GMP_NAIL_BITS;
 102           udiv_qrnnd (dummy, r, r, n0, d);
 103           r >>= GMP_NAIL_BITS;
 104         }
 105       return r;
 106     }
 107
 108   count_leading_zeros (cnt, d);
 109   d <<= cnt;
 110
 111   n1 = up[un - 1] << GMP_NAIL_BITS;
 112   r = (r << cnt) | (n1 >> (GMP_LIMB_BITS - cnt));
 113
 114   if (UDIV_NEEDS_NORMALIZATION
 115       && BELOW_THRESHOLD (un, MOD_1_UNNORM_THRESHOLD))
 116     {
 117       for (i = un - 2; i >= 0; i--)
 118         {
 119           n0 = up[i] << GMP_NAIL_BITS;
 120           udiv_qrnnd (dummy, r, r,
 121                       (n1 << cnt) | (n0 >> (GMP_NUMB_BITS - cnt)),
 122                       d);
 123           r >>= GMP_NAIL_BITS;
 124           n1 = n0;
 125         }
 126       udiv_qrnnd (dummy, r, r, n1 << cnt, d);
 127       r >>= GMP_NAIL_BITS;
 128       return r >> cnt;
 129     }
 130   else
 131     {
 132       mp_limb_t inv;
 133       invert_limb (inv, d);
 134
 135       for (i = un - 2; i >= 0; i--)
 136         {
 137           n0 = up[i] << GMP_NAIL_BITS;
 138           udiv_qrnnd_preinv (dummy, r, r,
 139                              (n1 << cnt) | (n0 >> (GMP_NUMB_BITS - cnt)),
 140                              d, inv);
 141           r >>= GMP_NAIL_BITS;
 142           n1 = n0;
 143         }
 144       udiv_qrnnd_preinv (dummy, r, r, n1 << cnt, d, inv);
 145       r >>= GMP_NAIL_BITS;
 146       return r >> cnt;
 147     }
 148 }
 149
 150 static mp_limb_t
 151 mpn_mod_1_norm (mp_srcptr up, mp_size_t un, mp_limb_t d)
 152 {
 153   mp_size_t  i;
 154   mp_limb_t  n0, r;
 155   mp_limb_t  dummy;
 156
 157   ASSERT (un > 0);
 158
 159   d <<= GMP_NAIL_BITS;
 160
 161   ASSERT (d & GMP_LIMB_HIGHBIT);
 162
 163   /* High limb is initial remainder, possibly with one subtract of
 164      d to get r<d.  */
 165   r = up[un - 1] << GMP_NAIL_BITS;
 166   if (r >= d)
 167     r -= d;
 168   r >>= GMP_NAIL_BITS;
 169   un--;
 170   if (un == 0)
 171     return r;
 172
 173   if (BELOW_THRESHOLD (un, MOD_1_NORM_THRESHOLD))
 174     {
 175       for (i = un - 1; i >= 0; i--)
 176         {
 177           n0 = up[i] << GMP_NAIL_BITS;
 178           udiv_qrnnd (dummy, r, r, n0, d);
 179           r >>= GMP_NAIL_BITS;
 180         }
 181       return r;
 182     }
 183   else
 184     {
 185       mp_limb_t  inv;
 186       invert_limb (inv, d);
 187       for (i = un - 1; i >= 0; i--)
 188         {
 189           n0 = up[i] << GMP_NAIL_BITS;
 190           udiv_qrnnd_preinv (dummy, r, r, n0, d, inv);
 191           r >>= GMP_NAIL_BITS;
 192         }
 193       return r;
 194     }
 195 }
 196
 197 mp_limb_t
 198 mpn_mod_1 (mp_srcptr ap, mp_size_t n, mp_limb_t b)
 199 {
 200   ASSERT (n >= 0);
 201   ASSERT (b != 0);
 202
 203   /* Should this be handled at all?  Rely on callers?  Note un==0 is currently
 204      required by mpz/fdiv_r_ui.c and possibly other places.  */
 205   if (n == 0)
 206     return 0;
 207
 208   if (UNLIKELY ((b & GMP_NUMB_HIGHBIT) != 0))
 209     {
 210       /* The functions below do not handle this large divisor.  */
 211       return mpn_mod_1_norm (ap, n, b);
 212     }
 213   else if (BELOW_THRESHOLD (n, MOD_1_1_THRESHOLD))
 214     {
 215       return mpn_mod_1_unnorm (ap, n, b);
 216     }
 217   else if (BELOW_THRESHOLD (n, MOD_1_2_THRESHOLD))
 218     {
 219       mp_limb_t pre[4];
 220       mpn_mod_1s_1p_cps (pre, b);
 221       return mpn_mod_1s_1p (ap, n, b << pre[1], pre);
 222     }
 223   else if (BELOW_THRESHOLD (n, MOD_1_4_THRESHOLD) || UNLIKELY (b > GMP_NUMB_MASK / 4))
 224     {
 225       mp_limb_t pre[5];
 226       mpn_mod_1s_2p_cps (pre, b);
 227       return mpn_mod_1s_2p (ap, n, b << pre[1], pre);
 228     }
 229   else
 230     {
 231       mp_limb_t pre[7];
 232       mpn_mod_1s_4p_cps (pre, b);
 233       return mpn_mod_1s_4p (ap, n, b << pre[1], pre);
 234     }
 235 }