--- /dev/null
+/* mpfr_pow -- power function x^y
+
+Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
+Contributed by the Arenaire and Cacao projects, INRIA.
+
+This file is part of the GNU MPFR Library.
+
+The GNU MPFR Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published by
+the Free Software Foundation; either version 2.1 of the License, or (at your
+option) any later version.
+
+The GNU MPFR Library is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
+License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with the GNU MPFR Library; see the file COPYING.LIB. If not, write to
+the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
+MA 02110-1301, USA. */
+
+#define MPFR_NEED_LONGLONG_H
+#include "mpfr-impl.h"
+
+/* return non zero iff x^y is exact.
+ Assumes x and y are ordinary numbers,
+ y is not an integer, x is not a power of 2 and x is positive
+
+ If x^y is exact, it computes it and sets *inexact.
+*/
+static int
+mpfr_pow_is_exact (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y,
+ mp_rnd_t rnd_mode, int *inexact)
+{
+ mpz_t a, c;
+ mp_exp_t d, b;
+ unsigned long i;
+ int res;
+
+ MPFR_ASSERTD (!MPFR_IS_SINGULAR (y));
+ MPFR_ASSERTD (!MPFR_IS_SINGULAR (x));
+ MPFR_ASSERTD (!mpfr_integer_p (y));
+ MPFR_ASSERTD (mpfr_cmp_si_2exp (x, MPFR_INT_SIGN (x),
+ MPFR_GET_EXP (x) - 1) != 0);
+ MPFR_ASSERTD (MPFR_IS_POS (x));
+
+ if (MPFR_IS_NEG (y))
+ return 0; /* x is not a power of two => x^-y is not exact */
+
+ /* compute d such that y = c*2^d with c odd integer */
+ mpz_init (c);
+ d = mpfr_get_z_exp (c, y);
+ i = mpz_scan1 (c, 0);
+ mpz_div_2exp (c, c, i);
+ d += i;
+ /* now y=c*2^d with c odd */
+ /* Since y is not an integer, d is necessarily < 0 */
+ MPFR_ASSERTD (d < 0);
+
+ /* Compute a,b such that x=a*2^b */
+ mpz_init (a);
+ b = mpfr_get_z_exp (a, x);
+ i = mpz_scan1 (a, 0);
+ mpz_div_2exp (a, a, i);
+ b += i;
+ /* now x=a*2^b with a is odd */
+
+ for (res = 1 ; d != 0 ; d++)
+ {
+ /* a*2^b is a square iff
+ (i) a is a square when b is even
+ (ii) 2*a is a square when b is odd */
+ if (b % 2 != 0)
+ {
+ mpz_mul_2exp (a, a, 1); /* 2*a */
+ b --;
+ }
+ MPFR_ASSERTD ((b % 2) == 0);
+ if (!mpz_perfect_square_p (a))
+ {
+ res = 0;
+ goto end;
+ }
+ mpz_sqrt (a, a);
+ b = b / 2;
+ }
+ /* Now x = (a'*2^b')^(2^-d) with d < 0
+ so x^y = ((a'*2^b')^(2^-d))^(c*2^d)
+ = ((a'*2^b')^c with c odd integer */
+ {
+ mpfr_t tmp;
+ mp_prec_t p;
+ MPFR_MPZ_SIZEINBASE2 (p, a);
+ mpfr_init2 (tmp, p); /* prec = 1 should not be possible */
+ res = mpfr_set_z (tmp, a, GMP_RNDN);
+ MPFR_ASSERTD (res == 0);
+ res = mpfr_mul_2si (tmp, tmp, b, GMP_RNDN);
+ MPFR_ASSERTD (res == 0);
+ *inexact = mpfr_pow_z (z, tmp, c, rnd_mode);
+ mpfr_clear (tmp);
+ res = 1;
+ }
+ end:
+ mpz_clear (a);
+ mpz_clear (c);
+ return res;
+}
+
+/* Return 1 if y is an odd integer, 0 otherwise. */
+static int
+is_odd (mpfr_srcptr y)
+{
+ mp_exp_t expo;
+ mp_prec_t prec;
+ mp_size_t yn;
+ mp_limb_t *yp;
+
+ /* NAN, INF or ZERO are not allowed */
+ MPFR_ASSERTD (!MPFR_IS_SINGULAR (y));
+
+ expo = MPFR_GET_EXP (y);
+ if (expo <= 0)
+ return 0; /* |y| < 1 and not 0 */
+
+ prec = MPFR_PREC(y);
+ if ((mpfr_prec_t) expo > prec)
+ return 0; /* y is a multiple of 2^(expo-prec), thus not odd */
+
+ /* 0 < expo <= prec:
+ y = 1xxxxxxxxxt.zzzzzzzzzzzzzzzzzz[000]
+ expo bits (prec-expo) bits
+
+ We have to check that:
+ (a) the bit 't' is set
+ (b) all the 'z' bits are zero
+ */
+
+ prec = ((prec - 1) / BITS_PER_MP_LIMB + 1) * BITS_PER_MP_LIMB - expo;
+ /* number of z+0 bits */
+
+ yn = prec / BITS_PER_MP_LIMB;
+ MPFR_ASSERTN(yn >= 0);
+ /* yn is the index of limb containing the 't' bit */
+
+ yp = MPFR_MANT(y);
+ /* if expo is a multiple of BITS_PER_MP_LIMB, t is bit 0 */
+ if (expo % BITS_PER_MP_LIMB == 0 ? (yp[yn] & 1) == 0
+ : yp[yn] << ((expo % BITS_PER_MP_LIMB) - 1) != MPFR_LIMB_HIGHBIT)
+ return 0;
+ while (--yn >= 0)
+ if (yp[yn] != 0)
+ return 0;
+ return 1;
+}
+
+/* Assumes that the exponent range has already been extended and if y is
+ an integer, then the result is not exact in unbounded exponent range. */
+int
+mpfr_pow_general (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y,
+ mp_rnd_t rnd_mode, int y_is_integer, mpfr_save_expo_t *expo)
+{
+ mpfr_t t, u, k, absx;
+ int k_non_zero = 0;
+ int check_exact_case = 0;
+ int inexact;
+ /* Declaration of the size variable */
+ mp_prec_t Nz = MPFR_PREC(z); /* target precision */
+ mp_prec_t Nt; /* working precision */
+ mp_exp_t err; /* error */
+ MPFR_ZIV_DECL (ziv_loop);
+
+
+ MPFR_LOG_FUNC (("x[%#R]=%R y[%#R]=%R rnd=%d", x, x, y, y, rnd_mode),
+ ("z[%#R]=%R inexact=%d", z, z, inexact));
+
+ /* We put the absolute value of x in absx, pointing to the significand
+ of x to avoid allocating memory for the significand of absx. */
+ MPFR_ALIAS(absx, x, /*sign=*/ 1, /*EXP=*/ MPFR_EXP(x));
+
+ /* We will compute the absolute value of the result. So, let's
+ invert the rounding mode if the result is negative. */
+ if (MPFR_IS_NEG (x) && is_odd (y))
+ rnd_mode = MPFR_INVERT_RND (rnd_mode);
+
+ /* compute the precision of intermediary variable */
+ /* the optimal number of bits : see algorithms.tex */
+ Nt = Nz + 5 + MPFR_INT_CEIL_LOG2 (Nz);
+
+ /* initialise of intermediary variable */
+ mpfr_init2 (t, Nt);
+
+ MPFR_ZIV_INIT (ziv_loop, Nt);
+ for (;;)
+ {
+ MPFR_BLOCK_DECL (flags1);
+
+ /* compute exp(y*ln|x|), using GMP_RNDU to get an upper bound, so
+ that we can detect underflows. */
+ mpfr_log (t, absx, MPFR_IS_NEG (y) ? GMP_RNDD : GMP_RNDU); /* ln|x| */
+ mpfr_mul (t, y, t, GMP_RNDU); /* y*ln|x| */
+ if (k_non_zero)
+ {
+ MPFR_LOG_MSG (("subtract k * ln(2)\n", 0));
+ mpfr_const_log2 (u, GMP_RNDD);
+ mpfr_mul (u, u, k, GMP_RNDD);
+ /* Error on u = k * log(2): < k * 2^(-Nt) < 1. */
+ mpfr_sub (t, t, u, GMP_RNDU);
+ MPFR_LOG_MSG (("t = y * ln|x| - k * ln(2)\n", 0));
+ MPFR_LOG_VAR (t);
+ }
+ /* estimate of the error -- see pow function in algorithms.tex.
+ The error on t is at most 1/2 + 3*2^(EXP(t)+1) ulps, which is
+ <= 2^(EXP(t)+3) for EXP(t) >= -1, and <= 2 ulps for EXP(t) <= -2.
+ Additional error if k_no_zero: treal = t * errk, with
+ 1 - |k| * 2^(-Nt) <= exp(-|k| * 2^(-Nt)) <= errk <= 1,
+ i.e., additional absolute error <= 2^(EXP(k)+EXP(t)-Nt).
+ Total error <= 2^err1 + 2^err2 <= 2^(max(err1,err2)+1). */
+ err = MPFR_NOTZERO (t) && MPFR_GET_EXP (t) >= -1 ?
+ MPFR_GET_EXP (t) + 3 : 1;
+ if (k_non_zero)
+ {
+ if (MPFR_GET_EXP (k) > err)
+ err = MPFR_GET_EXP (k);
+ err++;
+ }
+ MPFR_BLOCK (flags1, mpfr_exp (t, t, GMP_RNDN)); /* exp(y*ln|x|)*/
+ /* We need to test */
+ if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (t) || MPFR_UNDERFLOW (flags1)))
+ {
+ mp_prec_t Ntmin;
+ MPFR_BLOCK_DECL (flags2);
+
+ MPFR_ASSERTN (!k_non_zero);
+ MPFR_ASSERTN (!MPFR_IS_NAN (t));
+
+ /* Real underflow? */
+ if (MPFR_IS_ZERO (t))
+ {
+ /* Underflow. We computed rndn(exp(t)), where t >= y*ln|x|.
+ Therefore rndn(|x|^y) = 0, and we have a real underflow on
+ |x|^y. */
+ inexact = mpfr_underflow (z, rnd_mode == GMP_RNDN ? GMP_RNDZ
+ : rnd_mode, MPFR_SIGN_POS);
+ if (expo != NULL)
+ MPFR_SAVE_EXPO_UPDATE_FLAGS (*expo, MPFR_FLAGS_INEXACT
+ | MPFR_FLAGS_UNDERFLOW);
+ break;
+ }
+
+ /* Real overflow? */
+ if (MPFR_IS_INF (t))
+ {
+ /* Note: we can probably use a low precision for this test. */
+ mpfr_log (t, absx, MPFR_IS_NEG (y) ? GMP_RNDU : GMP_RNDD);
+ mpfr_mul (t, y, t, GMP_RNDD); /* y * ln|x| */
+ MPFR_BLOCK (flags2, mpfr_exp (t, t, GMP_RNDD));
+ /* t = lower bound on exp(y * ln|x|) */
+ if (MPFR_OVERFLOW (flags2))
+ {
+ /* We have computed a lower bound on |x|^y, and it
+ overflowed. Therefore we have a real overflow
+ on |x|^y. */
+ inexact = mpfr_overflow (z, rnd_mode, MPFR_SIGN_POS);
+ if (expo != NULL)
+ MPFR_SAVE_EXPO_UPDATE_FLAGS (*expo, MPFR_FLAGS_INEXACT
+ | MPFR_FLAGS_OVERFLOW);
+ break;
+ }
+ }
+
+ k_non_zero = 1;
+ Ntmin = sizeof(mp_exp_t) * CHAR_BIT;
+ if (Ntmin > Nt)
+ {
+ Nt = Ntmin;
+ mpfr_set_prec (t, Nt);
+ }
+ mpfr_init2 (u, Nt);
+ mpfr_init2 (k, Ntmin);
+ mpfr_log2 (k, absx, GMP_RNDN);
+ mpfr_mul (k, y, k, GMP_RNDN);
+ mpfr_round (k, k);
+ MPFR_LOG_VAR (k);
+ /* |y| < 2^Ntmin, therefore |k| < 2^Nt. */
+ continue;
+ }
+ if (MPFR_LIKELY (MPFR_CAN_ROUND (t, Nt - err, Nz, rnd_mode)))
+ {
+ inexact = mpfr_set (z, t, rnd_mode);
+ break;
+ }
+
+ /* check exact power, except when y is an integer (since the
+ exact cases for y integer have already been filtered out) */
+ if (check_exact_case == 0 && ! y_is_integer)
+ {
+ if (mpfr_pow_is_exact (z, absx, y, rnd_mode, &inexact))
+ break;
+ check_exact_case = 1;
+ }
+
+ /* reactualisation of the precision */
+ MPFR_ZIV_NEXT (ziv_loop, Nt);
+ mpfr_set_prec (t, Nt);
+ if (k_non_zero)
+ mpfr_set_prec (u, Nt);
+ }
+ MPFR_ZIV_FREE (ziv_loop);
+
+ if (k_non_zero)
+ {
+ int inex2;
+ long lk;
+
+ /* The rounded result in an unbounded exponent range is z * 2^k. As
+ * MPFR chooses underflow after rounding, the mpfr_mul_2si below will
+ * correctly detect underflows and overflows. However, in rounding to
+ * nearest, if z * 2^k = 2^(emin - 2), then the double rounding may
+ * affect the result. We need to cope with that before overwriting z.
+ * If inexact >= 0, then the real result is <= 2^(emin - 2), so that
+ * o(2^(emin - 2)) = +0 is correct. If inexact < 0, then the real
+ * result is > 2^(emin - 2) and we need to round to 2^(emin - 1).
+ */
+ MPFR_ASSERTN (MPFR_EMAX_MAX <= LONG_MAX);
+ lk = mpfr_get_si (k, GMP_RNDN);
+ if (rnd_mode == GMP_RNDN && inexact < 0 &&
+ MPFR_GET_EXP (z) + lk == __gmpfr_emin - 1 && mpfr_powerof2_raw (z))
+ {
+ /* Rounding to nearest, real result > z * 2^k = 2^(emin - 2),
+ * underflow case: as the minimum precision is > 1, we will
+ * obtain the correct result and exceptions by replacing z by
+ * nextabove(z).
+ */
+ MPFR_ASSERTN (MPFR_PREC_MIN > 1);
+ mpfr_nextabove (z);
+ }
+ mpfr_clear_flags ();
+ inex2 = mpfr_mul_2si (z, z, lk, rnd_mode);
+ if (inex2) /* underflow or overflow */
+ {
+ inexact = inex2;
+ if (expo != NULL)
+ MPFR_SAVE_EXPO_UPDATE_FLAGS (*expo, __gmpfr_flags);
+ }
+ mpfr_clears (u, k, (mpfr_ptr) 0);
+ }
+ mpfr_clear (t);
+
+ /* update the sign of the result if x was negative */
+ if (MPFR_IS_NEG (x) && is_odd (y))
+ {
+ MPFR_SET_NEG(z);
+ inexact = -inexact;
+ }
+
+ return inexact;
+}
+
+/* The computation of z = pow(x,y) is done by
+ z = exp(y * log(x)) = x^y
+ For the special cases, see Section F.9.4.4 of the C standard:
+ _ pow(±0, y) = ±inf for y an odd integer < 0.
+ _ pow(±0, y) = +inf for y < 0 and not an odd integer.
+ _ pow(±0, y) = ±0 for y an odd integer > 0.
+ _ pow(±0, y) = +0 for y > 0 and not an odd integer.
+ _ pow(-1, ±inf) = 1.
+ _ pow(+1, y) = 1 for any y, even a NaN.
+ _ pow(x, ±0) = 1 for any x, even a NaN.
+ _ pow(x, y) = NaN for finite x < 0 and finite non-integer y.
+ _ pow(x, -inf) = +inf for |x| < 1.
+ _ pow(x, -inf) = +0 for |x| > 1.
+ _ pow(x, +inf) = +0 for |x| < 1.
+ _ pow(x, +inf) = +inf for |x| > 1.
+ _ pow(-inf, y) = -0 for y an odd integer < 0.
+ _ pow(-inf, y) = +0 for y < 0 and not an odd integer.
+ _ pow(-inf, y) = -inf for y an odd integer > 0.
+ _ pow(-inf, y) = +inf for y > 0 and not an odd integer.
+ _ pow(+inf, y) = +0 for y < 0.
+ _ pow(+inf, y) = +inf for y > 0. */
+int
+mpfr_pow (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd_mode)
+{
+ int inexact;
+ int cmp_x_1;
+ int y_is_integer;
+ MPFR_SAVE_EXPO_DECL (expo);
+
+ MPFR_LOG_FUNC (("x[%#R]=%R y[%#R]=%R rnd=%d", x, x, y, y, rnd_mode),
+ ("z[%#R]=%R inexact=%d", z, z, inexact));
+
+ if (MPFR_ARE_SINGULAR (x, y))
+ {
+ /* pow(x, 0) returns 1 for any x, even a NaN. */
+ if (MPFR_UNLIKELY (MPFR_IS_ZERO (y)))
+ return mpfr_set_ui (z, 1, rnd_mode);
+ else if (MPFR_IS_NAN (x))
+ {
+ MPFR_SET_NAN (z);
+ MPFR_RET_NAN;
+ }
+ else if (MPFR_IS_NAN (y))
+ {
+ /* pow(+1, NaN) returns 1. */
+ if (mpfr_cmp_ui (x, 1) == 0)
+ return mpfr_set_ui (z, 1, rnd_mode);
+ MPFR_SET_NAN (z);
+ MPFR_RET_NAN;
+ }
+ else if (MPFR_IS_INF (y))
+ {
+ if (MPFR_IS_INF (x))
+ {
+ if (MPFR_IS_POS (y))
+ MPFR_SET_INF (z);
+ else
+ MPFR_SET_ZERO (z);
+ MPFR_SET_POS (z);
+ MPFR_RET (0);
+ }
+ else
+ {
+ int cmp;
+ cmp = mpfr_cmpabs (x, __gmpfr_one) * MPFR_INT_SIGN (y);
+ MPFR_SET_POS (z);
+ if (cmp > 0)
+ {
+ /* Return +inf. */
+ MPFR_SET_INF (z);
+ MPFR_RET (0);
+ }
+ else if (cmp < 0)
+ {
+ /* Return +0. */
+ MPFR_SET_ZERO (z);
+ MPFR_RET (0);
+ }
+ else
+ {
+ /* Return 1. */
+ return mpfr_set_ui (z, 1, rnd_mode);
+ }
+ }
+ }
+ else if (MPFR_IS_INF (x))
+ {
+ int negative;
+ /* Determine the sign now, in case y and z are the same object */
+ negative = MPFR_IS_NEG (x) && is_odd (y);
+ if (MPFR_IS_POS (y))
+ MPFR_SET_INF (z);
+ else
+ MPFR_SET_ZERO (z);
+ if (negative)
+ MPFR_SET_NEG (z);
+ else
+ MPFR_SET_POS (z);
+ MPFR_RET (0);
+ }
+ else
+ {
+ int negative;
+ MPFR_ASSERTD (MPFR_IS_ZERO (x));
+ /* Determine the sign now, in case y and z are the same object */
+ negative = MPFR_IS_NEG(x) && is_odd (y);
+ if (MPFR_IS_NEG (y))
+ MPFR_SET_INF (z);
+ else
+ MPFR_SET_ZERO (z);
+ if (negative)
+ MPFR_SET_NEG (z);
+ else
+ MPFR_SET_POS (z);
+ MPFR_RET (0);
+ }
+ }
+
+ /* x^y for x < 0 and y not an integer is not defined */
+ y_is_integer = mpfr_integer_p (y);
+ if (MPFR_IS_NEG (x) && ! y_is_integer)
+ {
+ MPFR_SET_NAN (z);
+ MPFR_RET_NAN;
+ }
+
+ /* now the result cannot be NaN:
+ (1) either x > 0
+ (2) or x < 0 and y is an integer */
+
+ cmp_x_1 = mpfr_cmpabs (x, __gmpfr_one);
+ if (cmp_x_1 == 0)
+ return mpfr_set_si (z, MPFR_IS_NEG (x) && is_odd (y) ? -1 : 1, rnd_mode);
+
+ /* now we have:
+ (1) either x > 0
+ (2) or x < 0 and y is an integer
+ and in addition |x| <> 1.
+ */
+
+ /* detect overflow: an overflow is possible if
+ (a) |x| > 1 and y > 0
+ (b) |x| < 1 and y < 0.
+ FIXME: this assumes 1 is always representable.
+
+ FIXME2: maybe we can test overflow and underflow simultaneously.
+ The idea is the following: first compute an approximation to
+ y * log2|x|, using rounding to nearest. If |x| is not too near from 1,
+ this approximation should be accurate enough, and in most cases enable
+ one to prove that there is no underflow nor overflow.
+ Otherwise, it should enable one to check only underflow or overflow,
+ instead of both cases as in the present case.
+ */
+ if (cmp_x_1 * MPFR_SIGN (y) > 0)
+ {
+ mpfr_t t;
+ int negative, overflow;
+
+ MPFR_SAVE_EXPO_MARK (expo);
+ mpfr_init2 (t, 53);
+ /* we want a lower bound on y*log2|x|:
+ (i) if x > 0, it suffices to round log2(x) towards zero, and
+ to round y*o(log2(x)) towards zero too;
+ (ii) if x < 0, we first compute t = o(-x), with rounding towards 1,
+ and then follow as in case (1). */
+ if (MPFR_SIGN (x) > 0)
+ mpfr_log2 (t, x, GMP_RNDZ);
+ else
+ {
+ mpfr_neg (t, x, (cmp_x_1 > 0) ? GMP_RNDZ : GMP_RNDU);
+ mpfr_log2 (t, t, GMP_RNDZ);
+ }
+ mpfr_mul (t, t, y, GMP_RNDZ);
+ overflow = mpfr_cmp_si (t, __gmpfr_emax) > 0;
+ mpfr_clear (t);
+ MPFR_SAVE_EXPO_FREE (expo);
+ if (overflow)
+ {
+ MPFR_LOG_MSG (("early overflow detection\n", 0));
+ negative = MPFR_SIGN(x) < 0 && is_odd (y);
+ return mpfr_overflow (z, rnd_mode, negative ? -1 : 1);
+ }
+ }
+
+ /* Basic underflow checking. One has:
+ * - if y > 0, |x^y| < 2^(EXP(x) * y);
+ * - if y < 0, |x^y| <= 2^((EXP(x) - 1) * y);
+ * so that one can compute a value ebound such that |x^y| < 2^ebound.
+ * If we have ebound <= emin - 2 (emin - 1 in directed rounding modes),
+ * then there is an underflow and we can decide the return value.
+ */
+ if (MPFR_IS_NEG (y) ? (MPFR_GET_EXP (x) > 1) : (MPFR_GET_EXP (x) < 0))
+ {
+ mpfr_t tmp;
+ mpfr_eexp_t ebound;
+ int inex2;
+
+ /* We must restore the flags. */
+ MPFR_SAVE_EXPO_MARK (expo);
+ mpfr_init2 (tmp, sizeof (mp_exp_t) * CHAR_BIT);
+ inex2 = mpfr_set_exp_t (tmp, MPFR_GET_EXP (x), GMP_RNDN);
+ MPFR_ASSERTN (inex2 == 0);
+ if (MPFR_IS_NEG (y))
+ {
+ inex2 = mpfr_sub_ui (tmp, tmp, 1, GMP_RNDN);
+ MPFR_ASSERTN (inex2 == 0);
+ }
+ mpfr_mul (tmp, tmp, y, GMP_RNDU);
+ if (MPFR_IS_NEG (y))
+ mpfr_nextabove (tmp);
+ /* tmp doesn't necessarily fit in ebound, but that doesn't matter
+ since we get the minimum value in such a case. */
+ ebound = mpfr_get_exp_t (tmp, GMP_RNDU);
+ mpfr_clear (tmp);
+ MPFR_SAVE_EXPO_FREE (expo);
+ if (MPFR_UNLIKELY (ebound <=
+ __gmpfr_emin - (rnd_mode == GMP_RNDN ? 2 : 1)))
+ {
+ /* warning: mpfr_underflow rounds away from 0 for GMP_RNDN */
+ MPFR_LOG_MSG (("early underflow detection\n", 0));
+ return mpfr_underflow (z,
+ rnd_mode == GMP_RNDN ? GMP_RNDZ : rnd_mode,
+ MPFR_SIGN (x) < 0 && is_odd (y) ? -1 : 1);
+ }
+ }
+
+ /* If y is an integer, we can use mpfr_pow_z (based on multiplications),
+ but if y is very large (I'm not sure about the best threshold -- VL),
+ we shouldn't use it, as it can be very slow and take a lot of memory
+ (and even crash or make other programs crash, as several hundred of
+ MBs may be necessary). Note that in such a case, either x = +/-2^b
+ (this case is handled below) or x^y cannot be represented exactly in
+ any precision supported by MPFR (the general case uses this property).
+ */
+ if (y_is_integer && (MPFR_GET_EXP (y) <= 256))
+ {
+ mpz_t zi;
+
+ MPFR_LOG_MSG (("special code for y not too large integer\n", 0));
+ mpz_init (zi);
+ mpfr_get_z (zi, y, GMP_RNDN);
+ inexact = mpfr_pow_z (z, x, zi, rnd_mode);
+ mpz_clear (zi);
+ return inexact;
+ }
+
+ /* Special case (+/-2^b)^Y which could be exact. If x is negative, then
+ necessarily y is a large integer. */
+ {
+ mp_exp_t b = MPFR_GET_EXP (x) - 1;
+
+ MPFR_ASSERTN (b >= LONG_MIN && b <= LONG_MAX); /* FIXME... */
+ if (mpfr_cmp_si_2exp (x, MPFR_SIGN(x), b) == 0)
+ {
+ mpfr_t tmp;
+ int sgnx = MPFR_SIGN (x);
+
+ MPFR_LOG_MSG (("special case (+/-2^b)^Y\n", 0));
+ /* now x = +/-2^b, so x^y = (+/-1)^y*2^(b*y) is exact whenever b*y is
+ an integer */
+ MPFR_SAVE_EXPO_MARK (expo);
+ mpfr_init2 (tmp, MPFR_PREC (y) + sizeof (long) * CHAR_BIT);
+ inexact = mpfr_mul_si (tmp, y, b, GMP_RNDN); /* exact */
+ MPFR_ASSERTN (inexact == 0);
+ /* Note: as the exponent range has been extended, an overflow is not
+ possible (due to basic overflow and underflow checking above, as
+ the result is ~ 2^tmp), and an underflow is not possible either
+ because b is an integer (thus either 0 or >= 1). */
+ mpfr_clear_flags ();
+ inexact = mpfr_exp2 (z, tmp, rnd_mode);
+ mpfr_clear (tmp);
+ if (sgnx < 0 && is_odd (y))
+ {
+ mpfr_neg (z, z, rnd_mode);
+ inexact = -inexact;
+ }
+ /* Without the following, the overflows3 test in tpow.c fails. */
+ MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);
+ MPFR_SAVE_EXPO_FREE (expo);
+ return mpfr_check_range (z, inexact, rnd_mode);
+ }
+ }
+
+ MPFR_SAVE_EXPO_MARK (expo);
+
+ /* Case where |y * log(x)| is very small. Warning: x can be negative, in
+ that case y is a large integer. */
+ {
+ mpfr_t t;
+ mp_exp_t err;
+
+ /* We need an upper bound on the exponent of y * log(x). */
+ mpfr_init2 (t, 16);
+ if (MPFR_IS_POS(x))
+ mpfr_log (t, x, cmp_x_1 < 0 ? GMP_RNDD : GMP_RNDU); /* away from 0 */
+ else
+ {
+ /* if x < -1, round to +Inf, else round to zero */
+ mpfr_neg (t, x, (mpfr_cmp_si (x, -1) < 0) ? GMP_RNDU : GMP_RNDD);
+ mpfr_log (t, t, (mpfr_cmp_ui (t, 1) < 0) ? GMP_RNDD : GMP_RNDU);
+ }
+ MPFR_ASSERTN (MPFR_IS_PURE_FP (t));
+ err = MPFR_GET_EXP (y) + MPFR_GET_EXP (t);
+ mpfr_clear (t);
+ mpfr_clear_flags ();
+ MPFR_SMALL_INPUT_AFTER_SAVE_EXPO (z, __gmpfr_one, - err, 0,
+ (MPFR_SIGN (y) > 0) ^ (cmp_x_1 < 0),
+ rnd_mode, expo, {});
+ }
+
+ /* General case */
+ inexact = mpfr_pow_general (z, x, y, rnd_mode, y_is_integer, &expo);
+
+ MPFR_SAVE_EXPO_FREE (expo);
+ return mpfr_check_range (z, inexact, rnd_mode);
+}
projects / msp430-gcc.git / blobdiff