ghc/rts/gmp/mpn/generic/divrem_1.c

   1 /* mpn_divrem_1(quot_ptr, qsize, dividend_ptr, dividend_size, divisor_limb) --
   2    Divide (DIVIDEND_PTR,,DIVIDEND_SIZE) by DIVISOR_LIMB.
   3    Write DIVIDEND_SIZE limbs of quotient at QUOT_PTR.
   4    Return the single-limb remainder.
   5    There are no constraints on the value of the divisor.
   6
   7    QUOT_PTR and DIVIDEND_PTR might point to the same limb.
   8
   9 Copyright (C) 1991, 1993, 1994, 1996, 1998, 1999, 2000 Free Software
  10 Foundation, Inc.
  11
  12 This file is part of the GNU MP Library.
  13
  14 The GNU MP Library is free software; you can redistribute it and/or modify
  15 it under the terms of the GNU Lesser General Public License as published by
  16 the Free Software Foundation; either version 2.1 of the License, or (at your
  17 option) any later version.
  18
  19 The GNU MP Library is distributed in the hope that it will be useful, but
  20 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  21 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  22 License for more details.
  23
  24 You should have received a copy of the GNU Lesser General Public License
  25 along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
  26 the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
  27 MA 02111-1307, USA. */
  28
  29 #include "gmp.h"
  30 #include "gmp-impl.h"
  31 #include "longlong.h"
  32
  33
  34
  35 /* __gmpn_divmod_1_internal(quot_ptr,dividend_ptr,dividend_size,divisor_limb)
  36    Divide (DIVIDEND_PTR,,DIVIDEND_SIZE) by DIVISOR_LIMB.
  37    Write DIVIDEND_SIZE limbs of quotient at QUOT_PTR.
  38    Return the single-limb remainder.
  39    There are no constraints on the value of the divisor.
  40
  41    QUOT_PTR and DIVIDEND_PTR might point to the same limb. */
  42
  43 #ifndef UMUL_TIME
  44 #define UMUL_TIME 1
  45 #endif
  46
  47 #ifndef UDIV_TIME
  48 #define UDIV_TIME UMUL_TIME
  49 #endif
  50
  51 static mp_limb_t
  52 #if __STDC__
  53 __gmpn_divmod_1_internal (mp_ptr quot_ptr,
  54               mp_srcptr dividend_ptr, mp_size_t dividend_size,
  55               mp_limb_t divisor_limb)
  56 #else
  57 __gmpn_divmod_1_internal (quot_ptr, dividend_ptr, dividend_size, divisor_limb)
  58      mp_ptr quot_ptr;
  59      mp_srcptr dividend_ptr;
  60      mp_size_t dividend_size;
  61      mp_limb_t divisor_limb;
  62 #endif
  63 {
  64   mp_size_t i;
  65   mp_limb_t n1, n0, r;
  66   int dummy;
  67
  68   /* ??? Should this be handled at all?  Rely on callers?  */
  69   if (dividend_size == 0)
  70     return 0;
  71
  72   /* If multiplication is much faster than division, and the
  73      dividend is large, pre-invert the divisor, and use
  74      only multiplications in the inner loop.  */
  75
  76   /* This test should be read:
  77        Does it ever help to use udiv_qrnnd_preinv?
  78          && Does what we save compensate for the inversion overhead?  */
  79   if (UDIV_TIME > (2 * UMUL_TIME + 6)
  80       && (UDIV_TIME - (2 * UMUL_TIME + 6)) * dividend_size > UDIV_TIME)
  81     {
  82       int normalization_steps;
  83
  84       count_leading_zeros (normalization_steps, divisor_limb);
  85       if (normalization_steps != 0)
  86         {
  87           mp_limb_t divisor_limb_inverted;
  88
  89           divisor_limb <<= normalization_steps;
  90           invert_limb (divisor_limb_inverted, divisor_limb);
  91
  92           n1 = dividend_ptr[dividend_size - 1];
  93           r = n1 >> (BITS_PER_MP_LIMB - normalization_steps);
  94
  95           /* Possible optimization:
  96              if (r == 0
  97              && divisor_limb > ((n1 << normalization_steps)
  98                              | (dividend_ptr[dividend_size - 2] >> ...)))
  99              ...one division less... */
 100
 101           for (i = dividend_size - 2; i >= 0; i--)
 102             {
 103               n0 = dividend_ptr[i];
 104               udiv_qrnnd_preinv (quot_ptr[i + 1], r, r,
 105                                  ((n1 << normalization_steps)
 106                                   | (n0 >> (BITS_PER_MP_LIMB - normalization_steps))),
 107                                  divisor_limb, divisor_limb_inverted);
 108               n1 = n0;
 109             }
 110           udiv_qrnnd_preinv (quot_ptr[0], r, r,
 111                              n1 << normalization_steps,
 112                              divisor_limb, divisor_limb_inverted);
 113           return r >> normalization_steps;
 114         }
 115       else
 116         {
 117           mp_limb_t divisor_limb_inverted;
 118
 119           invert_limb (divisor_limb_inverted, divisor_limb);
 120
 121           i = dividend_size - 1;
 122           r = dividend_ptr[i];
 123
 124           if (r >= divisor_limb)
 125             r = 0;
 126           else
 127             {
 128               quot_ptr[i] = 0;
 129               i--;
 130             }
 131
 132           for (; i >= 0; i--)
 133             {
 134               n0 = dividend_ptr[i];
 135               udiv_qrnnd_preinv (quot_ptr[i], r, r,
 136                                  n0, divisor_limb, divisor_limb_inverted);
 137             }
 138           return r;
 139         }
 140     }
 141   else
 142     {
 143       if (UDIV_NEEDS_NORMALIZATION)
 144         {
 145           int normalization_steps;
 146
 147           count_leading_zeros (normalization_steps, divisor_limb);
 148           if (normalization_steps != 0)
 149             {
 150               divisor_limb <<= normalization_steps;
 151
 152               n1 = dividend_ptr[dividend_size - 1];
 153               r = n1 >> (BITS_PER_MP_LIMB - normalization_steps);
 154
 155               /* Possible optimization:
 156                  if (r == 0
 157                  && divisor_limb > ((n1 << normalization_steps)
 158                                  | (dividend_ptr[dividend_size - 2] >> ...)))
 159                  ...one division less... */
 160
 161               for (i = dividend_size - 2; i >= 0; i--)
 162                 {
 163                   n0 = dividend_ptr[i];
 164                   udiv_qrnnd (quot_ptr[i + 1], r, r,
 165                               ((n1 << normalization_steps)
 166                                | (n0 >> (BITS_PER_MP_LIMB - normalization_steps))),
 167                               divisor_limb);
 168                   n1 = n0;
 169                 }
 170               udiv_qrnnd (quot_ptr[0], r, r,
 171                           n1 << normalization_steps,
 172                           divisor_limb);
 173               return r >> normalization_steps;
 174             }
 175         }
 176       /* No normalization needed, either because udiv_qrnnd doesn't require
 177          it, or because DIVISOR_LIMB is already normalized.  */
 178
 179       i = dividend_size - 1;
 180       r = dividend_ptr[i];
 181
 182       if (r >= divisor_limb)
 183         r = 0;
 184       else
 185         {
 186           quot_ptr[i] = 0;
 187           i--;
 188         }
 189
 190       for (; i >= 0; i--)
 191         {
 192           n0 = dividend_ptr[i];
 193           udiv_qrnnd (quot_ptr[i], r, r, n0, divisor_limb);
 194         }
 195       return r;
 196     }
 197 }
 198
 199
 200
 201 mp_limb_t
 202 #if __STDC__
 203 mpn_divrem_1 (mp_ptr qp, mp_size_t qxn,
 204               mp_srcptr np, mp_size_t nn,
 205               mp_limb_t d)
 206 #else
 207 mpn_divrem_1 (qp, qxn, np, nn, d)
 208      mp_ptr qp;
 209      mp_size_t qxn;
 210      mp_srcptr np;
 211      mp_size_t nn;
 212      mp_limb_t d;
 213 #endif
 214 {
 215   mp_limb_t rlimb;
 216   mp_size_t i;
 217
 218   /* Develop integer part of quotient.  */
 219   rlimb = __gmpn_divmod_1_internal (qp + qxn, np, nn, d);
 220
 221   /* Develop fraction part of quotient.  This is not as fast as it should;
 222      the preinvert stuff from __gmpn_divmod_1_internal ought to be used here
 223      too.  */
 224   if (UDIV_NEEDS_NORMALIZATION)
 225     {
 226       int normalization_steps;
 227
 228       count_leading_zeros (normalization_steps, d);
 229       if (normalization_steps != 0)
 230         {
 231           d <<= normalization_steps;
 232           rlimb <<= normalization_steps;
 233
 234           for (i = qxn - 1; i >= 0; i--)
 235             udiv_qrnnd (qp[i], rlimb, rlimb, 0, d);
 236
 237           return rlimb >> normalization_steps;
 238         }
 239       else
 240         /* fall out */
 241         ;
 242     }
 243
 244   for (i = qxn - 1; i >= 0; i--)
 245     udiv_qrnnd (qp[i], rlimb, rlimb, 0, d);
 246
 247   return rlimb;
 248 }