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 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
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_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
mpfr_pow(mpfr_ptr z,mpfr_srcptr x,mpfr_srcptr y,mpfr_rnd_t rnd_mode)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