[gcc] GCC 4.5.0=>4.5.1

gcc/mpfr/pow_z.c

Index: gcc/mpfr/pow_z.c
diff --git a/gcc/mpfr/pow_z.c b/gcc/mpfr/pow_z.c
deleted file mode 100644
index 3e88a9bf2533a76dcce33020ef94c1443f5e2c9b..0000000000000000000000000000000000000000
--- a/gcc/mpfr/pow_z.c
+++ /dev/null
@@ -1,365 +0,0 @@
-/* mpfr_pow_z -- power function x^z with z a MPZ
-
-Copyright 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"
-
-/* y <- x^|z| with z != 0
-   if cr=1: ensures correct rounding of y
-   if cr=0: does not ensure correct rounding, but avoid spurious overflow
-   or underflow, and uses the precision of y as working precision (warning,
-   y and x might be the same variable). */
-static int
-mpfr_pow_pos_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mp_rnd_t rnd, int cr)
-{
-  mpfr_t res;
-  mp_prec_t prec, err;
-  int inexact;
-  mp_rnd_t rnd1, rnd2;
-  mpz_t absz;
-  mp_size_t size_z;
-  MPFR_ZIV_DECL (loop);
-  MPFR_BLOCK_DECL (flags);
-
-  MPFR_LOG_FUNC (("x[%#R]=%R z=? rnd=%d cr=%d", x, x, rnd, cr),
-                 ("y[%#R]=%R inexact=%d", y, y, inexact));
-
-  MPFR_ASSERTD (mpz_sgn (z) != 0);
-
-  if (MPFR_UNLIKELY (mpz_cmpabs_ui (z, 1) == 0))
-    return mpfr_set (y, x, rnd);
-
-  absz[0] = z[0];
-  SIZ (absz) = ABS(SIZ(absz)); /* Hack to get abs(z) */
-  MPFR_MPZ_SIZEINBASE2 (size_z, z);
-
-  /* round towards 1 (or -1) to avoid spurious overflow or underflow,
-     i.e. if an overflow or underflow occurs, it is a real exception
-     and is not just due to the rounding error. */
-  rnd1 = (MPFR_EXP(x) >= 1) ? GMP_RNDZ
-    : (MPFR_IS_POS(x) ? GMP_RNDU : GMP_RNDD);
-  rnd2 = (MPFR_EXP(x) >= 1) ? GMP_RNDD : GMP_RNDU;
-
-  if (cr != 0)
-    prec = MPFR_PREC (y) + 3 + size_z + MPFR_INT_CEIL_LOG2 (MPFR_PREC (y));
-  else
-    prec = MPFR_PREC (y);
-  mpfr_init2 (res, prec);
-
-  MPFR_ZIV_INIT (loop, prec);
-  for (;;)
-    {
-      unsigned int inexmul;  /* will be non-zero if res may be inexact */
-      mp_size_t i = size_z;
-
-      /* now 2^(i-1) <= z < 2^i */
-      /* see below (case z < 0) for the error analysis, which is identical,
-         except if z=n, the maximal relative error is here 2(n-1)2^(-prec)
-         instead of 2(2n-1)2^(-prec) for z<0. */
-      MPFR_ASSERTD (prec > (mpfr_prec_t) i);
-      err = prec - 1 - (mpfr_prec_t) i;
-
-      MPFR_BLOCK (flags,
-                  inexmul = mpfr_mul (res, x, x, rnd2);
-                  MPFR_ASSERTD (i >= 2);
-                  if (mpz_tstbit (absz, i - 2))
-                    inexmul |= mpfr_mul (res, res, x, rnd1);
-                  for (i -= 3; i >= 0 && !MPFR_BLOCK_EXCEP; i--)
-                    {
-                      inexmul |= mpfr_mul (res, res, res, rnd2);
-                      if (mpz_tstbit (absz, i))
-                        inexmul |= mpfr_mul (res, res, x, rnd1);
-                    });
-      if (MPFR_LIKELY (inexmul == 0 || cr == 0
-                       || MPFR_OVERFLOW (flags) || MPFR_UNDERFLOW (flags)
-                       || MPFR_CAN_ROUND (res, err, MPFR_PREC (y), rnd)))
-        break;
-      /* Can't decide correct rounding, increase the precision */
-      MPFR_ZIV_NEXT (loop, prec);
-      mpfr_set_prec (res, prec);
-    }
-  MPFR_ZIV_FREE (loop);
-
-  /* Check Overflow */
-  if (MPFR_OVERFLOW (flags))
-    {
-      MPFR_LOG_MSG (("overflow\n", 0));
-      inexact = mpfr_overflow (y, rnd, mpz_odd_p (absz) ?
-                               MPFR_SIGN (x) : MPFR_SIGN_POS);
-    }
-  /* Check Underflow */
-  else if (MPFR_UNDERFLOW (flags))
-    {
-      MPFR_LOG_MSG (("underflow\n", 0));
-      if (rnd == GMP_RNDN)
-        {
-          mpfr_t y2, zz;
-
-          /* We cannot decide now whether the result should be rounded
-             toward zero or +Inf. So, let's use the general case of
-             mpfr_pow, which can do that. But the problem is that the
-             result can be exact! However, it is sufficient to try to
-             round on 2 bits (the precision does not matter in case of
-             underflow, since MPFR does not have subnormals), in which
-             case, the result cannot be exact due to previous filtering
-             of trivial cases. */
-          MPFR_ASSERTD (mpfr_cmp_si_2exp (x, MPFR_SIGN (x),
-                                          MPFR_EXP (x) - 1) != 0);
-          mpfr_init2 (y2, 2);
-          mpfr_init2 (zz, ABS (SIZ (z)) * BITS_PER_MP_LIMB);
-          inexact = mpfr_set_z (zz, z, GMP_RNDN);
-          MPFR_ASSERTN (inexact == 0);
-          inexact = mpfr_pow_general (y2, x, zz, rnd, 1,
-                                      (mpfr_save_expo_t *) NULL);
-          mpfr_clear (zz);
-          mpfr_set (y, y2, GMP_RNDN);
-          mpfr_clear (y2);
-          __gmpfr_flags = MPFR_FLAGS_INEXACT | MPFR_FLAGS_UNDERFLOW;
-        }
-      else
-        {
-          inexact = mpfr_underflow (y, rnd, mpz_odd_p (absz) ?
-                                    MPFR_SIGN (x) : MPFR_SIGN_POS);
-        }
-    }
-  else
-    inexact = mpfr_set (y, res, rnd);
-
-  mpfr_clear (res);
-  return inexact;
-}
-
-/* The computation of y = pow(x,z) is done by
- *    y = set_ui(1)      if z = 0
- *    y = pow_ui(x,z)    if z > 0
- *    y = pow_ui(1/x,-z) if z < 0
- *
- * Note: in case z < 0, we could also compute 1/pow_ui(x,-z). However, in
- * case MAX < 1/MIN, where MAX is the largest positive value, i.e.,
- * MAX = nextbelow(+Inf), and MIN is the smallest positive value, i.e.,
- * MIN = nextabove(+0), then x^(-z) might produce an overflow, whereas
- * x^z is representable.
- */
-
-int
-mpfr_pow_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mp_rnd_t rnd)
-{
-  int   inexact;
-  mpz_t tmp;
-  MPFR_SAVE_EXPO_DECL (expo);
-
-  MPFR_LOG_FUNC (("x[%#R]=%R z=? rnd=%d", x, x, rnd),
-                 ("y[%#R]=%R inexact=%d", y, y, inexact));
-
-  /* x^0 = 1 for any x, even a NaN */
-  if (MPFR_UNLIKELY (mpz_sgn (z) == 0))
-    return mpfr_set_ui (y, 1, rnd);
-
-  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
-    {
-      if (MPFR_IS_NAN (x))
-        {
-          MPFR_SET_NAN (y);
-          MPFR_RET_NAN;
-        }
-      else if (MPFR_IS_INF (x))
-        {
-          /* Inf^n = Inf, (-Inf)^n = Inf for n even, -Inf for n odd */
-          /* Inf ^(-n) = 0, sign = + if x>0 or z even */
-          if (mpz_sgn (z) > 0)
-            MPFR_SET_INF (y);
-          else
-            MPFR_SET_ZERO (y);
-          if (MPFR_UNLIKELY (MPFR_IS_NEG (x) && mpz_odd_p (z)))
-            MPFR_SET_NEG (y);
-          else
-            MPFR_SET_POS (y);
-          MPFR_RET (0);
-        }
-      else /* x is zero */
-        {
-          MPFR_ASSERTD (MPFR_IS_ZERO(x));
-          if (mpz_sgn (z) > 0)
-            /* 0^n = +/-0 for any n */
-            MPFR_SET_ZERO (y);
-          else
-            /* 0^(-n) if +/- INF */
-            MPFR_SET_INF (y);
-          if (MPFR_LIKELY (MPFR_IS_POS (x) || mpz_even_p (z)))
-            MPFR_SET_POS (y);
-          else
-            MPFR_SET_NEG (y);
-          MPFR_RET(0);
-        }
-    }
-
-  /* detect exact powers: x^-n is exact iff x is a power of 2
-     Do it if n > 0 too as this is faster and this filtering is
-     needed in case of underflow. */
-  if (MPFR_UNLIKELY (mpfr_cmp_si_2exp (x, MPFR_SIGN (x),
-                                       MPFR_EXP (x) - 1) == 0))
-    {
-      mp_exp_t expx = MPFR_EXP (x); /* warning: x and y may be the same
-                                       variable */
-
-      MPFR_LOG_MSG (("x^n with x power of two\n", 0));
-      mpfr_set_si (y, mpz_odd_p (z) ? MPFR_INT_SIGN(x) : 1, rnd);
-      MPFR_ASSERTD (MPFR_IS_FP (y));
-      mpz_init (tmp);
-      mpz_mul_si (tmp, z, expx - 1);
-      MPFR_ASSERTD (MPFR_GET_EXP (y) == 1);
-      mpz_add_ui (tmp, tmp, 1);
-      inexact = 0;
-      if (MPFR_UNLIKELY (mpz_cmp_si (tmp, __gmpfr_emin) < 0))
-        {
-          MPFR_LOG_MSG (("underflow\n", 0));
-          /* |y| is a power of two, thus |y| <= 2^(emin-2), and in
-             rounding to nearest, the value must be rounded to 0. */
-          if (rnd == GMP_RNDN)
-            rnd = GMP_RNDZ;
-          inexact = mpfr_underflow (y, rnd, MPFR_SIGN (y));
-        }
-      else if (MPFR_UNLIKELY (mpz_cmp_si (tmp, __gmpfr_emax) > 0))
-        {
-          MPFR_LOG_MSG (("overflow\n", 0));
-          inexact = mpfr_overflow (y, rnd, MPFR_SIGN (y));
-        }
-      else
-        MPFR_SET_EXP (y, mpz_get_si (tmp));
-      mpz_clear (tmp);
-      MPFR_RET (inexact);
-    }
-
-  MPFR_SAVE_EXPO_MARK (expo);
-
-  if (mpz_sgn (z) > 0)
-    {
-      inexact = mpfr_pow_pos_z (y, x, z, rnd, 1);
-      MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);
-    }
-  else
-    {
-      /* Declaration of the intermediary variable */
-      mpfr_t t;
-      mp_prec_t Nt;   /* Precision of the intermediary variable */
-      mp_rnd_t rnd1;
-      mp_size_t size_z;
-      MPFR_ZIV_DECL (loop);
-
-      MPFR_MPZ_SIZEINBASE2 (size_z, z);
-
-      /* initial working precision */
-      Nt = MPFR_PREC (y);
-      Nt = Nt + size_z + 3 + MPFR_INT_CEIL_LOG2 (Nt);
-      /* ensures Nt >= bits(z)+2 */
-
-      /* initialise of intermediary variable */
-      mpfr_init2 (t, Nt);
-
-      /* We will compute rnd(rnd1(1/x) ^ (-z)), where rnd1 is the rounding
-         toward sign(x), to avoid spurious overflow or underflow. */
-      rnd1 = MPFR_EXP (x) < 1 ? GMP_RNDZ :
-        (MPFR_SIGN (x) > 0 ? GMP_RNDU : GMP_RNDD);
-
-      MPFR_ZIV_INIT (loop, Nt);
-      for (;;)
-        {
-          MPFR_BLOCK_DECL (flags);
-
-          /* compute (1/x)^(-z), -z>0 */
-          /* As emin = -emax, an underflow cannot occur in the division.
-             And if an overflow occurs, then this means that x^z overflows
-             too (since we have rounded toward 1 or -1). */
-          MPFR_BLOCK (flags, mpfr_ui_div (t, 1, x, rnd1));
-          MPFR_ASSERTD (! MPFR_UNDERFLOW (flags));
-          /* t = (1/x)*(1+theta) where |theta| <= 2^(-Nt) */
-          if (MPFR_UNLIKELY (MPFR_OVERFLOW (flags)))
-            goto overflow;
-          MPFR_BLOCK (flags, mpfr_pow_pos_z (t, t, z, rnd, 0));
-          /* Now if z=-n, t = x^z*(1+theta)^(2n-1) where |theta| <= 2^(-Nt),
-             with theta maybe different from above. If (2n-1)*2^(-Nt) <= 1/2,
-             which is satisfied as soon as Nt >= bits(z)+2, then we can use
-             Lemma \ref{lemma_graillat} from algorithms.tex, which yields
-             t = x^z*(1+theta) with |theta| <= 2(2n-1)*2^(-Nt), thus the
-             error is bounded by 2(2n-1) ulps <= 2^(bits(z)+2) ulps. */
-          if (MPFR_UNLIKELY (MPFR_OVERFLOW (flags)))
-            {
-            overflow:
-              MPFR_ZIV_FREE (loop);
-              mpfr_clear (t);
-              MPFR_SAVE_EXPO_FREE (expo);
-              MPFR_LOG_MSG (("overflow\n", 0));
-              return mpfr_overflow (y, rnd,
-                                    mpz_odd_p (z) ? MPFR_SIGN (x) :
-                                    MPFR_SIGN_POS);
-            }
-          if (MPFR_UNLIKELY (MPFR_UNDERFLOW (flags)))
-            {
-              MPFR_ZIV_FREE (loop);
-              mpfr_clear (t);
-              MPFR_LOG_MSG (("underflow\n", 0));
-              if (rnd == GMP_RNDN)
-                {
-                  mpfr_t y2, zz;
-
-                  /* We cannot decide now whether the result should be
-                     rounded toward zero or away from zero. So, like
-                     in mpfr_pow_pos_z, let's use the general case of
-                     mpfr_pow in precision 2. */
-                  MPFR_ASSERTD (mpfr_cmp_si_2exp (x, MPFR_SIGN (x),
-                                                  MPFR_EXP (x) - 1) != 0);
-                  mpfr_init2 (y2, 2);
-                  mpfr_init2 (zz, ABS (SIZ (z)) * BITS_PER_MP_LIMB);
-                  inexact = mpfr_set_z (zz, z, GMP_RNDN);
-                  MPFR_ASSERTN (inexact == 0);
-                  inexact = mpfr_pow_general (y2, x, zz, rnd, 1,
-                                              (mpfr_save_expo_t *) NULL);
-                  mpfr_clear (zz);
-                  mpfr_set (y, y2, GMP_RNDN);
-                  mpfr_clear (y2);
-                  MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, MPFR_FLAGS_UNDERFLOW);
-                  goto end;
-                }
-              else
-                {
-                  MPFR_SAVE_EXPO_FREE (expo);
-                  return mpfr_underflow (y, rnd, mpz_odd_p (z) ?
-                                         MPFR_SIGN (x) : MPFR_SIGN_POS);
-                }
-            }
-          if (MPFR_LIKELY (MPFR_CAN_ROUND (t, Nt - size_z - 2, MPFR_PREC (y),
-                                           rnd)))
-            break;
-          /* actualisation of the precision */
-          MPFR_ZIV_NEXT (loop, Nt);
-          mpfr_set_prec (t, Nt);
-        }
-      MPFR_ZIV_FREE (loop);
-
-      inexact = mpfr_set (y, t, rnd);
-      mpfr_clear (t);
-    }
-
- end:
-  MPFR_SAVE_EXPO_FREE (expo);
-  return mpfr_check_range (y, inexact, rnd);
-}