contrib/mpfr/src/div_ui.c

   1 /* mpfr_div_{ui,si} -- divide a floating-point number by a machine integer
   2
   3 Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
   4 Contributed by the AriC and Caramel projects, INRIA.
   5
   6 This file is part of the GNU MPFR Library.
   7
   8 The GNU MPFR Library is free software; you can redistribute it and/or modify
   9 it under the terms of the GNU Lesser General Public License as published by
  10 the Free Software Foundation; either version 3 of the License, or (at your
  11 option) any later version.
  12
  13 The GNU MPFR Library is distributed in the hope that it will be useful, but
  14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  15 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  16 License for more details.
  17
  18 You should have received a copy of the GNU Lesser General Public License
  19 along with the GNU MPFR Library; see the file COPYING.LESSER.  If not, see
  20 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
  21 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
  22
  23 #define MPFR_NEED_LONGLONG_H
  24 #include "mpfr-impl.h"
  25
  26 /* returns 0 if result exact, non-zero otherwise */
  27 int
  28 mpfr_div_ui (mpfr_ptr y, mpfr_srcptr x, unsigned long int u, mpfr_rnd_t rnd_mode)
  29 {
  30   long i;
  31   int sh;
  32   mp_size_t xn, yn, dif;
  33   mp_limb_t *xp, *yp, *tmp, c, d;
  34   mpfr_exp_t exp;
  35   int inexact, middle = 1, nexttoinf;
  36   MPFR_TMP_DECL(marker);
  37
  38   MPFR_LOG_FUNC
  39     (("x[%Pu]=%.*Rg u=%lu rnd=%d",
  40       mpfr_get_prec(x), mpfr_log_prec, x, u, rnd_mode),
  41      ("y[%Pu]=%.*Rg inexact=%d",
  42       mpfr_get_prec(y), mpfr_log_prec, y, inexact));
  43
  44   if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
  45     {
  46       if (MPFR_IS_NAN (x))
  47         {
  48           MPFR_SET_NAN (y);
  49           MPFR_RET_NAN;
  50         }
  51       else if (MPFR_IS_INF (x))
  52         {
  53           MPFR_SET_INF (y);
  54           MPFR_SET_SAME_SIGN (y, x);
  55           MPFR_RET (0);
  56         }
  57       else
  58         {
  59           MPFR_ASSERTD (MPFR_IS_ZERO(x));
  60           if (u == 0) /* 0/0 is NaN */
  61             {
  62               MPFR_SET_NAN(y);
  63               MPFR_RET_NAN;
  64             }
  65           else
  66             {
  67               MPFR_SET_ZERO(y);
  68               MPFR_SET_SAME_SIGN (y, x);
  69               MPFR_RET(0);
  70             }
  71         }
  72     }
  73   else if (MPFR_UNLIKELY (u <= 1))
  74     {
  75       if (u < 1)
  76         {
  77           /* x/0 is Inf since x != 0*/
  78           MPFR_SET_INF (y);
  79           MPFR_SET_SAME_SIGN (y, x);
  80           mpfr_set_divby0 ();
  81           MPFR_RET (0);
  82         }
  83       else /* y = x/1 = x */
  84         return mpfr_set (y, x, rnd_mode);
  85     }
  86   else if (MPFR_UNLIKELY (IS_POW2 (u)))
  87     return mpfr_div_2si (y, x, MPFR_INT_CEIL_LOG2 (u), rnd_mode);
  88
  89   MPFR_SET_SAME_SIGN (y, x);
  90
  91   MPFR_TMP_MARK (marker);
  92   xn = MPFR_LIMB_SIZE (x);
  93   yn = MPFR_LIMB_SIZE (y);
  94
  95   xp = MPFR_MANT (x);
  96   yp = MPFR_MANT (y);
  97   exp = MPFR_GET_EXP (x);
  98
  99   dif = yn + 1 - xn;
 100
 101   /* we need to store yn+1 = xn + dif limbs of the quotient */
 102   /* don't use tmp=yp since the mpn_lshift call below requires yp >= tmp+1 */
 103   tmp = MPFR_TMP_LIMBS_ALLOC (yn + 1);
 104
 105   c = (mp_limb_t) u;
 106   MPFR_ASSERTN (u == c);
 107   if (dif >= 0)
 108     c = mpn_divrem_1 (tmp, dif, xp, xn, c); /* used all the dividend */
 109   else /* dif < 0 i.e. xn > yn, don't use the (-dif) low limbs from x */
 110     c = mpn_divrem_1 (tmp, 0, xp - dif, yn + 1, c);
 111
 112   inexact = (c != 0);
 113
 114   /* First pass in estimating next bit of the quotient, in case of RNDN    *
 115    * In case we just have the right number of bits (postpone this ?),      *
 116    * we need to check whether the remainder is more or less than half      *
 117    * the divisor. The test must be performed with a subtraction, so as     *
 118    * to prevent carries.                                                   */
 119
 120   if (MPFR_LIKELY (rnd_mode == MPFR_RNDN))
 121     {
 122       if (c < (mp_limb_t) u - c) /* We have u > c */
 123         middle = -1;
 124       else if (c > (mp_limb_t) u - c)
 125         middle = 1;
 126       else
 127         middle = 0; /* exactly in the middle */
 128     }
 129
 130   /* If we believe that we are right in the middle or exact, we should check
 131      that we did not neglect any word of x (division large / 1 -> small). */
 132
 133   for (i=0; ((inexact == 0) || (middle == 0)) && (i < -dif); i++)
 134     if (xp[i])
 135       inexact = middle = 1; /* larger than middle */
 136
 137   /*
 138      If the high limb of the result is 0 (xp[xn-1] < u), remove it.
 139      Otherwise, compute the left shift to be performed to normalize.
 140      In the latter case, we discard some low bits computed. They
 141      contain information useful for the rounding, hence the updating
 142      of middle and inexact.
 143   */
 144
 145   if (tmp[yn] == 0)
 146     {
 147       MPN_COPY(yp, tmp, yn);
 148       exp -= GMP_NUMB_BITS;
 149     }
 150   else
 151     {
 152       int shlz;
 153
 154       count_leading_zeros (shlz, tmp[yn]);
 155
 156       /* shift left to normalize */
 157       if (MPFR_LIKELY (shlz != 0))
 158         {
 159           mp_limb_t w = tmp[0] << shlz;
 160
 161           mpn_lshift (yp, tmp + 1, yn, shlz);
 162           yp[0] += tmp[0] >> (GMP_NUMB_BITS - shlz);
 163
 164           if (w > (MPFR_LIMB_ONE << (GMP_NUMB_BITS - 1)))
 165             { middle = 1; }
 166           else if (w < (MPFR_LIMB_ONE << (GMP_NUMB_BITS - 1)))
 167             { middle = -1; }
 168           else
 169             { middle = (c != 0); }
 170
 171           inexact = inexact || (w != 0);
 172           exp -= shlz;
 173         }
 174       else
 175         { /* this happens only if u == 1 and xp[xn-1] >=
 176              1<<(GMP_NUMB_BITS-1). It might be better to handle the
 177              u == 1 case seperately ?
 178           */
 179
 180              MPN_COPY (yp, tmp + 1, yn);
 181         }
 182     }
 183
 184   MPFR_UNSIGNED_MINUS_MODULO (sh, MPFR_PREC (y));
 185   /* it remains sh bits in less significant limb of y */
 186
 187   d = *yp & MPFR_LIMB_MASK (sh);
 188   *yp ^= d; /* set to zero lowest sh bits */
 189
 190   MPFR_TMP_FREE (marker);
 191
 192   if (exp < __gmpfr_emin - 1)
 193     return mpfr_underflow (y, rnd_mode == MPFR_RNDN ? MPFR_RNDZ : rnd_mode,
 194                            MPFR_SIGN (y));
 195
 196   if (MPFR_UNLIKELY (d == 0 && inexact == 0))
 197     nexttoinf = 0;  /* result is exact */
 198   else
 199     {
 200       MPFR_UPDATE2_RND_MODE(rnd_mode, MPFR_SIGN (y));
 201       switch (rnd_mode)
 202         {
 203         case MPFR_RNDZ:
 204           inexact = - MPFR_INT_SIGN (y);  /* result is inexact */
 205           nexttoinf = 0;
 206           break;
 207
 208         case MPFR_RNDA:
 209           inexact = MPFR_INT_SIGN (y);
 210           nexttoinf = 1;
 211           break;
 212
 213         default: /* should be MPFR_RNDN */
 214           MPFR_ASSERTD (rnd_mode == MPFR_RNDN);
 215           /* We have one more significant bit in yn. */
 216           if (sh && d < (MPFR_LIMB_ONE << (sh - 1)))
 217             {
 218               inexact = - MPFR_INT_SIGN (y);
 219               nexttoinf = 0;
 220             }
 221           else if (sh && d > (MPFR_LIMB_ONE << (sh - 1)))
 222             {
 223               inexact = MPFR_INT_SIGN (y);
 224               nexttoinf = 1;
 225             }
 226           else /* sh = 0 or d = 1 << (sh-1) */
 227             {
 228               /* The first case is "false" even rounding (significant bits
 229                  indicate even rounding, but the result is inexact, so up) ;
 230                  The second case is the case where middle should be used to
 231                  decide the direction of rounding (no further bit computed) ;
 232                  The third is the true even rounding.
 233               */
 234               if ((sh && inexact) || (!sh && middle > 0) ||
 235                   (!inexact && *yp & (MPFR_LIMB_ONE << sh)))
 236                 {
 237                   inexact = MPFR_INT_SIGN (y);
 238                   nexttoinf = 1;
 239                 }
 240               else
 241                 {
 242                   inexact = - MPFR_INT_SIGN (y);
 243                   nexttoinf = 0;
 244                 }
 245             }
 246         }
 247     }
 248
 249   if (nexttoinf &&
 250       MPFR_UNLIKELY (mpn_add_1 (yp, yp, yn, MPFR_LIMB_ONE << sh)))
 251     {
 252       exp++;
 253       yp[yn-1] = MPFR_LIMB_HIGHBIT;
 254     }
 255
 256   /* Set the exponent. Warning! One may still have an underflow. */
 257   MPFR_EXP (y) = exp;
 258
 259   return mpfr_check_range (y, inexact, rnd_mode);
 260 }
 261
 262 int
 263 mpfr_div_si (mpfr_ptr y, mpfr_srcptr x, long int u, mpfr_rnd_t rnd_mode)
 264 {
 265   int res;
 266
 267   MPFR_LOG_FUNC
 268     (("x[%Pu]=%.*Rg u=%ld rnd=%d",
 269       mpfr_get_prec(x), mpfr_log_prec, x, u, rnd_mode),
 270      ("y[%Pu]=%.*Rg inexact=%d",
 271       mpfr_get_prec(y), mpfr_log_prec, y, res));
 272
 273   if (u >= 0)
 274     res = mpfr_div_ui (y, x, u, rnd_mode);
 275   else
 276     {
 277       res = -mpfr_div_ui (y, x, -u, MPFR_INVERT_RND (rnd_mode));
 278       MPFR_CHANGE_SIGN (y);
 279     }
 280   return res;
 281 }