1/* mpfr_div_{ui,si} -- divide a floating-point number by a machine integer 2 3Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc. 4Contributed by the AriC and Caramel projects, INRIA. 5 6This file is part of the GNU MPFR Library. 7 8The GNU MPFR Library is free software; you can redistribute it and/or modify 9it under the terms of the GNU Lesser General Public License as published by 10the Free Software Foundation; either version 3 of the License, or (at your 11option) any later version. 12 13The GNU MPFR Library is distributed in the hope that it will be useful, but 14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16License for more details. 17 18You should have received a copy of the GNU Lesser General Public License 19along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 20http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 2151 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 */ 27int 28mpfr_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 262int 263mpfr_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