Back to index

glibc  2.9
mod_1.c
Go to the documentation of this file.
00001 /* mpn_mod_1(dividend_ptr, dividend_size, divisor_limb) --
00002    Divide (DIVIDEND_PTR,,DIVIDEND_SIZE) by DIVISOR_LIMB.
00003    Return the single-limb remainder.
00004    There are no constraints on the value of the divisor.
00005 
00006 Copyright (C) 1991, 1993, 1994, Free Software Foundation, Inc.
00007 
00008 This file is part of the GNU MP Library.
00009 
00010 The GNU MP Library is free software; you can redistribute it and/or modify
00011 it under the terms of the GNU Lesser General Public License as published by
00012 the Free Software Foundation; either version 2.1 of the License, or (at your
00013 option) any later version.
00014 
00015 The GNU MP Library is distributed in the hope that it will be useful, but
00016 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
00017 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
00018 License for more details.
00019 
00020 You should have received a copy of the GNU Lesser General Public License
00021 along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
00022 the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
00023 MA 02111-1307, USA. */
00024 
00025 #include <gmp.h>
00026 #include "gmp-impl.h"
00027 #include "longlong.h"
00028 
00029 #ifndef UMUL_TIME
00030 #define UMUL_TIME 1
00031 #endif
00032 
00033 #ifndef UDIV_TIME
00034 #define UDIV_TIME UMUL_TIME
00035 #endif
00036 
00037 /* FIXME: We should be using invert_limb (or invert_normalized_limb)
00038    here (not udiv_qrnnd).  */
00039 
00040 mp_limb_t
00041 #if __STDC__
00042 mpn_mod_1 (mp_srcptr dividend_ptr, mp_size_t dividend_size,
00043           mp_limb_t divisor_limb)
00044 #else
00045 mpn_mod_1 (dividend_ptr, dividend_size, divisor_limb)
00046      mp_srcptr dividend_ptr;
00047      mp_size_t dividend_size;
00048      mp_limb_t divisor_limb;
00049 #endif
00050 {
00051   mp_size_t i;
00052   mp_limb_t n1, n0, r;
00053   mp_limb_t dummy;
00054 
00055   /* Botch: Should this be handled at all?  Rely on callers?  */
00056   if (dividend_size == 0)
00057     return 0;
00058 
00059   /* If multiplication is much faster than division, and the
00060      dividend is large, pre-invert the divisor, and use
00061      only multiplications in the inner loop.  */
00062 
00063   /* This test should be read:
00064        Does it ever help to use udiv_qrnnd_preinv?
00065         && Does what we save compensate for the inversion overhead?  */
00066   if (UDIV_TIME > (2 * UMUL_TIME + 6)
00067       && (UDIV_TIME - (2 * UMUL_TIME + 6)) * dividend_size > UDIV_TIME)
00068     {
00069       int normalization_steps;
00070 
00071       count_leading_zeros (normalization_steps, divisor_limb);
00072       if (normalization_steps != 0)
00073        {
00074          mp_limb_t divisor_limb_inverted;
00075 
00076          divisor_limb <<= normalization_steps;
00077 
00078          /* Compute (2**2N - 2**N * DIVISOR_LIMB) / DIVISOR_LIMB.  The
00079             result is a (N+1)-bit approximation to 1/DIVISOR_LIMB, with the
00080             most significant bit (with weight 2**N) implicit.  */
00081 
00082          /* Special case for DIVISOR_LIMB == 100...000.  */
00083          if (divisor_limb << 1 == 0)
00084            divisor_limb_inverted = ~(mp_limb_t) 0;
00085          else
00086            udiv_qrnnd (divisor_limb_inverted, dummy,
00087                      -divisor_limb, 0, divisor_limb);
00088 
00089          n1 = dividend_ptr[dividend_size - 1];
00090          r = n1 >> (BITS_PER_MP_LIMB - normalization_steps);
00091 
00092          /* Possible optimization:
00093             if (r == 0
00094             && divisor_limb > ((n1 << normalization_steps)
00095                           | (dividend_ptr[dividend_size - 2] >> ...)))
00096             ...one division less... */
00097 
00098          for (i = dividend_size - 2; i >= 0; i--)
00099            {
00100              n0 = dividend_ptr[i];
00101              udiv_qrnnd_preinv (dummy, r, r,
00102                              ((n1 << normalization_steps)
00103                               | (n0 >> (BITS_PER_MP_LIMB - normalization_steps))),
00104                              divisor_limb, divisor_limb_inverted);
00105              n1 = n0;
00106            }
00107          udiv_qrnnd_preinv (dummy, r, r,
00108                           n1 << normalization_steps,
00109                           divisor_limb, divisor_limb_inverted);
00110          return r >> normalization_steps;
00111        }
00112       else
00113        {
00114          mp_limb_t divisor_limb_inverted;
00115 
00116          /* Compute (2**2N - 2**N * DIVISOR_LIMB) / DIVISOR_LIMB.  The
00117             result is a (N+1)-bit approximation to 1/DIVISOR_LIMB, with the
00118             most significant bit (with weight 2**N) implicit.  */
00119 
00120          /* Special case for DIVISOR_LIMB == 100...000.  */
00121          if (divisor_limb << 1 == 0)
00122            divisor_limb_inverted = ~(mp_limb_t) 0;
00123          else
00124            udiv_qrnnd (divisor_limb_inverted, dummy,
00125                      -divisor_limb, 0, divisor_limb);
00126 
00127          i = dividend_size - 1;
00128          r = dividend_ptr[i];
00129 
00130          if (r >= divisor_limb)
00131            r = 0;
00132          else
00133            i--;
00134 
00135          for (; i >= 0; i--)
00136            {
00137              n0 = dividend_ptr[i];
00138              udiv_qrnnd_preinv (dummy, r, r,
00139                              n0, divisor_limb, divisor_limb_inverted);
00140            }
00141          return r;
00142        }
00143     }
00144   else
00145     {
00146       if (UDIV_NEEDS_NORMALIZATION)
00147        {
00148          int normalization_steps;
00149 
00150          count_leading_zeros (normalization_steps, divisor_limb);
00151          if (normalization_steps != 0)
00152            {
00153              divisor_limb <<= normalization_steps;
00154 
00155              n1 = dividend_ptr[dividend_size - 1];
00156              r = n1 >> (BITS_PER_MP_LIMB - normalization_steps);
00157 
00158              /* Possible optimization:
00159                if (r == 0
00160                && divisor_limb > ((n1 << normalization_steps)
00161                              | (dividend_ptr[dividend_size - 2] >> ...)))
00162                ...one division less... */
00163 
00164              for (i = dividend_size - 2; i >= 0; i--)
00165               {
00166                 n0 = dividend_ptr[i];
00167                 udiv_qrnnd (dummy, r, r,
00168                            ((n1 << normalization_steps)
00169                             | (n0 >> (BITS_PER_MP_LIMB - normalization_steps))),
00170                            divisor_limb);
00171                 n1 = n0;
00172               }
00173              udiv_qrnnd (dummy, r, r,
00174                        n1 << normalization_steps,
00175                        divisor_limb);
00176              return r >> normalization_steps;
00177            }
00178        }
00179       /* No normalization needed, either because udiv_qrnnd doesn't require
00180         it, or because DIVISOR_LIMB is already normalized.  */
00181 
00182       i = dividend_size - 1;
00183       r = dividend_ptr[i];
00184 
00185       if (r >= divisor_limb)
00186        r = 0;
00187       else
00188        i--;
00189 
00190       for (; i >= 0; i--)
00191        {
00192          n0 = dividend_ptr[i];
00193          udiv_qrnnd (dummy, r, r, n0, divisor_limb);
00194        }
00195       return r;
00196     }
00197 }