1 /* mpfr_mul -- multiply two floating-point numbers 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 27 /********* BEGINNING CHECK *************/ 28 29 /* Check if we have to check the result of mpfr_mul. 30 TODO: Find a better (and faster?) check than using old implementation */ 31 #if MPFR_WANT_ASSERT >= 2 32 33 int mpfr_mul2 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode); 34 static int 35 mpfr_mul3 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) 36 { 37 /* Old implementation */ 38 int sign_product, cc, inexact; 39 mpfr_exp_t ax; 40 mp_limb_t *tmp; 41 mp_limb_t b1; 42 mpfr_prec_t bq, cq; 43 mp_size_t bn, cn, tn, k; 44 MPFR_TMP_DECL(marker); 45 46 /* deal with special cases */ 47 if (MPFR_ARE_SINGULAR(b,c)) 48 { 49 if (MPFR_IS_NAN(b) || MPFR_IS_NAN(c)) 50 { 51 MPFR_SET_NAN(a); 52 MPFR_RET_NAN; 53 } 54 sign_product = MPFR_MULT_SIGN( MPFR_SIGN(b) , MPFR_SIGN(c) ); 55 if (MPFR_IS_INF(b)) 56 { 57 if (MPFR_IS_INF(c) || MPFR_NOTZERO(c)) 58 { 59 MPFR_SET_SIGN(a, sign_product); 60 MPFR_SET_INF(a); 61 MPFR_RET(0); /* exact */ 62 } 63 else 64 { 65 MPFR_SET_NAN(a); 66 MPFR_RET_NAN; 67 } 68 } 69 else if (MPFR_IS_INF(c)) 70 { 71 if (MPFR_NOTZERO(b)) 72 { 73 MPFR_SET_SIGN(a, sign_product); 74 MPFR_SET_INF(a); 75 MPFR_RET(0); /* exact */ 76 } 77 else 78 { 79 MPFR_SET_NAN(a); 80 MPFR_RET_NAN; 81 } 82 } 83 else 84 { 85 MPFR_ASSERTD(MPFR_IS_ZERO(b) || MPFR_IS_ZERO(c)); 86 MPFR_SET_SIGN(a, sign_product); 87 MPFR_SET_ZERO(a); 88 MPFR_RET(0); /* 0 * 0 is exact */ 89 } 90 } 91 sign_product = MPFR_MULT_SIGN( MPFR_SIGN(b) , MPFR_SIGN(c) ); 92 93 ax = MPFR_GET_EXP (b) + MPFR_GET_EXP (c); 94 95 bq = MPFR_PREC (b); 96 cq = MPFR_PREC (c); 97 98 MPFR_ASSERTN ((mpfr_uprec_t) bq + cq <= MPFR_PREC_MAX); 99 100 bn = MPFR_PREC2LIMBS (bq); /* number of limbs of b */ 101 cn = MPFR_PREC2LIMBS (cq); /* number of limbs of c */ 102 k = bn + cn; /* effective nb of limbs used by b*c (= tn or tn+1) below */ 103 tn = MPFR_PREC2LIMBS (bq + cq); 104 /* <= k, thus no int overflow */ 105 MPFR_ASSERTD(tn <= k); 106 107 /* Check for no size_t overflow*/ 108 MPFR_ASSERTD((size_t) k <= ((size_t) -1) / MPFR_BYTES_PER_MP_LIMB); 109 MPFR_TMP_MARK(marker); 110 tmp = MPFR_TMP_LIMBS_ALLOC (k); 111 112 /* multiplies two mantissa in temporary allocated space */ 113 b1 = (MPFR_LIKELY(bn >= cn)) ? 114 mpn_mul (tmp, MPFR_MANT(b), bn, MPFR_MANT(c), cn) 115 : mpn_mul (tmp, MPFR_MANT(c), cn, MPFR_MANT(b), bn); 116 117 /* now tmp[0]..tmp[k-1] contains the product of both mantissa, 118 with tmp[k-1]>=2^(GMP_NUMB_BITS-2) */ 119 b1 >>= GMP_NUMB_BITS - 1; /* msb from the product */ 120 121 /* if the mantissas of b and c are uniformly distributed in ]1/2, 1], 122 then their product is in ]1/4, 1/2] with probability 2*ln(2)-1 ~ 0.386 123 and in [1/2, 1] with probability 2-2*ln(2) ~ 0.614 */ 124 tmp += k - tn; 125 if (MPFR_UNLIKELY(b1 == 0)) 126 mpn_lshift (tmp, tmp, tn, 1); /* tn <= k, so no stack corruption */ 127 cc = mpfr_round_raw (MPFR_MANT (a), tmp, bq + cq, 128 MPFR_IS_NEG_SIGN(sign_product), 129 MPFR_PREC (a), rnd_mode, &inexact); 130 131 /* cc = 1 ==> result is a power of two */ 132 if (MPFR_UNLIKELY(cc)) 133 MPFR_MANT(a)[MPFR_LIMB_SIZE(a)-1] = MPFR_LIMB_HIGHBIT; 134 135 MPFR_TMP_FREE(marker); 136 137 { 138 mpfr_exp_t ax2 = ax + (mpfr_exp_t) (b1 - 1 + cc); 139 if (MPFR_UNLIKELY( ax2 > __gmpfr_emax)) 140 return mpfr_overflow (a, rnd_mode, sign_product); 141 if (MPFR_UNLIKELY( ax2 < __gmpfr_emin)) 142 { 143 /* In the rounding to the nearest mode, if the exponent of the exact 144 result (i.e. before rounding, i.e. without taking cc into account) 145 is < __gmpfr_emin - 1 or the exact result is a power of 2 (i.e. if 146 both arguments are powers of 2) in absolute value, then round to 147 zero. */ 148 if (rnd_mode == MPFR_RNDN && 149 (ax + (mpfr_exp_t) b1 < __gmpfr_emin || 150 (mpfr_powerof2_raw (b) && mpfr_powerof2_raw (c)))) 151 rnd_mode = MPFR_RNDZ; 152 return mpfr_underflow (a, rnd_mode, sign_product); 153 } 154 MPFR_SET_EXP (a, ax2); 155 MPFR_SET_SIGN(a, sign_product); 156 } 157 MPFR_RET (inexact); 158 } 159 160 int 161 mpfr_mul (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) 162 { 163 mpfr_t ta, tb, tc; 164 int inexact1, inexact2; 165 166 if (rnd_mode == MPFR_RNDF) 167 return mpfr_mul2 (a, b, c, rnd_mode); 168 169 mpfr_init2 (ta, MPFR_PREC (a)); 170 mpfr_init2 (tb, MPFR_PREC (b)); 171 mpfr_init2 (tc, MPFR_PREC (c)); 172 MPFR_ASSERTN (mpfr_set (tb, b, MPFR_RNDN) == 0); 173 MPFR_ASSERTN (mpfr_set (tc, c, MPFR_RNDN) == 0); 174 175 inexact2 = mpfr_mul3 (ta, tb, tc, rnd_mode); 176 inexact1 = mpfr_mul2 (a, b, c, rnd_mode); 177 if (MPFR_IS_NAN (ta) && MPFR_IS_NAN (a)) 178 { 179 /* Getting both NaN is OK. */ 180 } 181 else if (! mpfr_equal_p (ta, a) || ! SAME_SIGN (inexact1, inexact2)) 182 { 183 fprintf (stderr, "mpfr_mul return different values for %s\n" 184 "Prec_a = %lu, Prec_b = %lu, Prec_c = %lu\nb = ", 185 mpfr_print_rnd_mode (rnd_mode), 186 MPFR_PREC (a), MPFR_PREC (b), MPFR_PREC (c)); 187 mpfr_fdump (stderr, b); 188 fprintf (stderr, "c = "); 189 mpfr_fdump (stderr, c); 190 fprintf (stderr, "OldMul: "); 191 mpfr_fdump (stderr, ta); 192 fprintf (stderr, "NewMul: "); 193 mpfr_fdump (stderr, a); 194 fprintf (stderr, "NewInexact = %d | OldInexact = %d\n", 195 inexact1, inexact2); 196 MPFR_ASSERTN(0); 197 } 198 199 mpfr_clears (ta, tb, tc, (mpfr_ptr) 0); 200 return inexact1; 201 } 202 203 # define mpfr_mul mpfr_mul2 204 205 #endif /* MPFR_WANT_ASSERT >= 2 */ 206 207 /****** END OF CHECK *******/ 208 209 /* Multiply 2 mpfr_t */ 210 211 #if !defined(MPFR_GENERIC_ABI) 212 213 /* Special code for prec(a) < GMP_NUMB_BITS and 214 prec(b), prec(c) <= GMP_NUMB_BITS. 215 Note: this code was copied in sqr.c, function mpfr_sqr_1 (this saves a few cycles 216 with respect to have this function exported). As a consequence, any change here 217 should be reported in mpfr_sqr_1. */ 218 static int 219 mpfr_mul_1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode, 220 mpfr_prec_t p) 221 { 222 mp_limb_t a0; 223 mpfr_limb_ptr ap = MPFR_MANT(a); 224 mp_limb_t b0 = MPFR_MANT(b)[0]; 225 mp_limb_t c0 = MPFR_MANT(c)[0]; 226 mpfr_exp_t ax; 227 mpfr_prec_t sh = GMP_NUMB_BITS - p; 228 mp_limb_t rb, sb, mask = MPFR_LIMB_MASK(sh); 229 230 /* When prec(b), prec(c) <= GMP_NUMB_BITS / 2, we could replace umul_ppmm 231 by a limb multiplication as follows, but we assume umul_ppmm is as fast 232 as a limb multiplication on modern processors: 233 a0 = (b0 >> (GMP_NUMB_BITS / 2)) * (c0 >> (GMP_NUMB_BITS / 2)); 234 sb = 0; 235 */ 236 ax = MPFR_GET_EXP(b) + MPFR_GET_EXP(c); 237 umul_ppmm (a0, sb, b0, c0); 238 if (a0 < MPFR_LIMB_HIGHBIT) 239 { 240 ax --; 241 /* TODO: This is actually an addition with carry (no shifts and no OR 242 needed in asm). Make sure that GCC generates optimized code once 243 it supports carry-in. */ 244 a0 = (a0 << 1) | (sb >> (GMP_NUMB_BITS - 1)); 245 sb <<= 1; 246 } 247 rb = a0 & (MPFR_LIMB_ONE << (sh - 1)); 248 sb |= (a0 & mask) ^ rb; 249 ap[0] = a0 & ~mask; 250 251 MPFR_SIGN(a) = MPFR_MULT_SIGN (MPFR_SIGN (b), MPFR_SIGN (c)); 252 253 /* rounding */ 254 if (MPFR_UNLIKELY(ax > __gmpfr_emax)) 255 return mpfr_overflow (a, rnd_mode, MPFR_SIGN(a)); 256 257 /* Warning: underflow should be checked *after* rounding, thus when rounding 258 away and when a > 0.111...111*2^(emin-1), or when rounding to nearest and 259 a >= 0.111...111[1]*2^(emin-1), there is no underflow. */ 260 if (MPFR_UNLIKELY(ax < __gmpfr_emin)) 261 { 262 if (ax == __gmpfr_emin - 1 && ap[0] == ~mask && 263 ((rnd_mode == MPFR_RNDN && rb) || 264 (MPFR_IS_LIKE_RNDA(rnd_mode, MPFR_IS_NEG (a)) && (rb | sb)))) 265 goto rounding; /* no underflow */ 266 /* For RNDN, mpfr_underflow always rounds away, thus for |a| <= 2^(emin-2) 267 we have to change to RNDZ. This corresponds to: 268 (a) either ax < emin - 1 269 (b) or ax = emin - 1 and ap[0] = 1000....000 and rb = sb = 0 */ 270 if (rnd_mode == MPFR_RNDN && 271 (ax < __gmpfr_emin - 1 || 272 (ap[0] == MPFR_LIMB_HIGHBIT && (rb | sb) == 0))) 273 rnd_mode = MPFR_RNDZ; 274 return mpfr_underflow (a, rnd_mode, MPFR_SIGN(a)); 275 } 276 277 rounding: 278 MPFR_EXP (a) = ax; /* Don't use MPFR_SET_EXP since ax might be < __gmpfr_emin 279 in the cases "goto rounding" above. */ 280 if ((rb == 0 && sb == 0) || rnd_mode == MPFR_RNDF) 281 { 282 MPFR_ASSERTD(ax >= __gmpfr_emin); 283 MPFR_RET (0); 284 } 285 else if (rnd_mode == MPFR_RNDN) 286 { 287 if (rb == 0 || (sb == 0 && (ap[0] & (MPFR_LIMB_ONE << sh)) == 0)) 288 goto truncate; 289 else 290 goto add_one_ulp; 291 } 292 else if (MPFR_IS_LIKE_RNDZ(rnd_mode, MPFR_IS_NEG(a))) 293 { 294 truncate: 295 MPFR_ASSERTD(ax >= __gmpfr_emin); 296 MPFR_RET(-MPFR_SIGN(a)); 297 } 298 else /* round away from zero */ 299 { 300 add_one_ulp: 301 ap[0] += MPFR_LIMB_ONE << sh; 302 if (ap[0] == 0) 303 { 304 ap[0] = MPFR_LIMB_HIGHBIT; 305 if (MPFR_UNLIKELY(ax + 1 > __gmpfr_emax)) 306 return mpfr_overflow (a, rnd_mode, MPFR_SIGN(a)); 307 MPFR_ASSERTD(ax + 1 <= __gmpfr_emax); 308 MPFR_ASSERTD(ax + 1 >= __gmpfr_emin); 309 MPFR_SET_EXP (a, ax + 1); 310 } 311 MPFR_RET(MPFR_SIGN(a)); 312 } 313 } 314 315 /* Special code for prec(a) = GMP_NUMB_BITS and 316 prec(b), prec(c) <= GMP_NUMB_BITS. */ 317 static int 318 mpfr_mul_1n (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) 319 { 320 mp_limb_t a0; 321 mpfr_limb_ptr ap = MPFR_MANT(a); 322 mp_limb_t b0 = MPFR_MANT(b)[0]; 323 mp_limb_t c0 = MPFR_MANT(c)[0]; 324 mpfr_exp_t ax; 325 mp_limb_t rb, sb; 326 327 ax = MPFR_GET_EXP(b) + MPFR_GET_EXP(c); 328 umul_ppmm (a0, sb, b0, c0); 329 if (a0 < MPFR_LIMB_HIGHBIT) 330 { 331 ax --; 332 /* TODO: This is actually an addition with carry (no shifts and no OR 333 needed in asm). Make sure that GCC generates optimized code once 334 it supports carry-in. */ 335 a0 = (a0 << 1) | (sb >> (GMP_NUMB_BITS - 1)); 336 sb <<= 1; 337 } 338 rb = sb & MPFR_LIMB_HIGHBIT; 339 sb = sb & ~MPFR_LIMB_HIGHBIT; 340 ap[0] = a0; 341 342 MPFR_SIGN(a) = MPFR_MULT_SIGN (MPFR_SIGN (b), MPFR_SIGN (c)); 343 344 /* rounding */ 345 if (MPFR_UNLIKELY(ax > __gmpfr_emax)) 346 return mpfr_overflow (a, rnd_mode, MPFR_SIGN(a)); 347 348 /* Warning: underflow should be checked *after* rounding, thus when rounding 349 away and when a > 0.111...111*2^(emin-1), or when rounding to nearest and 350 a >= 0.111...111[1]*2^(emin-1), there is no underflow. 351 Note: this case can only occur when the initial a0 (after the umul_ppmm 352 call above) had its most significant bit 0, since the largest a0 is 353 obtained for b0 = c0 = B-1 where B=2^GMP_NUMB_BITS, thus b0*c0 <= (B-1)^2 354 thus a0 <= B-2. */ 355 if (MPFR_UNLIKELY(ax < __gmpfr_emin)) 356 { 357 if (ax == __gmpfr_emin - 1 && ap[0] == ~MPFR_LIMB_ZERO && 358 ((rnd_mode == MPFR_RNDN && rb) || 359 (MPFR_IS_LIKE_RNDA(rnd_mode, MPFR_IS_NEG (a)) && (rb | sb)))) 360 goto rounding; /* no underflow */ 361 /* For RNDN, mpfr_underflow always rounds away, thus for |a| <= 2^(emin-2) 362 we have to change to RNDZ. This corresponds to: 363 (a) either ax < emin - 1 364 (b) or ax = emin - 1 and ap[0] = 1000....000 and rb = sb = 0 */ 365 if (rnd_mode == MPFR_RNDN && 366 (ax < __gmpfr_emin - 1 || 367 (ap[0] == MPFR_LIMB_HIGHBIT && (rb | sb) == 0))) 368 rnd_mode = MPFR_RNDZ; 369 return mpfr_underflow (a, rnd_mode, MPFR_SIGN(a)); 370 } 371 372 rounding: 373 MPFR_EXP (a) = ax; /* Don't use MPFR_SET_EXP since ax might be < __gmpfr_emin 374 in the cases "goto rounding" above. */ 375 if ((rb == 0 && sb == 0) || rnd_mode == MPFR_RNDF) 376 { 377 MPFR_ASSERTD(ax >= __gmpfr_emin); 378 MPFR_RET (0); 379 } 380 else if (rnd_mode == MPFR_RNDN) 381 { 382 if (rb == 0 || (sb == 0 && (ap[0] & MPFR_LIMB_ONE) == 0)) 383 goto truncate; 384 else 385 goto add_one_ulp; 386 } 387 else if (MPFR_IS_LIKE_RNDZ(rnd_mode, MPFR_IS_NEG(a))) 388 { 389 truncate: 390 MPFR_ASSERTD(ax >= __gmpfr_emin); 391 MPFR_RET(-MPFR_SIGN(a)); 392 } 393 else /* round away from zero */ 394 { 395 add_one_ulp: 396 ap[0] += MPFR_LIMB_ONE; 397 if (ap[0] == 0) 398 { 399 ap[0] = MPFR_LIMB_HIGHBIT; 400 if (MPFR_UNLIKELY(ax + 1 > __gmpfr_emax)) 401 return mpfr_overflow (a, rnd_mode, MPFR_SIGN(a)); 402 MPFR_ASSERTD(ax + 1 <= __gmpfr_emax); 403 MPFR_ASSERTD(ax + 1 >= __gmpfr_emin); 404 MPFR_SET_EXP (a, ax + 1); 405 } 406 MPFR_RET(MPFR_SIGN(a)); 407 } 408 } 409 410 /* Special code for GMP_NUMB_BITS < prec(a) < 2*GMP_NUMB_BITS and 411 GMP_NUMB_BITS < prec(b), prec(c) <= 2*GMP_NUMB_BITS. 412 Note: this code was copied in sqr.c, function mpfr_sqr_2 (this saves a few cycles 413 with respect to have this function exported). As a consequence, any change here 414 should be reported in mpfr_sqr_2. */ 415 static int 416 mpfr_mul_2 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode, 417 mpfr_prec_t p) 418 { 419 mp_limb_t h, l, u, v, w; 420 mpfr_limb_ptr ap = MPFR_MANT(a); 421 mpfr_exp_t ax = MPFR_GET_EXP(b) + MPFR_GET_EXP(c); 422 mpfr_prec_t sh = 2 * GMP_NUMB_BITS - p; 423 mp_limb_t rb, sb, sb2, mask = MPFR_LIMB_MASK(sh); 424 mp_limb_t *bp = MPFR_MANT(b), *cp = MPFR_MANT(c); 425 426 /* we store the 4-limb product in h=ap[1], l=ap[0], sb=ap[-1], sb2=ap[-2] */ 427 umul_ppmm (h, l, bp[1], cp[1]); 428 umul_ppmm (u, v, bp[1], cp[0]); 429 l += u; 430 h += (l < u); 431 umul_ppmm (u, w, bp[0], cp[1]); 432 l += u; 433 h += (l < u); 434 435 /* now the full product is {h, l, v + w + high(b0*c0), low(b0*c0)}, 436 where the lower part contributes to less than 3 ulps to {h, l} */ 437 438 /* If h has its most significant bit set and the low sh-1 bits of l are not 439 000...000 nor 111...111 nor 111...110, then we can round correctly; 440 if h has zero as most significant bit, we have to shift left h and l, 441 thus if the low sh-2 bits are not 000...000 nor 111...111 nor 111...110, 442 then we can round correctly. To avoid an extra test we consider the latter 443 case (if we can round, we can also round in the former case). 444 For sh <= 3, we have mask <= 7, thus (mask>>2) <= 1, and the approximation 445 cannot be enough. */ 446 if (MPFR_LIKELY(((l + 2) & (mask >> 2)) > 2)) 447 sb = sb2 = 1; /* result cannot be exact in that case */ 448 else 449 { 450 umul_ppmm (sb, sb2, bp[0], cp[0]); 451 /* the full product is {h, l, sb + v + w, sb2} */ 452 sb += v; 453 l += (sb < v); 454 h += (l == 0) && (sb < v); 455 sb += w; 456 l += (sb < w); 457 h += (l == 0) && (sb < w); 458 } 459 if (h < MPFR_LIMB_HIGHBIT) 460 { 461 ax --; 462 h = (h << 1) | (l >> (GMP_NUMB_BITS - 1)); 463 l = (l << 1) | (sb >> (GMP_NUMB_BITS - 1)); 464 sb <<= 1; 465 /* no need to shift sb2 since we only want to know if it is zero or not */ 466 } 467 ap[1] = h; 468 rb = l & (MPFR_LIMB_ONE << (sh - 1)); 469 sb |= ((l & mask) ^ rb) | sb2; 470 ap[0] = l & ~mask; 471 472 MPFR_SIGN(a) = MPFR_MULT_SIGN (MPFR_SIGN (b), MPFR_SIGN (c)); 473 474 /* rounding */ 475 if (MPFR_UNLIKELY(ax > __gmpfr_emax)) 476 return mpfr_overflow (a, rnd_mode, MPFR_SIGN(a)); 477 478 /* Warning: underflow should be checked *after* rounding, thus when rounding 479 away and when a > 0.111...111*2^(emin-1), or when rounding to nearest and 480 a >= 0.111...111[1]*2^(emin-1), there is no underflow. */ 481 if (MPFR_UNLIKELY(ax < __gmpfr_emin)) 482 { 483 if (ax == __gmpfr_emin - 1 && 484 ap[1] == MPFR_LIMB_MAX && 485 ap[0] == ~mask && 486 ((rnd_mode == MPFR_RNDN && rb) || 487 (MPFR_IS_LIKE_RNDA(rnd_mode, MPFR_IS_NEG (a)) && (rb | sb)))) 488 goto rounding; /* no underflow */ 489 /* for RNDN, mpfr_underflow always rounds away, thus for |a| <= 2^(emin-2) 490 we have to change to RNDZ */ 491 if (rnd_mode == MPFR_RNDN && 492 (ax < __gmpfr_emin - 1 || 493 (ap[1] == MPFR_LIMB_HIGHBIT && ap[0] == 0 && (rb | sb) == 0))) 494 rnd_mode = MPFR_RNDZ; 495 return mpfr_underflow (a, rnd_mode, MPFR_SIGN(a)); 496 } 497 498 rounding: 499 MPFR_EXP (a) = ax; /* Don't use MPFR_SET_EXP since ax might be < __gmpfr_emin 500 in the cases "goto rounding" above. */ 501 if ((rb == 0 && sb == 0) || rnd_mode == MPFR_RNDF) 502 { 503 MPFR_ASSERTD(ax >= __gmpfr_emin); 504 MPFR_RET (0); 505 } 506 else if (rnd_mode == MPFR_RNDN) 507 { 508 if (rb == 0 || (sb == 0 && (ap[0] & (MPFR_LIMB_ONE << sh)) == 0)) 509 goto truncate; 510 else 511 goto add_one_ulp; 512 } 513 else if (MPFR_IS_LIKE_RNDZ(rnd_mode, MPFR_IS_NEG(a))) 514 { 515 truncate: 516 MPFR_ASSERTD(ax >= __gmpfr_emin); 517 MPFR_RET(-MPFR_SIGN(a)); 518 } 519 else /* round away from zero */ 520 { 521 add_one_ulp: 522 ap[0] += MPFR_LIMB_ONE << sh; 523 ap[1] += (ap[0] == 0); 524 if (ap[1] == 0) 525 { 526 ap[1] = MPFR_LIMB_HIGHBIT; 527 if (MPFR_UNLIKELY(ax + 1 > __gmpfr_emax)) 528 return mpfr_overflow (a, rnd_mode, MPFR_SIGN(a)); 529 MPFR_ASSERTD(ax + 1 <= __gmpfr_emax); 530 MPFR_ASSERTD(ax + 1 >= __gmpfr_emin); 531 MPFR_SET_EXP (a, ax + 1); 532 } 533 MPFR_RET(MPFR_SIGN(a)); 534 } 535 } 536 537 /* Special code for 2*GMP_NUMB_BITS < prec(a) < 3*GMP_NUMB_BITS and 538 2*GMP_NUMB_BITS < prec(b), prec(c) <= 3*GMP_NUMB_BITS. */ 539 static int 540 mpfr_mul_3 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode, 541 mpfr_prec_t p) 542 { 543 mp_limb_t a0, a1, a2, h, l, cy; 544 mpfr_limb_ptr ap = MPFR_MANT(a); 545 mpfr_exp_t ax = MPFR_GET_EXP(b) + MPFR_GET_EXP(c); 546 mpfr_prec_t sh = 3 * GMP_NUMB_BITS - p; 547 mp_limb_t rb, sb, sb2, mask = MPFR_LIMB_MASK(sh); 548 mp_limb_t *bp = MPFR_MANT(b), *cp = MPFR_MANT(c); 549 550 /* we store the upper 3-limb product in a2, a1, a0: 551 b2*c2, b2*c1+b1*c2, b2*c0+b1*c1+b0*c2 */ 552 umul_ppmm (a2, a1, bp[2], cp[2]); 553 umul_ppmm (h, a0, bp[2], cp[1]); 554 a1 += h; 555 a2 += (a1 < h); 556 umul_ppmm (h, l, bp[1], cp[2]); 557 a1 += h; 558 a2 += (a1 < h); 559 a0 += l; 560 cy = a0 < l; /* carry in a1 */ 561 umul_ppmm (h, l, bp[2], cp[0]); 562 a0 += h; 563 cy += (a0 < h); 564 umul_ppmm (h, l, bp[1], cp[1]); 565 a0 += h; 566 cy += (a0 < h); 567 umul_ppmm (h, l, bp[0], cp[2]); 568 a0 += h; 569 cy += (a0 < h); 570 /* now propagate cy */ 571 a1 += cy; 572 a2 += (a1 < cy); 573 574 /* Now the approximate product {a2, a1, a0} has an error of less than 575 5 ulps (3 ulps for the ignored low limbs of b2*c0+b1*c1+b0*c2, 576 plus 2 ulps for the ignored b1*c0+b0*c1 (plus b0*c0)). 577 Since we might shift by 1 bit, we make sure the low sh-2 bits of a0 578 are not 0, -1, -2, -3 or -4. */ 579 580 if (MPFR_LIKELY(((a0 + 4) & (mask >> 2)) > 4)) 581 sb = sb2 = 1; /* result cannot be exact in that case */ 582 else 583 { 584 mp_limb_t p[6]; 585 mpn_mul_n (p, bp, cp, 3); 586 a2 = p[5]; 587 a1 = p[4]; 588 a0 = p[3]; 589 sb = p[2]; 590 sb2 = p[1] | p[0]; 591 } 592 if (a2 < MPFR_LIMB_HIGHBIT) 593 { 594 ax --; 595 a2 = (a2 << 1) | (a1 >> (GMP_NUMB_BITS - 1)); 596 a1 = (a1 << 1) | (a0 >> (GMP_NUMB_BITS - 1)); 597 a0 = (a0 << 1) | (sb >> (GMP_NUMB_BITS - 1)); 598 sb <<= 1; 599 /* no need to shift sb2: we only need to know if it is zero or not */ 600 } 601 ap[2] = a2; 602 ap[1] = a1; 603 rb = a0 & (MPFR_LIMB_ONE << (sh - 1)); 604 sb |= ((a0 & mask) ^ rb) | sb2; 605 ap[0] = a0 & ~mask; 606 607 MPFR_SIGN(a) = MPFR_MULT_SIGN (MPFR_SIGN (b), MPFR_SIGN (c)); 608 609 /* rounding */ 610 if (MPFR_UNLIKELY(ax > __gmpfr_emax)) 611 return mpfr_overflow (a, rnd_mode, MPFR_SIGN(a)); 612 613 /* Warning: underflow should be checked *after* rounding, thus when rounding 614 away and when a > 0.111...111*2^(emin-1), or when rounding to nearest and 615 a >= 0.111...111[1]*2^(emin-1), there is no underflow. */ 616 if (MPFR_UNLIKELY(ax < __gmpfr_emin)) 617 { 618 if (ax == __gmpfr_emin - 1 && 619 ap[2] == MPFR_LIMB_MAX && 620 ap[1] == MPFR_LIMB_MAX && 621 ap[0] == ~mask && 622 ((rnd_mode == MPFR_RNDN && rb) || 623 (MPFR_IS_LIKE_RNDA(rnd_mode, MPFR_IS_NEG (a)) && (rb | sb)))) 624 goto rounding; /* no underflow */ 625 /* for RNDN, mpfr_underflow always rounds away, thus for |a| <= 2^(emin-2) 626 we have to change to RNDZ */ 627 if (rnd_mode == MPFR_RNDN && 628 (ax < __gmpfr_emin - 1 || 629 (ap[2] == MPFR_LIMB_HIGHBIT && ap[1] == 0 && ap[0] == 0 630 && (rb | sb) == 0))) 631 rnd_mode = MPFR_RNDZ; 632 return mpfr_underflow (a, rnd_mode, MPFR_SIGN(a)); 633 } 634 635 rounding: 636 MPFR_EXP (a) = ax; /* Don't use MPFR_SET_EXP since ax might be < __gmpfr_emin 637 in the cases "goto rounding" above. */ 638 if ((rb == 0 && sb == 0) || rnd_mode == MPFR_RNDF) 639 { 640 MPFR_ASSERTD(ax >= __gmpfr_emin); 641 MPFR_RET (0); 642 } 643 else if (rnd_mode == MPFR_RNDN) 644 { 645 if (rb == 0 || (sb == 0 && (ap[0] & (MPFR_LIMB_ONE << sh)) == 0)) 646 goto truncate; 647 else 648 goto add_one_ulp; 649 } 650 else if (MPFR_IS_LIKE_RNDZ(rnd_mode, MPFR_IS_NEG(a))) 651 { 652 truncate: 653 MPFR_ASSERTD(ax >= __gmpfr_emin); 654 MPFR_RET(-MPFR_SIGN(a)); 655 } 656 else /* round away from zero */ 657 { 658 add_one_ulp: 659 ap[0] += MPFR_LIMB_ONE << sh; 660 ap[1] += (ap[0] == 0); 661 ap[2] += (ap[1] == 0) && (ap[0] == 0); 662 if (ap[2] == 0) 663 { 664 ap[2] = MPFR_LIMB_HIGHBIT; 665 if (MPFR_UNLIKELY(ax + 1 > __gmpfr_emax)) 666 return mpfr_overflow (a, rnd_mode, MPFR_SIGN(a)); 667 MPFR_ASSERTD(ax + 1 <= __gmpfr_emax); 668 MPFR_ASSERTD(ax + 1 >= __gmpfr_emin); 669 MPFR_SET_EXP (a, ax + 1); 670 } 671 MPFR_RET(MPFR_SIGN(a)); 672 } 673 } 674 675 #endif /* !defined(MPFR_GENERIC_ABI) */ 676 677 /* Note: mpfr_sqr will call mpfr_mul if bn > MPFR_SQR_THRESHOLD, 678 in order to use Mulders' mulhigh, which is handled only here 679 to avoid partial code duplication. There is some overhead due 680 to the additional tests, but slowdown should not be noticeable 681 as this code is not executed in very small precisions. */ 682 683 MPFR_HOT_FUNCTION_ATTR int 684 mpfr_mul (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode) 685 { 686 int sign, inexact; 687 mpfr_exp_t ax, ax2; 688 mp_limb_t *tmp; 689 mp_limb_t b1; 690 mpfr_prec_t aq, bq, cq; 691 mp_size_t bn, cn, tn, k, threshold; 692 MPFR_TMP_DECL (marker); 693 694 MPFR_LOG_FUNC 695 (("b[%Pu]=%.*Rg c[%Pu]=%.*Rg rnd=%d", 696 mpfr_get_prec (b), mpfr_log_prec, b, 697 mpfr_get_prec (c), mpfr_log_prec, c, rnd_mode), 698 ("a[%Pu]=%.*Rg inexact=%d", 699 mpfr_get_prec (a), mpfr_log_prec, a, inexact)); 700 701 /* deal with special cases */ 702 if (MPFR_ARE_SINGULAR (b, c)) 703 { 704 if (MPFR_IS_NAN (b) || MPFR_IS_NAN (c)) 705 { 706 MPFR_SET_NAN (a); 707 MPFR_RET_NAN; 708 } 709 sign = MPFR_MULT_SIGN (MPFR_SIGN (b), MPFR_SIGN (c)); 710 if (MPFR_IS_INF (b)) 711 { 712 if (!MPFR_IS_ZERO (c)) 713 { 714 MPFR_SET_SIGN (a, sign); 715 MPFR_SET_INF (a); 716 MPFR_RET (0); 717 } 718 else 719 { 720 MPFR_SET_NAN (a); 721 MPFR_RET_NAN; 722 } 723 } 724 else if (MPFR_IS_INF (c)) 725 { 726 if (!MPFR_IS_ZERO (b)) 727 { 728 MPFR_SET_SIGN (a, sign); 729 MPFR_SET_INF (a); 730 MPFR_RET(0); 731 } 732 else 733 { 734 MPFR_SET_NAN (a); 735 MPFR_RET_NAN; 736 } 737 } 738 else 739 { 740 MPFR_ASSERTD (MPFR_IS_ZERO(b) || MPFR_IS_ZERO(c)); 741 MPFR_SET_SIGN (a, sign); 742 MPFR_SET_ZERO (a); 743 MPFR_RET (0); 744 } 745 } 746 747 aq = MPFR_GET_PREC (a); 748 bq = MPFR_GET_PREC (b); 749 cq = MPFR_GET_PREC (c); 750 751 #if !defined(MPFR_GENERIC_ABI) 752 if (aq == bq && aq == cq) 753 { 754 if (aq < GMP_NUMB_BITS) 755 return mpfr_mul_1 (a, b, c, rnd_mode, aq); 756 757 if (GMP_NUMB_BITS < aq && aq < 2 * GMP_NUMB_BITS) 758 return mpfr_mul_2 (a, b, c, rnd_mode, aq); 759 760 if (aq == GMP_NUMB_BITS) 761 return mpfr_mul_1n (a, b, c, rnd_mode); 762 763 if (2 * GMP_NUMB_BITS < aq && aq < 3 * GMP_NUMB_BITS) 764 return mpfr_mul_3 (a, b, c, rnd_mode, aq); 765 } 766 #endif 767 768 sign = MPFR_MULT_SIGN (MPFR_SIGN (b), MPFR_SIGN (c)); 769 770 ax = MPFR_GET_EXP (b) + MPFR_GET_EXP (c); 771 /* Note: the exponent of the exact result will be e = bx + cx + ec with 772 ec in {-1,0,1} and the following assumes that e is representable. */ 773 774 /* FIXME: Useful since we do an exponent check after? 775 * It is useful iff the precision is big, there is an overflow 776 * and we are doing further mults...*/ 777 #ifdef HUGE 778 if (MPFR_UNLIKELY (ax > __gmpfr_emax + 1)) 779 return mpfr_overflow (a, rnd_mode, sign); 780 if (MPFR_UNLIKELY (ax < __gmpfr_emin - 2)) 781 return mpfr_underflow (a, rnd_mode == MPFR_RNDN ? MPFR_RNDZ : rnd_mode, 782 sign); 783 #endif 784 785 MPFR_ASSERTN ((mpfr_uprec_t) bq + cq <= MPFR_PREC_MAX); 786 787 bn = MPFR_PREC2LIMBS (bq); /* number of limbs of b */ 788 cn = MPFR_PREC2LIMBS (cq); /* number of limbs of c */ 789 k = bn + cn; /* effective nb of limbs used by b*c (= tn or tn+1) below */ 790 tn = MPFR_PREC2LIMBS (bq + cq); 791 MPFR_ASSERTD (tn <= k); /* tn <= k, thus no int overflow */ 792 793 /* Check for no size_t overflow. */ 794 MPFR_ASSERTD ((size_t) k <= ((size_t) -1) / MPFR_BYTES_PER_MP_LIMB); 795 MPFR_TMP_MARK (marker); 796 tmp = MPFR_TMP_LIMBS_ALLOC (k); 797 798 /* multiplies two mantissa in temporary allocated space */ 799 if (MPFR_UNLIKELY (bn < cn)) 800 { 801 mpfr_srcptr z = b; 802 mp_size_t zn = bn; 803 b = c; 804 bn = cn; 805 c = z; 806 cn = zn; 807 } 808 MPFR_ASSERTD (bn >= cn); 809 if (bn <= 2) 810 { 811 /* The 3 cases perform the same first operation. */ 812 umul_ppmm (tmp[1], tmp[0], MPFR_MANT (b)[0], MPFR_MANT (c)[0]); 813 if (bn == 1) 814 { 815 /* 1 limb * 1 limb */ 816 b1 = tmp[1]; 817 } 818 else if (MPFR_UNLIKELY (cn == 1)) 819 { 820 /* 2 limbs * 1 limb */ 821 mp_limb_t t; 822 umul_ppmm (tmp[2], t, MPFR_MANT (b)[1], MPFR_MANT (c)[0]); 823 add_ssaaaa (tmp[2], tmp[1], tmp[2], tmp[1], 0, t); 824 b1 = tmp[2]; 825 } 826 else 827 { 828 /* 2 limbs * 2 limbs */ 829 mp_limb_t t1, t2, t3; 830 /* First 2 limbs * 1 limb */ 831 umul_ppmm (tmp[2], t1, MPFR_MANT (b)[1], MPFR_MANT (c)[0]); 832 add_ssaaaa (tmp[2], tmp[1], tmp[2], tmp[1], 0, t1); 833 /* Second, the other 2 limbs * 1 limb product */ 834 umul_ppmm (t1, t2, MPFR_MANT (b)[0], MPFR_MANT (c)[1]); 835 umul_ppmm (tmp[3], t3, MPFR_MANT (b)[1], MPFR_MANT (c)[1]); 836 add_ssaaaa (tmp[3], t1, tmp[3], t1, 0, t3); 837 /* Sum those two partial products */ 838 add_ssaaaa (tmp[2], tmp[1], tmp[2], tmp[1], t1, t2); 839 tmp[3] += (tmp[2] < t1); 840 b1 = tmp[3]; 841 } 842 b1 >>= (GMP_NUMB_BITS - 1); 843 tmp += k - tn; 844 if (MPFR_UNLIKELY (b1 == 0)) 845 mpn_lshift (tmp, tmp, tn, 1); /* tn <= k, so no stack corruption */ 846 } 847 else 848 /* Mulders' mulhigh. This code can also be used via mpfr_sqr, 849 hence the tests b != c. */ 850 if (MPFR_UNLIKELY (cn > (threshold = b != c ? 851 MPFR_MUL_THRESHOLD : MPFR_SQR_THRESHOLD))) 852 { 853 mp_limb_t *bp, *cp; 854 mp_size_t n; 855 mpfr_prec_t p; 856 857 /* First check if we can reduce the precision of b or c: 858 exact values are a nightmare for the short product trick */ 859 bp = MPFR_MANT (b); 860 cp = MPFR_MANT (c); 861 MPFR_STAT_STATIC_ASSERT (MPFR_MUL_THRESHOLD >= 1 && 862 MPFR_SQR_THRESHOLD >= 1); 863 if (MPFR_UNLIKELY ((bp[0] == 0 && bp[1] == 0) || 864 (cp[0] == 0 && cp[1] == 0))) 865 { 866 mpfr_t b_tmp, c_tmp; 867 868 MPFR_TMP_FREE (marker); 869 /* Check for b */ 870 while (*bp == 0) 871 { 872 bp++; 873 bn--; 874 MPFR_ASSERTD (bn > 0); 875 } /* This must end since the most significant limb is != 0 */ 876 877 /* Check for c too: if b == c, this will do nothing */ 878 while (*cp == 0) 879 { 880 cp++; 881 cn--; 882 MPFR_ASSERTD (cn > 0); 883 } /* This must end since the most significant limb is != 0 */ 884 885 /* It is not the fastest way, but it is safer. */ 886 MPFR_SET_SAME_SIGN (b_tmp, b); 887 MPFR_SET_EXP (b_tmp, MPFR_GET_EXP (b)); 888 MPFR_PREC (b_tmp) = bn * GMP_NUMB_BITS; 889 MPFR_MANT (b_tmp) = bp; 890 891 if (b != c) 892 { 893 MPFR_SET_SAME_SIGN (c_tmp, c); 894 MPFR_SET_EXP (c_tmp, MPFR_GET_EXP (c)); 895 MPFR_PREC (c_tmp) = cn * GMP_NUMB_BITS; 896 MPFR_MANT (c_tmp) = cp; 897 898 /* Call again mpfr_mul with the fixed arguments */ 899 return mpfr_mul (a, b_tmp, c_tmp, rnd_mode); 900 } 901 else 902 /* Call mpfr_mul instead of mpfr_sqr as the precision 903 is probably still high enough. */ 904 return mpfr_mul (a, b_tmp, b_tmp, rnd_mode); 905 } 906 907 /* Compute estimated precision of mulhigh. 908 We could use `+ (n < cn) + (n < bn)' instead of `+ 2', 909 but does it worth it? */ 910 n = MPFR_LIMB_SIZE (a) + 1; 911 n = MIN (n, cn); 912 MPFR_ASSERTD (n >= 1 && 2*n <= k && n <= cn && n <= bn); 913 p = n * GMP_NUMB_BITS - MPFR_INT_CEIL_LOG2 (n + 2); 914 bp += bn - n; 915 cp += cn - n; 916 917 /* Check if MulHigh can produce a roundable result. 918 We may lose 1 bit due to RNDN, 1 due to final shift. */ 919 if (MPFR_UNLIKELY (aq > p - 5)) 920 { 921 if (MPFR_UNLIKELY (aq > p - 5 + GMP_NUMB_BITS 922 || bn <= threshold + 1)) 923 { 924 /* MulHigh can't produce a roundable result. */ 925 MPFR_LOG_MSG (("mpfr_mulhigh can't be used (%lu VS %lu)\n", 926 aq, p)); 927 goto full_multiply; 928 } 929 /* Add one extra limb to mantissa of b and c. */ 930 if (bn > n) 931 bp --; 932 else 933 { 934 bp = MPFR_TMP_LIMBS_ALLOC (n + 1); 935 bp[0] = 0; 936 MPN_COPY (bp + 1, MPFR_MANT (b) + bn - n, n); 937 } 938 if (b != c) 939 { 940 if (cn > n) 941 cp --; /* FIXME: Could this happen? */ 942 else 943 { 944 cp = MPFR_TMP_LIMBS_ALLOC (n + 1); 945 cp[0] = 0; 946 MPN_COPY (cp + 1, MPFR_MANT (c) + cn - n, n); 947 } 948 } 949 /* We will compute with one extra limb */ 950 n++; 951 /* ceil(log2(n+2)) takes into account the lost bits due to 952 Mulders' short product */ 953 p = n * GMP_NUMB_BITS - MPFR_INT_CEIL_LOG2 (n + 2); 954 /* Due to some nasty reasons we can have only 4 bits */ 955 MPFR_ASSERTD (aq <= p - 4); 956 957 if (MPFR_LIKELY (k < 2*n)) 958 { 959 tmp = MPFR_TMP_LIMBS_ALLOC (2 * n); 960 tmp += 2*n-k; /* `tmp' still points to an area of `k' limbs */ 961 } 962 } 963 MPFR_LOG_MSG (("Use mpfr_mulhigh (%lu VS %lu)\n", aq, p)); 964 /* Compute an approximation of the product of b and c */ 965 if (b != c) 966 mpfr_mulhigh_n (tmp + k - 2 * n, bp, cp, n); 967 else 968 mpfr_sqrhigh_n (tmp + k - 2 * n, bp, n); 969 /* now tmp[k-n]..tmp[k-1] contains an approximation of the n upper 970 limbs of the product, with tmp[k-1] >= 2^(GMP_NUMB_BITS-2) */ 971 b1 = tmp[k-1] >> (GMP_NUMB_BITS - 1); /* msb from the product */ 972 973 /* If the mantissas of b and c are uniformly distributed in (1/2, 1], 974 then their product is in (1/4, 1/2] with probability 2*ln(2)-1 975 ~ 0.386 and in [1/2, 1] with probability 2-2*ln(2) ~ 0.614 */ 976 if (MPFR_UNLIKELY (b1 == 0)) 977 /* Warning: the mpfr_mulhigh_n call above only surely affects 978 tmp[k-n-1..k-1], thus we shift only those limbs */ 979 mpn_lshift (tmp + k - n - 1, tmp + k - n - 1, n + 1, 1); 980 tmp += k - tn; 981 /* now the approximation is in tmp[tn-n]...tmp[tn-1] */ 982 MPFR_ASSERTD (MPFR_LIMB_MSB (tmp[tn-1]) != 0); 983 984 /* for RNDF, we simply use RNDZ, since anyway here we multiply numbers 985 with large precisions, thus the overhead of RNDZ is small */ 986 if (rnd_mode == MPFR_RNDF) 987 rnd_mode = MPFR_RNDZ; 988 989 /* if the most significant bit b1 is zero, we have only p-1 correct 990 bits */ 991 if (MPFR_UNLIKELY (!mpfr_round_p (tmp, tn, p + b1 - 1, 992 aq + (rnd_mode == MPFR_RNDN)))) 993 { 994 tmp -= k - tn; /* tmp may have changed, FIX IT!!!!! */ 995 goto full_multiply; 996 } 997 } 998 else 999 { 1000 full_multiply: 1001 MPFR_LOG_MSG (("Use mpn_mul\n", 0)); 1002 b1 = mpn_mul (tmp, MPFR_MANT (b), bn, MPFR_MANT (c), cn); 1003 1004 /* now tmp[0]..tmp[k-1] contains the product of both mantissa, 1005 with tmp[k-1]>=2^(GMP_NUMB_BITS-2) */ 1006 b1 >>= GMP_NUMB_BITS - 1; /* msb from the product */ 1007 1008 /* if the mantissas of b and c are uniformly distributed in (1/2, 1], 1009 then their product is in (1/4, 1/2] with probability 2*ln(2)-1 1010 ~ 0.386 and in [1/2, 1] with probability 2-2*ln(2) ~ 0.614 */ 1011 tmp += k - tn; 1012 if (MPFR_UNLIKELY (b1 == 0)) 1013 mpn_lshift (tmp, tmp, tn, 1); /* tn <= k, so no stack corruption */ 1014 } 1015 1016 ax2 = ax + (mpfr_exp_t) (b1 - 1); 1017 MPFR_RNDRAW (inexact, a, tmp, bq + cq, rnd_mode, sign, ax2++); 1018 MPFR_TMP_FREE (marker); 1019 MPFR_EXP (a) = ax2; /* Can't use MPFR_SET_EXP: Expo may be out of range */ 1020 MPFR_SET_SIGN (a, sign); 1021 if (MPFR_UNLIKELY (ax2 > __gmpfr_emax)) 1022 return mpfr_overflow (a, rnd_mode, sign); 1023 if (MPFR_UNLIKELY (ax2 < __gmpfr_emin)) 1024 { 1025 /* In the rounding to the nearest mode, if the exponent of the exact 1026 result (i.e. before rounding, i.e. without taking cc into account) 1027 is < __gmpfr_emin - 1 or the exact result is a power of 2 (i.e. if 1028 both arguments are powers of 2), then round to zero. */ 1029 if (rnd_mode == MPFR_RNDN 1030 && (ax + (mpfr_exp_t) b1 < __gmpfr_emin 1031 || (mpfr_powerof2_raw (b) && mpfr_powerof2_raw (c)))) 1032 rnd_mode = MPFR_RNDZ; 1033 return mpfr_underflow (a, rnd_mode, sign); 1034 } 1035 MPFR_RET (inexact); 1036 } 1037