ghc/runtime/gmp/mpn_dm_1.c

   1 /* mpn_divmod_1(quot_ptr, 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 Free Software Foundation, Inc.
  10
  11 This file is part of the GNU MP Library.
  12
  13 The GNU MP Library is free software; you can redistribute it and/or modify
  14 it under the terms of the GNU General Public License as published by
  15 the Free Software Foundation; either version 2, or (at your option)
  16 any later version.
  17
  18 The GNU MP Library is distributed in the hope that it will be useful,
  19 but WITHOUT ANY WARRANTY; without even the implied warranty of
  20 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  21 GNU General Public License for more details.
  22
  23 You should have received a copy of the GNU General Public License
  24 along with the GNU MP Library; see the file COPYING.  If not, write to
  25 the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.  */
  26
  27 #include "gmp.h"
  28 #include "gmp-impl.h"
  29 #include "longlong.h"
  30
  31 #ifndef UMUL_TIME
  32 #define UMUL_TIME 1
  33 #endif
  34
  35 #ifndef UDIV_TIME
  36 #define UDIV_TIME UMUL_TIME
  37 #endif
  38
  39 #if UDIV_TIME > 2 * UMUL_TIME
  40 #undef UDIV_NEEDS_NORMALIZATION
  41 #define UDIV_NEEDS_NORMALIZATION 1
  42 #endif
  43
  44 #define udiv_qrnnd_preinv(q, r, nh, nl, d, di) \
  45   do {                                                                  \
  46     unsigned long int _q, _ql, _r;                                      \
  47     unsigned long int _xh, _xl;                                         \
  48     umul_ppmm (_q, _ql, (nh), (di));                                    \
  49     _q += (nh);                 /* DI is 2**BITS_PER_MP_LIMB too small.  */\
  50     umul_ppmm (_xh, _xl, _q, (d));                                      \
  51     sub_ddmmss (_xh, _r, (nh), (nl), _xh, _xl);                         \
  52     if (_xh != 0)                                                       \
  53       {                                                                 \
  54         sub_ddmmss (_xh, _r, _xh, _r, 0, (d));                          \
  55         _q += 1;                                                        \
  56         if (_xh != 0)                                                   \
  57           {                                                             \
  58             sub_ddmmss (_xh, _r, _xh, _r, 0, (d));                      \
  59             _q += 1;                                                    \
  60           }                                                             \
  61       }                                                                 \
  62     if (_r >= (d))                                                      \
  63       {                                                                 \
  64         _r -= (d);                                                      \
  65         _q += 1;                                                        \
  66       }                                                                 \
  67     (r) = _r;                                                           \
  68     (q) = _q;                                                           \
  69   } while (0)
  70
  71 mp_limb
  72 #ifdef __STDC__
  73 mpn_divmod_1 (mp_ptr quot_ptr,
  74               mp_srcptr dividend_ptr, mp_size dividend_size,
  75               unsigned long int divisor_limb)
  76 #else
  77 mpn_divmod_1 (quot_ptr, dividend_ptr, dividend_size, divisor_limb)
  78      mp_ptr quot_ptr;
  79      mp_srcptr dividend_ptr;
  80      mp_size dividend_size;
  81      unsigned long int divisor_limb;
  82 #endif
  83 {
  84   mp_size i;
  85   mp_limb n1, n0, r;
  86
  87   /* Botch: Should this be handled at all?  Rely on callers?  */
  88   if (dividend_size == 0)
  89     return 0;
  90
  91   if (UDIV_NEEDS_NORMALIZATION)
  92     {
  93       int normalization_steps;
  94
  95       count_leading_zeros (normalization_steps, divisor_limb);
  96       if (normalization_steps != 0)
  97         {
  98           divisor_limb <<= normalization_steps;
  99
 100           n1 = dividend_ptr[dividend_size - 1];
 101           r = n1 >> (BITS_PER_MP_LIMB - normalization_steps);
 102
 103           /* Possible optimization:
 104              if (r == 0
 105              && divisor_limb > ((n1 << normalization_steps)
 106                              | (dividend_ptr[dividend_size - 2] >> ...)))
 107              ...one division less...
 108              [Don't forget to zero most sign. quotient limb!]  */
 109
 110           /* If multiplication is much faster than division, and the
 111              dividend is large, pre-invert the divisor, and use
 112              only multiplications in the inner loop.  */
 113           if (UDIV_TIME > 2 * UMUL_TIME && dividend_size >= 4)
 114             {
 115               mp_limb divisor_limb_inverted;
 116               int dummy;
 117               /* Compute (2**64 - 2**32 * DIVISOR_LIMB) / DIVISOR_LIMB.
 118                  The result is an 33-bit approximation to 1/DIVISOR_LIMB,
 119                  with the most significant bit (weight 2**32) implicit.  */
 120
 121               /* Special case for DIVISOR_LIMB == 100...000.  */
 122               if (divisor_limb << 1 == 0)
 123                 divisor_limb_inverted = ~0;
 124               else
 125                 udiv_qrnnd (divisor_limb_inverted, dummy,
 126                             -divisor_limb, 0, divisor_limb);
 127
 128               for (i = dividend_size - 2; i >= 0; i--)
 129                 {
 130                   n0 = dividend_ptr[i];
 131                   udiv_qrnnd_preinv (quot_ptr[i + 1], r, r,
 132                                      ((n1 << normalization_steps)
 133                                       | (n0 >> (BITS_PER_MP_LIMB - normalization_steps))),
 134                                      divisor_limb, divisor_limb_inverted);
 135                   n1 = n0;
 136                 }
 137               udiv_qrnnd_preinv (quot_ptr[0], r, r,
 138                           n1 << normalization_steps,
 139                           divisor_limb, divisor_limb_inverted);
 140               return r >> normalization_steps;
 141             }
 142           else
 143             {
 144               for (i = dividend_size - 2; i >= 0; i--)
 145                 {
 146                   n0 = dividend_ptr[i];
 147                   udiv_qrnnd (quot_ptr[i + 1], r, r,
 148                               ((n1 << normalization_steps)
 149                                | (n0 >> (BITS_PER_MP_LIMB - normalization_steps))),
 150                               divisor_limb);
 151                   n1 = n0;
 152                 }
 153               udiv_qrnnd (quot_ptr[0], r, r,
 154                           n1 << normalization_steps,
 155                           divisor_limb);
 156               return r >> normalization_steps;
 157             }
 158         }
 159     }
 160
 161   /* No normalization needed, either because udiv_qrnnd doesn't require
 162      it, or because DIVISOR_LIMB is already normalized.  */
 163
 164   i = dividend_size - 1;
 165   r = dividend_ptr[i];
 166
 167   if (r >= divisor_limb)
 168     {
 169       r = 0;
 170     }
 171   else
 172     {
 173       /* Callers expect the quotient to be DIVIDEND_SIZE limbs.  Store
 174          a leading zero to make that expectation come true.  */
 175       quot_ptr[i] = 0;
 176       i--;
 177     }
 178
 179   for (; i >= 0; i--)
 180     {
 181       n0 = dividend_ptr[i];
 182       udiv_qrnnd (quot_ptr[i], r, r, n0, divisor_limb);
 183     }
 184   return r;
 185 }