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