contrib/mpfr/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 Free Software Foundation, Inc.
   4 Contributed by the Arenaire and Cacao 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 2.1 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.LIB.  If not, write to
  20 the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
  21 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, mp_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   mp_exp_t exp;
  35   int inexact, middle = 1, nexttoinf;
  36   MPFR_TMP_DECL(marker);
  37
  38   if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
  39     {
  40       if (MPFR_IS_NAN (x))
  41         {
  42           MPFR_SET_NAN (y);
  43           MPFR_RET_NAN;
  44         }
  45       else if (MPFR_IS_INF (x))
  46         {
  47           MPFR_SET_INF (y);
  48           MPFR_SET_SAME_SIGN (y, x);
  49           MPFR_RET (0);
  50         }
  51       else
  52         {
  53           MPFR_ASSERTD (MPFR_IS_ZERO(x));
  54           if (u == 0) /* 0/0 is NaN */
  55             {
  56               MPFR_SET_NAN(y);
  57               MPFR_RET_NAN;
  58             }
  59           else
  60             {
  61               MPFR_SET_ZERO(y);
  62               MPFR_SET_SAME_SIGN (y, x);
  63               MPFR_RET(0);
  64             }
  65         }
  66     }
  67   else if (MPFR_UNLIKELY (u <= 1))
  68     {
  69       if (u < 1)
  70         {
  71           /* x/0 is Inf since x != 0*/
  72           MPFR_SET_INF (y);
  73           MPFR_SET_SAME_SIGN (y, x);
  74           MPFR_RET (0);
  75         }
  76       else /* y = x/1 = x */
  77         return mpfr_set (y, x, rnd_mode);
  78     }
  79   else if (MPFR_UNLIKELY (IS_POW2 (u)))
  80     return mpfr_div_2si (y, x, MPFR_INT_CEIL_LOG2 (u), rnd_mode);
  81
  82   MPFR_CLEAR_FLAGS (y);
  83
  84   MPFR_SET_SAME_SIGN (y, x);
  85
  86   MPFR_TMP_MARK (marker);
  87   xn = MPFR_LIMB_SIZE (x);
  88   yn = MPFR_LIMB_SIZE (y);
  89
  90   xp = MPFR_MANT (x);
  91   yp = MPFR_MANT (y);
  92   exp = MPFR_GET_EXP (x);
  93
  94   dif = yn + 1 - xn;
  95
  96   /* we need to store yn+1 = xn + dif limbs of the quotient */
  97   /* don't use tmp=yp since the mpn_lshift call below requires yp >= tmp+1 */
  98   tmp = (mp_limb_t*) MPFR_TMP_ALLOC ((yn + 1) * BYTES_PER_MP_LIMB);
  99
 100   c = (mp_limb_t) u;
 101   MPFR_ASSERTN (u == c);
 102   if (dif >= 0)
 103     c = mpn_divrem_1 (tmp, dif, xp, xn, c); /* used all the dividend */
 104   else /* dif < 0 i.e. xn > yn, don't use the (-dif) low limbs from x */
 105     c = mpn_divrem_1 (tmp, 0, xp - dif, yn + 1, c);
 106
 107   inexact = (c != 0);
 108
 109   /* First pass in estimating next bit of the quotient, in case of RNDN    *
 110    * In case we just have the right number of bits (postpone this ?),      *
 111    * we need to check whether the remainder is more or less than half      *
 112    * the divisor. The test must be performed with a subtraction, so as     *
 113    * to prevent carries.                                                   */
 114
 115   if (MPFR_LIKELY (rnd_mode == GMP_RNDN))
 116     {
 117       if (c < (mp_limb_t) u - c) /* We have u > c */
 118         middle = -1;
 119       else if (c > (mp_limb_t) u - c)
 120         middle = 1;
 121       else
 122         middle = 0; /* exactly in the middle */
 123     }
 124
 125   /* If we believe that we are right in the middle or exact, we should check
 126      that we did not neglect any word of x (division large / 1 -> small). */
 127
 128   for (i=0; ((inexact == 0) || (middle == 0)) && (i < -dif); i++)
 129     if (xp[i])
 130       inexact = middle = 1; /* larger than middle */
 131
 132   /*
 133      If the high limb of the result is 0 (xp[xn-1] < u), remove it.
 134      Otherwise, compute the left shift to be performed to normalize.
 135      In the latter case, we discard some low bits computed. They
 136      contain information useful for the rounding, hence the updating
 137      of middle and inexact.
 138   */
 139
 140   if (tmp[yn] == 0)
 141     {
 142       MPN_COPY(yp, tmp, yn);
 143       exp -= BITS_PER_MP_LIMB;
 144     }
 145   else
 146     {
 147       int shlz;
 148
 149       count_leading_zeros (shlz, tmp[yn]);
 150
 151       /* shift left to normalize */
 152       if (MPFR_LIKELY (shlz != 0))
 153         {
 154           mp_limb_t w = tmp[0] << shlz;
 155
 156           mpn_lshift (yp, tmp + 1, yn, shlz);
 157           yp[0] += tmp[0] >> (BITS_PER_MP_LIMB - shlz);
 158
 159           if (w > (MPFR_LIMB_ONE << (BITS_PER_MP_LIMB - 1)))
 160             { middle = 1; }
 161           else if (w < (MPFR_LIMB_ONE << (BITS_PER_MP_LIMB - 1)))
 162             { middle = -1; }
 163           else
 164             { middle = (c != 0); }
 165
 166           inexact = inexact || (w != 0);
 167           exp -= shlz;
 168         }
 169       else
 170         { /* this happens only if u == 1 and xp[xn-1] >=
 171              1<<(BITS_PER_MP_LIMB-1). It might be better to handle the
 172              u == 1 case seperately ?
 173           */
 174
 175              MPN_COPY (yp, tmp + 1, yn);
 176         }
 177     }
 178
 179   MPFR_UNSIGNED_MINUS_MODULO (sh, MPFR_PREC (y));
 180   /* it remains sh bits in less significant limb of y */
 181
 182   d = *yp & MPFR_LIMB_MASK (sh);
 183   *yp ^= d; /* set to zero lowest sh bits */
 184
 185   MPFR_TMP_FREE (marker);
 186
 187   if (exp < __gmpfr_emin - 1)
 188     return mpfr_underflow (y, rnd_mode == GMP_RNDN ? GMP_RNDZ : rnd_mode,
 189                            MPFR_SIGN (y));
 190
 191   if (MPFR_UNLIKELY (d == 0 && inexact == 0))
 192     nexttoinf = 0;  /* result is exact */
 193   else
 194     switch (rnd_mode)
 195       {
 196       case GMP_RNDZ:
 197         inexact = - MPFR_INT_SIGN (y);  /* result is inexact */
 198         nexttoinf = 0;
 199         break;
 200
 201       case GMP_RNDU:
 202         inexact = 1;
 203         nexttoinf = MPFR_IS_POS (y);
 204         break;
 205
 206       case GMP_RNDD:
 207         inexact = -1;
 208         nexttoinf = MPFR_IS_NEG (y);
 209         break;
 210
 211       default:
 212         MPFR_ASSERTD (rnd_mode == GMP_RNDN);
 213         /* We have one more significant bit in yn. */
 214         if (sh && d < (MPFR_LIMB_ONE << (sh - 1)))
 215           {
 216             inexact = - MPFR_INT_SIGN (y);
 217             nexttoinf = 0;
 218           }
 219         else if (sh && d > (MPFR_LIMB_ONE << (sh - 1)))
 220           {
 221             inexact = MPFR_INT_SIGN (y);
 222             nexttoinf = 1;
 223           }
 224         else /* sh = 0 or d = 1 << (sh-1) */
 225           {
 226             /* The first case is "false" even rounding (significant bits
 227                indicate even rounding, but the result is inexact, so up) ;
 228                The second case is the case where middle should be used to
 229                decide the direction of rounding (no further bit computed) ;
 230                The third is the true even rounding.
 231             */
 232             if ((sh && inexact) || (!sh && middle > 0) ||
 233                 (!inexact && *yp & (MPFR_LIMB_ONE << sh)))
 234               {
 235                 inexact = MPFR_INT_SIGN (y);
 236                 nexttoinf = 1;
 237               }
 238             else
 239               {
 240                 inexact = - MPFR_INT_SIGN (y);
 241                 nexttoinf = 0;
 242               }
 243           }
 244       }
 245
 246   if (nexttoinf &&
 247       MPFR_UNLIKELY (mpn_add_1 (yp, yp, yn, MPFR_LIMB_ONE << sh)))
 248     {
 249       exp++;
 250       yp[yn-1] = MPFR_LIMB_HIGHBIT;
 251     }
 252
 253   /* Set the exponent. Warning! One may still have an underflow. */
 254   MPFR_EXP (y) = exp;
 255
 256   return mpfr_check_range (y, inexact, rnd_mode);
 257 }
 258
 259 int mpfr_div_si (mpfr_ptr y, mpfr_srcptr x, long int u, mp_rnd_t rnd_mode)
 260 {
 261   int res;
 262
 263   if (u >= 0)
 264     res = mpfr_div_ui (y, x, u, rnd_mode);
 265   else
 266     {
 267       res = -mpfr_div_ui (y, x, -u, MPFR_INVERT_RND (rnd_mode));
 268       MPFR_CHANGE_SIGN (y);
 269     }
 270   return res;
 271 }