1 /* mpfr_pow -- power function x^y 2 3 Copyright 2001-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 /* return non zero iff x^y is exact. 27 Assumes x and y are ordinary numbers, 28 y is not an integer, x is not a power of 2 and x is positive 29 30 If x^y is exact, it computes it and sets *inexact. 31 */ 32 static int 33 mpfr_pow_is_exact (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, 34 mpfr_rnd_t rnd_mode, int *inexact) 35 { 36 mpz_t a, c; 37 mpfr_exp_t d, b; 38 unsigned long i; 39 int res; 40 41 MPFR_ASSERTD (!MPFR_IS_SINGULAR (y)); 42 MPFR_ASSERTD (!MPFR_IS_SINGULAR (x)); 43 MPFR_ASSERTD (!mpfr_integer_p (y)); 44 MPFR_ASSERTD (mpfr_cmp_si_2exp (x, MPFR_INT_SIGN (x), 45 MPFR_GET_EXP (x) - 1) != 0); 46 MPFR_ASSERTD (MPFR_IS_POS (x)); 47 48 if (MPFR_IS_NEG (y)) 49 return 0; /* x is not a power of two => x^-y is not exact */ 50 51 /* compute d such that y = c*2^d with c odd integer */ 52 mpz_init (c); 53 d = mpfr_get_z_2exp (c, y); 54 i = mpz_scan1 (c, 0); 55 mpz_fdiv_q_2exp (c, c, i); 56 d += i; 57 /* now y=c*2^d with c odd */ 58 /* Since y is not an integer, d is necessarily < 0 */ 59 MPFR_ASSERTD (d < 0); 60 61 /* Compute a,b such that x=a*2^b */ 62 mpz_init (a); 63 b = mpfr_get_z_2exp (a, x); 64 i = mpz_scan1 (a, 0); 65 mpz_fdiv_q_2exp (a, a, i); 66 b += i; 67 /* now x=a*2^b with a is odd */ 68 69 for (res = 1 ; d != 0 ; d++) 70 { 71 /* a*2^b is a square iff 72 (i) a is a square when b is even 73 (ii) 2*a is a square when b is odd */ 74 if (b % 2 != 0) 75 { 76 mpz_mul_2exp (a, a, 1); /* 2*a */ 77 b --; 78 } 79 MPFR_ASSERTD ((b % 2) == 0); 80 if (!mpz_perfect_square_p (a)) 81 { 82 res = 0; 83 goto end; 84 } 85 mpz_sqrt (a, a); 86 b = b / 2; 87 } 88 /* Now x = (a'*2^b')^(2^-d) with d < 0 89 so x^y = ((a'*2^b')^(2^-d))^(c*2^d) 90 = ((a'*2^b')^c with c odd integer */ 91 { 92 mpfr_t tmp; 93 mpfr_prec_t p; 94 MPFR_MPZ_SIZEINBASE2 (p, a); 95 mpfr_init2 (tmp, p); /* prec = 1 should not be possible */ 96 res = mpfr_set_z (tmp, a, MPFR_RNDN); 97 MPFR_ASSERTD (res == 0); 98 res = mpfr_mul_2si (tmp, tmp, b, MPFR_RNDN); 99 MPFR_ASSERTD (res == 0); 100 *inexact = mpfr_pow_z (z, tmp, c, rnd_mode); 101 mpfr_clear (tmp); 102 res = 1; 103 } 104 end: 105 mpz_clear (a); 106 mpz_clear (c); 107 return res; 108 } 109 110 /* Assumes that the exponent range has already been extended and if y is 111 an integer, then the result is not exact in unbounded exponent range. 112 If x < 0, assumes y is an integer. 113 */ 114 int 115 mpfr_pow_general (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, 116 mpfr_rnd_t rnd_mode, int y_is_integer, mpfr_save_expo_t *expo) 117 { 118 mpfr_t t, u, k, absx; 119 int neg_result = 0; 120 int k_non_zero = 0; 121 int check_exact_case = 0; 122 int inexact; 123 /* Declaration of the size variable */ 124 mpfr_prec_t Nz = MPFR_PREC(z); /* target precision */ 125 mpfr_prec_t Nt; /* working precision */ 126 mpfr_exp_t err; /* error */ 127 MPFR_ZIV_DECL (ziv_loop); 128 129 130 MPFR_LOG_FUNC 131 (("x[%Pu]=%.*Rg y[%Pu]=%.*Rg rnd=%d", 132 mpfr_get_prec (x), mpfr_log_prec, x, 133 mpfr_get_prec (y), mpfr_log_prec, y, rnd_mode), 134 ("z[%Pu]=%.*Rg inexact=%d", 135 mpfr_get_prec (z), mpfr_log_prec, z, inexact)); 136 137 /* We put the absolute value of x in absx, pointing to the significand 138 of x to avoid allocating memory for the significand of absx. */ 139 MPFR_ALIAS(absx, x, /*sign=*/ 1, /*EXP=*/ MPFR_EXP(x)); 140 141 /* We will compute the absolute value of the result. So, let's 142 invert the rounding mode if the result is negative (in which case 143 y not an integer was already filtered out). */ 144 if (MPFR_IS_NEG (x)) 145 { 146 MPFR_ASSERTD (y_is_integer); 147 if (mpfr_odd_p (y)) 148 { 149 neg_result = 1; 150 rnd_mode = MPFR_INVERT_RND (rnd_mode); 151 } 152 } 153 154 /* Compute the precision of intermediary variable. */ 155 /* The increment 9 + MPFR_INT_CEIL_LOG2 (Nz) gives few Ziv failures 156 in binary64 and binary128 formats: 157 mfv5 -p53 -e1 mpfr_pow: 5903 / 6469.59 / 6686 158 mfv5 -p113 -e1 mpfr_pow: 10913 / 11989.46 / 12321 */ 159 Nt = Nz + 9 + MPFR_INT_CEIL_LOG2 (Nz); 160 161 /* initialize of intermediary variable */ 162 mpfr_init2 (t, Nt); 163 164 MPFR_ZIV_INIT (ziv_loop, Nt); 165 for (;;) 166 { 167 MPFR_BLOCK_DECL (flags1); 168 169 /* compute exp(y*ln|x|), using MPFR_RNDU to get an upper bound, so 170 that we can detect underflows. */ 171 mpfr_log (t, absx, MPFR_IS_NEG (y) ? MPFR_RNDD : MPFR_RNDU); /* ln|x| */ 172 mpfr_mul (t, y, t, MPFR_RNDU); /* y*ln|x| */ 173 if (k_non_zero) 174 { 175 MPFR_LOG_MSG (("subtract k * ln(2)\n", 0)); 176 mpfr_const_log2 (u, MPFR_RNDD); 177 mpfr_mul (u, u, k, MPFR_RNDD); 178 /* Error on u = k * log(2): < k * 2^(-Nt) < 1. */ 179 mpfr_sub (t, t, u, MPFR_RNDU); 180 MPFR_LOG_MSG (("t = y * ln|x| - k * ln(2)\n", 0)); 181 MPFR_LOG_VAR (t); 182 } 183 /* estimate of the error -- see pow function in algorithms.tex. 184 The error on t is at most 1/2 + 3*2^(EXP(t)+1) ulps, which is 185 <= 2^(EXP(t)+3) for EXP(t) >= -1, and <= 2 ulps for EXP(t) <= -2. 186 Additional error if k_no_zero: treal = t * errk, with 187 1 - |k| * 2^(-Nt) <= exp(-|k| * 2^(-Nt)) <= errk <= 1, 188 i.e., additional absolute error <= 2^(EXP(k)+EXP(t)-Nt). 189 Total error <= 2^err1 + 2^err2 <= 2^(max(err1,err2)+1). */ 190 err = MPFR_NOTZERO (t) && MPFR_GET_EXP (t) >= -1 ? 191 MPFR_GET_EXP (t) + 3 : 1; 192 if (k_non_zero) 193 { 194 if (MPFR_GET_EXP (k) > err) 195 err = MPFR_GET_EXP (k); 196 err++; 197 } 198 MPFR_BLOCK (flags1, mpfr_exp (t, t, MPFR_RNDN)); /* exp(y*ln|x|)*/ 199 /* We need to test */ 200 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (t) || MPFR_UNDERFLOW (flags1))) 201 { 202 mpfr_prec_t Ntmin; 203 MPFR_BLOCK_DECL (flags2); 204 205 MPFR_ASSERTN (!k_non_zero); 206 MPFR_ASSERTN (!MPFR_IS_NAN (t)); 207 208 /* Real underflow? */ 209 if (MPFR_IS_ZERO (t)) 210 { 211 /* Underflow. We computed rndn(exp(t)), where t >= y*ln|x|. 212 Therefore rndn(|x|^y) = 0, and we have a real underflow on 213 |x|^y. */ 214 inexact = mpfr_underflow (z, rnd_mode == MPFR_RNDN ? MPFR_RNDZ 215 : rnd_mode, MPFR_SIGN_POS); 216 if (expo != NULL) 217 MPFR_SAVE_EXPO_UPDATE_FLAGS (*expo, MPFR_FLAGS_INEXACT 218 | MPFR_FLAGS_UNDERFLOW); 219 break; 220 } 221 222 /* Real overflow? */ 223 if (MPFR_IS_INF (t)) 224 { 225 /* Note: we can probably use a low precision for this test. */ 226 mpfr_log (t, absx, MPFR_IS_NEG (y) ? MPFR_RNDU : MPFR_RNDD); 227 mpfr_mul (t, y, t, MPFR_RNDD); /* y * ln|x| */ 228 MPFR_BLOCK (flags2, mpfr_exp (t, t, MPFR_RNDD)); 229 /* t = lower bound on exp(y * ln|x|) */ 230 if (MPFR_OVERFLOW (flags2)) 231 { 232 /* We have computed a lower bound on |x|^y, and it 233 overflowed. Therefore we have a real overflow 234 on |x|^y. */ 235 inexact = mpfr_overflow (z, rnd_mode, MPFR_SIGN_POS); 236 if (expo != NULL) 237 MPFR_SAVE_EXPO_UPDATE_FLAGS (*expo, MPFR_FLAGS_INEXACT 238 | MPFR_FLAGS_OVERFLOW); 239 break; 240 } 241 } 242 243 k_non_zero = 1; 244 Ntmin = sizeof(mpfr_exp_t) * CHAR_BIT; 245 if (Ntmin > Nt) 246 { 247 Nt = Ntmin; 248 mpfr_set_prec (t, Nt); 249 } 250 mpfr_init2 (u, Nt); 251 mpfr_init2 (k, Ntmin); 252 mpfr_log2 (k, absx, MPFR_RNDN); 253 mpfr_mul (k, y, k, MPFR_RNDN); 254 mpfr_round (k, k); 255 MPFR_LOG_VAR (k); 256 /* |y| < 2^Ntmin, therefore |k| < 2^Nt. */ 257 continue; 258 } 259 if (MPFR_LIKELY (MPFR_CAN_ROUND (t, Nt - err, Nz, rnd_mode))) 260 { 261 inexact = mpfr_set (z, t, rnd_mode); 262 break; 263 } 264 265 /* check exact power, except when y is an integer (since the 266 exact cases for y integer have already been filtered out) */ 267 if (check_exact_case == 0 && ! y_is_integer) 268 { 269 if (mpfr_pow_is_exact (z, absx, y, rnd_mode, &inexact)) 270 break; 271 check_exact_case = 1; 272 } 273 274 /* reactualisation of the precision */ 275 MPFR_ZIV_NEXT (ziv_loop, Nt); 276 mpfr_set_prec (t, Nt); 277 if (k_non_zero) 278 mpfr_set_prec (u, Nt); 279 } 280 MPFR_ZIV_FREE (ziv_loop); 281 282 if (k_non_zero) 283 { 284 int inex2; 285 long lk; 286 287 /* The rounded result in an unbounded exponent range is z * 2^k. As 288 * MPFR chooses underflow after rounding, the mpfr_mul_2si below will 289 * correctly detect underflows and overflows. However, in rounding to 290 * nearest, if z * 2^k = 2^(emin - 2), then the double rounding may 291 * affect the result. We need to cope with that before overwriting z. 292 * This can occur only if k < 0 (this test is necessary to avoid a 293 * potential integer overflow). 294 * If inexact >= 0, then the real result is <= 2^(emin - 2), so that 295 * o(2^(emin - 2)) = +0 is correct. If inexact < 0, then the real 296 * result is > 2^(emin - 2) and we need to round to 2^(emin - 1). 297 */ 298 MPFR_ASSERTN (MPFR_EXP_MAX <= LONG_MAX); 299 lk = mpfr_get_si (k, MPFR_RNDN); 300 /* Due to early overflow detection, |k| should not be much larger than 301 * MPFR_EMAX_MAX, and as MPFR_EMAX_MAX <= MPFR_EXP_MAX/2 <= LONG_MAX/2, 302 * an overflow should not be possible in mpfr_get_si (and lk is exact). 303 * And one even has the following assertion. TODO: complete proof. 304 */ 305 MPFR_ASSERTD (lk > LONG_MIN && lk < LONG_MAX); 306 /* Note: even in case of overflow (lk inexact), the code is correct. 307 * Indeed, for the 3 occurrences of lk: 308 * - The test lk < 0 is correct as sign(lk) = sign(k). 309 * - In the test MPFR_GET_EXP (z) == __gmpfr_emin - 1 - lk, 310 * if lk is inexact, then lk = LONG_MIN <= MPFR_EXP_MIN 311 * (the minimum value of the mpfr_exp_t type), and 312 * __gmpfr_emin - 1 - lk >= MPFR_EMIN_MIN - 1 - 2 * MPFR_EMIN_MIN 313 * >= - MPFR_EMIN_MIN - 1 = MPFR_EMAX_MAX - 1. However, from the 314 * choice of k, z has been chosen to be around 1, so that the 315 * result of the test is false, as if lk were exact. 316 * - In the mpfr_mul_2si (z, z, lk, rnd_mode), if lk is inexact, 317 * then |lk| >= LONG_MAX >= MPFR_EXP_MAX, and as z is around 1, 318 * mpfr_mul_2si underflows or overflows in the same way as if 319 * lk were exact. 320 * TODO: give a bound on |t|, then on |EXP(z)|. 321 */ 322 if (rnd_mode == MPFR_RNDN && inexact < 0 && lk < 0 && 323 MPFR_GET_EXP (z) == __gmpfr_emin - 1 - lk && mpfr_powerof2_raw (z)) 324 { 325 /* Rounding to nearest, real result > z * 2^k = 2^(emin - 2), 326 * underflow case: 327 * (a) if the precision of z is > 1, we will obtain the correct 328 * result and exceptions by replacing z by nextabove(z). 329 * (b) if the precision of z is 1, we first copy z to zcopy of 330 * precision 2 bits and perform nextabove(zcopy). 331 */ 332 if (MPFR_PREC(z) >= 2) 333 mpfr_nextabove (z); 334 else 335 { 336 mpfr_t zcopy; 337 mpfr_init2 (zcopy, MPFR_PREC(z) + 1); 338 mpfr_set (zcopy, z, MPFR_RNDZ); 339 mpfr_nextabove (zcopy); 340 inex2 = mpfr_mul_2si (z, zcopy, lk, rnd_mode); 341 mpfr_clear (zcopy); 342 goto under_over; 343 } 344 } 345 MPFR_CLEAR_FLAGS (); 346 inex2 = mpfr_mul_2si (z, z, lk, rnd_mode); 347 under_over: 348 if (inex2) /* underflow or overflow */ 349 { 350 inexact = inex2; 351 if (expo != NULL) 352 MPFR_SAVE_EXPO_UPDATE_FLAGS (*expo, __gmpfr_flags); 353 } 354 mpfr_clears (u, k, (mpfr_ptr) 0); 355 } 356 mpfr_clear (t); 357 358 /* update the sign of the result if x was negative */ 359 if (neg_result) 360 { 361 MPFR_SET_NEG(z); 362 inexact = -inexact; 363 } 364 365 return inexact; 366 } 367 368 /* The computation of z = pow(x,y) is done by 369 z = exp(y * log(x)) = x^y 370 For the special cases, see Section F.9.4.4 of the C standard: 371 _ pow(±0, y) = ±inf for y an odd integer < 0. 372 _ pow(±0, y) = +inf for y < 0 and not an odd integer. 373 _ pow(±0, y) = ±0 for y an odd integer > 0. 374 _ pow(±0, y) = +0 for y > 0 and not an odd integer. 375 _ pow(-1, ±inf) = 1. 376 _ pow(+1, y) = 1 for any y, even a NaN. 377 _ pow(x, ±0) = 1 for any x, even a NaN. 378 _ pow(x, y) = NaN for finite x < 0 and finite non-integer y. 379 _ pow(x, -inf) = +inf for |x| < 1. 380 _ pow(x, -inf) = +0 for |x| > 1. 381 _ pow(x, +inf) = +0 for |x| < 1. 382 _ pow(x, +inf) = +inf for |x| > 1. 383 _ pow(-inf, y) = -0 for y an odd integer < 0. 384 _ pow(-inf, y) = +0 for y < 0 and not an odd integer. 385 _ pow(-inf, y) = -inf for y an odd integer > 0. 386 _ pow(-inf, y) = +inf for y > 0 and not an odd integer. 387 _ pow(+inf, y) = +0 for y < 0. 388 _ pow(+inf, y) = +inf for y > 0. */ 389 int 390 mpfr_pow (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mpfr_rnd_t rnd_mode) 391 { 392 int inexact; 393 int cmp_x_1; 394 int y_is_integer; 395 MPFR_SAVE_EXPO_DECL (expo); 396 397 MPFR_LOG_FUNC 398 (("x[%Pu]=%.*Rg y[%Pu]=%.*Rg rnd=%d", 399 mpfr_get_prec (x), mpfr_log_prec, x, 400 mpfr_get_prec (y), mpfr_log_prec, y, rnd_mode), 401 ("z[%Pu]=%.*Rg inexact=%d", 402 mpfr_get_prec (z), mpfr_log_prec, z, inexact)); 403 404 if (MPFR_ARE_SINGULAR (x, y)) 405 { 406 /* pow(x, 0) returns 1 for any x, even a NaN. */ 407 if (MPFR_UNLIKELY (MPFR_IS_ZERO (y))) 408 return mpfr_set_ui (z, 1, rnd_mode); 409 else if (MPFR_IS_NAN (x)) 410 { 411 MPFR_SET_NAN (z); 412 MPFR_RET_NAN; 413 } 414 else if (MPFR_IS_NAN (y)) 415 { 416 /* pow(+1, NaN) returns 1. */ 417 if (mpfr_cmp_ui (x, 1) == 0) 418 return mpfr_set_ui (z, 1, rnd_mode); 419 MPFR_SET_NAN (z); 420 MPFR_RET_NAN; 421 } 422 else if (MPFR_IS_INF (y)) 423 { 424 if (MPFR_IS_INF (x)) 425 { 426 if (MPFR_IS_POS (y)) 427 MPFR_SET_INF (z); 428 else 429 MPFR_SET_ZERO (z); 430 MPFR_SET_POS (z); 431 MPFR_RET (0); 432 } 433 else 434 { 435 int cmp; 436 cmp = mpfr_cmpabs (x, __gmpfr_one) * MPFR_INT_SIGN (y); 437 MPFR_SET_POS (z); 438 if (cmp > 0) 439 { 440 /* Return +inf. */ 441 MPFR_SET_INF (z); 442 MPFR_RET (0); 443 } 444 else if (cmp < 0) 445 { 446 /* Return +0. */ 447 MPFR_SET_ZERO (z); 448 MPFR_RET (0); 449 } 450 else 451 { 452 /* Return 1. */ 453 return mpfr_set_ui (z, 1, rnd_mode); 454 } 455 } 456 } 457 else if (MPFR_IS_INF (x)) 458 { 459 int negative; 460 /* Determine the sign now, in case y and z are the same object */ 461 negative = MPFR_IS_NEG (x) && mpfr_odd_p (y); 462 if (MPFR_IS_POS (y)) 463 MPFR_SET_INF (z); 464 else 465 MPFR_SET_ZERO (z); 466 if (negative) 467 MPFR_SET_NEG (z); 468 else 469 MPFR_SET_POS (z); 470 MPFR_RET (0); 471 } 472 else 473 { 474 int negative; 475 MPFR_ASSERTD (MPFR_IS_ZERO (x)); 476 /* Determine the sign now, in case y and z are the same object */ 477 negative = MPFR_IS_NEG(x) && mpfr_odd_p (y); 478 if (MPFR_IS_NEG (y)) 479 { 480 MPFR_ASSERTD (! MPFR_IS_INF (y)); 481 MPFR_SET_INF (z); 482 MPFR_SET_DIVBY0 (); 483 } 484 else 485 MPFR_SET_ZERO (z); 486 if (negative) 487 MPFR_SET_NEG (z); 488 else 489 MPFR_SET_POS (z); 490 MPFR_RET (0); 491 } 492 } 493 494 /* x^y for x < 0 and y not an integer is not defined */ 495 y_is_integer = mpfr_integer_p (y); 496 if (MPFR_IS_NEG (x) && ! y_is_integer) 497 { 498 MPFR_SET_NAN (z); 499 MPFR_RET_NAN; 500 } 501 502 /* now the result cannot be NaN: 503 (1) either x > 0 504 (2) or x < 0 and y is an integer */ 505 506 cmp_x_1 = mpfr_cmpabs (x, __gmpfr_one); 507 if (cmp_x_1 == 0) 508 return mpfr_set_si (z, MPFR_IS_NEG (x) && mpfr_odd_p (y) ? -1 : 1, rnd_mode); 509 510 /* now we have: 511 (1) either x > 0 512 (2) or x < 0 and y is an integer 513 and in addition |x| <> 1. 514 */ 515 516 /* detect overflow: an overflow is possible if 517 (a) |x| > 1 and y > 0 518 (b) |x| < 1 and y < 0. 519 FIXME: this assumes 1 is always representable. 520 521 FIXME2: maybe we can test overflow and underflow simultaneously. 522 The idea is the following: first compute an approximation to 523 y * log2|x|, using rounding to nearest. If |x| is not too near from 1, 524 this approximation should be accurate enough, and in most cases enable 525 one to prove that there is no underflow nor overflow. 526 Otherwise, it should enable one to check only underflow or overflow, 527 instead of both cases as in the present case. 528 */ 529 530 /* fast check for cases where no overflow nor underflow is possible: 531 if |y| <= 2^15, and -32767 < EXP(x) <= 32767, then 532 |y*log2(x)| <= 2^15*32767 < 1073741823, thus for the default 533 emax=1073741823 and emin=-emax there can be no overflow nor underflow */ 534 if (__gmpfr_emax >= 1073741823 && __gmpfr_emin <= -1073741823 && 535 MPFR_EXP(y) <= 15 && -32767 < MPFR_EXP(x) && MPFR_EXP(x) <= 32767) 536 goto no_overflow_nor_underflow; 537 538 if (cmp_x_1 * MPFR_SIGN (y) > 0) 539 { 540 mpfr_t t; 541 int negative, overflow; 542 543 MPFR_SAVE_EXPO_MARK (expo); 544 mpfr_init2 (t, 53); 545 /* we want a lower bound on y*log2|x|: 546 (i) if x > 0, it suffices to round log2(x) toward zero, and 547 to round y*o(log2(x)) toward zero too; 548 (ii) if x < 0, we first compute t = o(-x), with rounding toward 1, 549 and then follow as in case (1). */ 550 if (MPFR_IS_POS (x)) 551 mpfr_log2 (t, x, MPFR_RNDZ); 552 else 553 { 554 mpfr_neg (t, x, (cmp_x_1 > 0) ? MPFR_RNDZ : MPFR_RNDU); 555 mpfr_log2 (t, t, MPFR_RNDZ); 556 } 557 mpfr_mul (t, t, y, MPFR_RNDZ); 558 overflow = mpfr_cmp_si (t, __gmpfr_emax) > 0; 559 mpfr_clear (t); 560 MPFR_SAVE_EXPO_FREE (expo); 561 if (overflow) 562 { 563 MPFR_LOG_MSG (("early overflow detection\n", 0)); 564 negative = MPFR_IS_NEG (x) && mpfr_odd_p (y); 565 return mpfr_overflow (z, rnd_mode, negative ? -1 : 1); 566 } 567 } 568 569 /* Basic underflow checking. One has: 570 * - if y > 0, |x^y| < 2^(EXP(x) * y); 571 * - if y < 0, |x^y| <= 2^((EXP(x) - 1) * y); 572 * so that one can compute a value ebound such that |x^y| < 2^ebound. 573 * If we have ebound <= emin - 2 (emin - 1 in directed rounding modes), 574 * then there is an underflow and we can decide the return value. 575 */ 576 if (MPFR_IS_NEG (y) ? (MPFR_GET_EXP (x) > 1) : (MPFR_GET_EXP (x) < 0)) 577 { 578 mp_limb_t tmp_limb[MPFR_EXP_LIMB_SIZE]; 579 mpfr_t tmp; 580 mpfr_eexp_t ebound; 581 int inex2; 582 583 /* We must restore the flags. */ 584 MPFR_SAVE_EXPO_MARK (expo); 585 MPFR_TMP_INIT1 (tmp_limb, tmp, sizeof (mpfr_exp_t) * CHAR_BIT); 586 inex2 = mpfr_set_exp_t (tmp, MPFR_GET_EXP (x), MPFR_RNDN); 587 MPFR_ASSERTN (inex2 == 0); 588 if (MPFR_IS_NEG (y)) 589 { 590 inex2 = mpfr_sub_ui (tmp, tmp, 1, MPFR_RNDN); 591 MPFR_ASSERTN (inex2 == 0); 592 } 593 mpfr_mul (tmp, tmp, y, MPFR_RNDU); 594 if (MPFR_IS_NEG (y)) 595 mpfr_nextabove (tmp); 596 /* tmp doesn't necessarily fit in ebound, but that doesn't matter 597 since we get the minimum value in such a case. */ 598 ebound = mpfr_get_exp_t (tmp, MPFR_RNDU); 599 MPFR_SAVE_EXPO_FREE (expo); 600 if (MPFR_UNLIKELY (ebound <= 601 __gmpfr_emin - (rnd_mode == MPFR_RNDN ? 2 : 1))) 602 { 603 /* warning: mpfr_underflow rounds away from 0 for MPFR_RNDN */ 604 MPFR_LOG_MSG (("early underflow detection\n", 0)); 605 return mpfr_underflow (z, 606 rnd_mode == MPFR_RNDN ? MPFR_RNDZ : rnd_mode, 607 MPFR_IS_NEG (x) && mpfr_odd_p (y) ? -1 : 1); 608 } 609 } 610 611 no_overflow_nor_underflow: 612 613 /* If y is an integer, we can use mpfr_pow_z (based on multiplications), 614 but if y is very large (I'm not sure about the best threshold -- VL), 615 we shouldn't use it, as it can be very slow and take a lot of memory 616 (and even crash or make other programs crash, as several hundred of 617 MBs may be necessary). Note that in such a case, either x = +/-2^b 618 (this case is handled below) or x^y cannot be represented exactly in 619 any precision supported by MPFR (the general case uses this property). 620 */ 621 if (y_is_integer && (MPFR_GET_EXP (y) <= 256)) 622 { 623 mpz_t zi; 624 625 MPFR_LOG_MSG (("special code for y not too large integer\n", 0)); 626 mpz_init (zi); 627 mpfr_get_z (zi, y, MPFR_RNDN); 628 inexact = mpfr_pow_z (z, x, zi, rnd_mode); 629 mpz_clear (zi); 630 return inexact; 631 } 632 633 /* Special case (+/-2^b)^Y which could be exact. If x is negative, then 634 necessarily y is a large integer. */ 635 if (mpfr_powerof2_raw (x)) 636 { 637 mpfr_exp_t b = MPFR_GET_EXP (x) - 1; 638 mpfr_t tmp; 639 int sgnx = MPFR_SIGN (x); 640 641 MPFR_ASSERTN (b >= LONG_MIN && b <= LONG_MAX); /* FIXME... */ 642 643 MPFR_LOG_MSG (("special case (+/-2^b)^Y\n", 0)); 644 /* now x = +/-2^b, so x^y = (+/-1)^y*2^(b*y) is exact whenever b*y is 645 an integer */ 646 MPFR_SAVE_EXPO_MARK (expo); 647 mpfr_init2 (tmp, MPFR_PREC (y) + sizeof (long) * CHAR_BIT); 648 inexact = mpfr_mul_si (tmp, y, b, MPFR_RNDN); /* exact */ 649 MPFR_ASSERTN (inexact == 0); 650 /* Note: as the exponent range has been extended, an overflow is not 651 possible (due to basic overflow and underflow checking above, as 652 the result is ~ 2^tmp), and an underflow is not possible either 653 because b is an integer (thus either 0 or >= 1). */ 654 MPFR_CLEAR_FLAGS (); 655 inexact = mpfr_exp2 (z, tmp, rnd_mode); 656 mpfr_clear (tmp); 657 if (sgnx < 0 && mpfr_odd_p (y)) 658 { 659 mpfr_neg (z, z, rnd_mode); 660 inexact = -inexact; 661 } 662 /* Without the following, the overflows3 test in tpow.c fails. */ 663 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags); 664 MPFR_SAVE_EXPO_FREE (expo); 665 return mpfr_check_range (z, inexact, rnd_mode); 666 } 667 668 MPFR_SAVE_EXPO_MARK (expo); 669 670 /* Case where y * log(x) is very small. Warning: x can be negative, in 671 that case y is a large integer. */ 672 { 673 mpfr_exp_t err, expx, logt; 674 675 /* We need an upper bound on the exponent of y * log(x). */ 676 if (MPFR_IS_POS(x)) 677 expx = cmp_x_1 > 0 ? MPFR_EXP(x) : 1 - MPFR_EXP(x); 678 else 679 expx = mpfr_cmp_si (x, -1) > 0 ? 1 - MPFR_EXP(x) : MPFR_EXP(x); 680 MPFR_ASSERTD(expx >= 0); 681 /* now |log(x)| < expx */ 682 logt = MPFR_INT_CEIL_LOG2 (expx); 683 /* now expx <= 2^logt */ 684 err = MPFR_GET_EXP (y) + logt; 685 MPFR_CLEAR_FLAGS (); 686 MPFR_SMALL_INPUT_AFTER_SAVE_EXPO (z, __gmpfr_one, - err, 0, 687 (MPFR_IS_POS (y)) ^ (cmp_x_1 < 0), 688 rnd_mode, expo, {}); 689 } 690 691 /* General case */ 692 inexact = mpfr_pow_general (z, x, y, rnd_mode, y_is_integer, &expo); 693 694 MPFR_SAVE_EXPO_FREE (expo); 695 return mpfr_check_range (z, inexact, rnd_mode); 696 } 697