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