1 /* mpfr_div_ui -- divide a floating-point number by a machine integer 2 3 Copyright 1999-2018 Free Software Foundation, Inc. 4 Contributed by the AriC and Caramba 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 #ifdef MPFR_COV_CHECK 27 int __gmpfr_cov_div_ui_sb[10][2] = { 0 }; 28 #endif 29 30 /* returns 0 if result exact, non-zero otherwise */ 31 #undef mpfr_div_ui 32 MPFR_HOT_FUNCTION_ATTR int 33 mpfr_div_ui (mpfr_ptr y, mpfr_srcptr x, unsigned long int u, 34 mpfr_rnd_t rnd_mode) 35 { 36 int sh; 37 mp_size_t i, xn, yn, dif; 38 mp_limb_t *xp, *yp, *tmp, c, d; 39 mpfr_exp_t exp; 40 int inexact; 41 mp_limb_t rb; /* round bit */ 42 mp_limb_t sb; /* sticky bit */ 43 MPFR_TMP_DECL(marker); 44 45 MPFR_LOG_FUNC 46 (("x[%Pu]=%.*Rg u=%lu rnd=%d", 47 mpfr_get_prec(x), mpfr_log_prec, x, u, rnd_mode), 48 ("y[%Pu]=%.*Rg inexact=%d", 49 mpfr_get_prec(y), mpfr_log_prec, y, inexact)); 50 51 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x))) 52 { 53 if (MPFR_IS_NAN (x)) 54 { 55 MPFR_SET_NAN (y); 56 MPFR_RET_NAN; 57 } 58 else if (MPFR_IS_INF (x)) 59 { 60 MPFR_SET_INF (y); 61 MPFR_SET_SAME_SIGN (y, x); 62 MPFR_RET (0); 63 } 64 else 65 { 66 MPFR_ASSERTD (MPFR_IS_ZERO (x)); 67 if (u == 0) /* 0/0 is NaN */ 68 { 69 MPFR_SET_NAN (y); 70 MPFR_RET_NAN; 71 } 72 else 73 { 74 MPFR_SET_ZERO(y); 75 MPFR_SET_SAME_SIGN (y, x); 76 MPFR_RET(0); 77 } 78 } 79 } 80 else if (MPFR_UNLIKELY (u <= 1)) 81 { 82 if (u < 1) 83 { 84 /* x/0 is Inf since x != 0 */ 85 MPFR_SET_INF (y); 86 MPFR_SET_SAME_SIGN (y, x); 87 MPFR_SET_DIVBY0 (); 88 MPFR_RET (0); 89 } 90 else /* y = x/1 = x */ 91 return mpfr_set (y, x, rnd_mode); 92 } 93 else if (MPFR_UNLIKELY (IS_POW2 (u))) 94 return mpfr_div_2si (y, x, MPFR_INT_CEIL_LOG2 (u), rnd_mode); 95 96 MPFR_SET_SAME_SIGN (y, x); 97 98 MPFR_TMP_MARK (marker); 99 100 xn = MPFR_LIMB_SIZE (x); 101 yn = MPFR_LIMB_SIZE (y); 102 103 xp = MPFR_MANT (x); 104 yp = MPFR_MANT (y); 105 exp = MPFR_GET_EXP (x); 106 107 dif = yn + 1 - xn; 108 109 /* we need to store yn + 1 = xn + dif limbs of the quotient */ 110 tmp = MPFR_TMP_LIMBS_ALLOC (yn + 1); 111 112 /* Notation: {p, n} denotes the integer formed by the n limbs 113 from p[0] to p[n-1]. Let B = 2^GMP_NUMB_BITS. 114 One has: 0 <= {p, n} < B^n. */ 115 116 MPFR_STAT_STATIC_ASSERT (MPFR_LIMB_MAX >= ULONG_MAX); 117 if (dif >= 0) 118 { 119 c = mpn_divrem_1 (tmp, dif, xp, xn, u); /* used all the dividend */ 120 /* {xp, xn} = ({tmp, xn+dif} * u + c) * B^(-dif) 121 = ({tmp, yn+1} * u + c) * B^(-dif) */ 122 } 123 else /* dif < 0, i.e. xn > yn+1; ignore the (-dif) low limbs from x */ 124 { 125 c = mpn_divrem_1 (tmp, 0, xp - dif, yn + 1, u); 126 /* {xp-dif, yn+1} = {tmp, yn+1} * u + c 127 thus 128 {xp, xn} = {xp, -dif} + {xp-dif, yn+1} * B^(-dif) 129 = {xp, -dif} + ({tmp, yn+1} * u + c) * B^(-dif) */ 130 } 131 132 /* Let r = {xp, -dif} / B^(-dif) if dif < 0, r = 0 otherwise; 0 <= r < 1. 133 Then {xp, xn} = ({tmp, yn+1} * u + c + r) * B^(-dif). 134 x / u = ({xp, xn} / u) * B^(-xn) * 2^exp 135 = ({tmp, yn+1} + (c + r) / u) * B^(-(yn+1)) * 2^exp 136 where 0 <= (c + r) / u < 1. */ 137 138 for (sb = 0, i = 0; sb == 0 && i < -dif; i++) 139 if (xp[i]) 140 sb = 1; 141 /* sb != 0 iff r != 0 */ 142 143 /* 144 If the highest limb of the result is 0 (xp[xn-1] < u), remove it. 145 Otherwise, compute the left shift to be performed to normalize. 146 In the latter case, we discard some low bits computed. They 147 contain information useful for the rounding, hence the updating 148 of middle and inexact. 149 */ 150 151 MPFR_UNSIGNED_MINUS_MODULO (sh, MPFR_PREC (y)); 152 /* sh: number of the trailing bits of y */ 153 154 if (tmp[yn] == 0) 155 { 156 MPN_COPY(yp, tmp, yn); 157 exp -= GMP_NUMB_BITS; 158 if (sh == 0) /* round bit is 1 iff (c + r) / u >= 1/2 */ 159 { 160 /* In this case tmp[yn]=0 and sh=0, the round bit is not in 161 {tmp,yn+1}. It is 1 iff 2*(c+r) - u >= 0. This means that in 162 some cases, we should look at the most significant bit of r. */ 163 if (c >= u - c) /* i.e. 2c >= u: round bit is always 1 */ 164 { 165 rb = 1; 166 /* The sticky bit is 1 unless 2c-u = 0 and r = 0. */ 167 sb |= 2 * c - u; 168 MPFR_COV_SET (div_ui_sb[0][!!sb]); 169 } 170 else /* 2*c < u */ 171 { 172 /* The round bit is 1 iff r >= 1/2 and 2*(c+1/2) = u. */ 173 rb = (c == u/2) && (dif < 0) && (xp[-dif-1] & MPFR_LIMB_HIGHBIT); 174 /* If rb is set, we need to recompute sb, since it might have 175 taken into account the msb of xp[-dif-1]. */ 176 if (rb) 177 { 178 sb = xp[-dif-1] << 1; /* discard the most significant bit */ 179 for (i = 0; sb == 0 && i < -dif-1; i++) 180 if (xp[i]) 181 sb = 1; 182 /* The dif < -1 case with sb = 0, i.e. [2][0], will 183 ensure that the body of the loop is covered. */ 184 MPFR_COV_SET (div_ui_sb[1 + (dif < -1)][!!sb]); 185 } 186 else 187 { 188 sb |= c; 189 MPFR_COV_SET (div_ui_sb[3][!!sb]); 190 } 191 } 192 } 193 else 194 { 195 /* round bit is in tmp[0] */ 196 rb = tmp[0] & (MPFR_LIMB_ONE << (sh - 1)); 197 sb |= (tmp[0] & MPFR_LIMB_MASK(sh - 1)) | c; 198 MPFR_COV_SET (div_ui_sb[4+!!rb][!!sb]); 199 } 200 } 201 else /* tmp[yn] != 0 */ 202 { 203 int shlz; 204 mp_limb_t w; 205 206 MPFR_ASSERTD (tmp[yn] != 0); 207 count_leading_zeros (shlz, tmp[yn]); 208 209 MPFR_ASSERTD (u >= 2); /* see special cases at the beginning */ 210 MPFR_ASSERTD (shlz > 0); /* since u >= 2 */ 211 212 /* shift left to normalize */ 213 w = tmp[0] << shlz; 214 mpn_lshift (yp, tmp + 1, yn, shlz); 215 yp[0] |= tmp[0] >> (GMP_NUMB_BITS - shlz); 216 /* now {yp, yn} is the approximate quotient, w is the next limb */ 217 218 if (sh == 0) /* round bit is upper bit from w */ 219 { 220 rb = w & MPFR_LIMB_HIGHBIT; 221 sb |= (w - rb) | c; 222 MPFR_COV_SET (div_ui_sb[6+!!rb][!!sb]); 223 } 224 else 225 { 226 rb = yp[0] & (MPFR_LIMB_ONE << (sh - 1)); 227 sb |= (yp[0] & MPFR_LIMB_MASK(sh - 1)) | w | c; 228 MPFR_COV_SET (div_ui_sb[8+!!rb][!!sb]); 229 } 230 231 exp -= shlz; 232 } 233 234 d = yp[0] & MPFR_LIMB_MASK (sh); 235 yp[0] ^= d; /* clear the lowest sh bits */ 236 237 MPFR_TMP_FREE (marker); 238 239 if (MPFR_UNLIKELY (exp < __gmpfr_emin - 1)) 240 return mpfr_underflow (y, rnd_mode == MPFR_RNDN ? MPFR_RNDZ : rnd_mode, 241 MPFR_SIGN (y)); 242 243 if (MPFR_UNLIKELY (rb == 0 && sb == 0)) 244 inexact = 0; /* result is exact */ 245 else 246 { 247 int nexttoinf; 248 249 MPFR_UPDATE2_RND_MODE(rnd_mode, MPFR_SIGN (y)); 250 switch (rnd_mode) 251 { 252 case MPFR_RNDZ: 253 case MPFR_RNDF: 254 inexact = - MPFR_INT_SIGN (y); /* result is inexact */ 255 nexttoinf = 0; 256 break; 257 258 case MPFR_RNDA: 259 inexact = MPFR_INT_SIGN (y); 260 nexttoinf = 1; 261 break; 262 263 default: /* should be MPFR_RNDN */ 264 MPFR_ASSERTD (rnd_mode == MPFR_RNDN); 265 /* We have one more significant bit in yn. */ 266 if (rb == 0) 267 { 268 inexact = - MPFR_INT_SIGN (y); 269 nexttoinf = 0; 270 } 271 else if (sb != 0) /* necessarily rb != 0 */ 272 { 273 inexact = MPFR_INT_SIGN (y); 274 nexttoinf = 1; 275 } 276 else /* middle case */ 277 { 278 if (yp[0] & (MPFR_LIMB_ONE << sh)) 279 { 280 inexact = MPFR_INT_SIGN (y); 281 nexttoinf = 1; 282 } 283 else 284 { 285 inexact = - MPFR_INT_SIGN (y); 286 nexttoinf = 0; 287 } 288 } 289 } 290 if (nexttoinf && 291 MPFR_UNLIKELY (mpn_add_1 (yp, yp, yn, MPFR_LIMB_ONE << sh))) 292 { 293 exp++; 294 yp[yn-1] = MPFR_LIMB_HIGHBIT; 295 } 296 } 297 298 /* Set the exponent. Warning! One may still have an underflow. */ 299 MPFR_EXP (y) = exp; 300 301 return mpfr_check_range (y, inexact, rnd_mode); 302 } 303