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