Back to index

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