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 } 282