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