Merged vlefevre branch:

    svn merge -c-5445 .
    svn merge -r5436:HEAD .../mpfr/branches/vlefevre
* pow.c:
  - Moved the general case from mpfr_pow() to a new internal function
    mpfr_pow_general().
  - In this function (from old code), avoid unnecessary overflow test
    if the intermediate result is not an infinity (which was the case
    of underflow with non-zero result, thus not an overflow).
  - Fixed a double-rounding problem that occurred in this function in
    some underflow cases when rescaling the result.
  - Added log messages.
* mpfr-impl.h: added mpfr_pow_general prototype.
* pow_z.c:
  - The underflow case of mpfr_pow_pos_z() in rounding to nearest,
    which was incorrect, is now handled by calling mpfr_pow_general(),
    which can scale the result thus decide whether the rounded result
    should be 0 or nextabove(0). To avoid the exact cases of x^y with
    y integer (not supported by mpfr_pow_general()), rounding is done
    in precision 2 (this is also faster!).
  - Fixed underflow-related bug (case exact result = 2^(emin-2), in
    rounding to nearest).
  - Added log messages.
* pow_ui.c:
  - Swapped parameters x and y for consistency (-> y = x^n).
  - Fixed the internal overflows and underflows (which could yield
    spurious overflows/underflows and incorrect results) by using
    mpfr_pow_z.
* tests/tpow_all.c:
  - Test flags in test_others and cmpres; cmpres argument z1 can now
    be a null pointer (if unknown pure FP value, thus not tested).
  - Added a test of 2^(emin - i/4) with 0 <= i <= 12, that triggered
    the bugs mentioned above (and now fixed).

git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@5505 280ebfd0-de03-0410-8827-d642c229c3f4

5 files changed, 474 insertions, 258 deletions

diff --git a/mpfr-impl.h b/mpfr-impl.h
index a1f974c69..a23a4f4ff 100644
--- a/mpfr-impl.h
+++ b/mpfr-impl.h
@@ -1565,6 +1565,9 @@ __MPFR_DECLSPEC int mpfr_exp_2 _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr,mp_rnd_t));
 __MPFR_DECLSPEC int mpfr_exp_3 _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr,mp_rnd_t));
 __MPFR_DECLSPEC int mpfr_powerof2_raw _MPFR_PROTO ((mpfr_srcptr));
 
+__MPFR_DECLSPEC int mpfr_pow_general _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr,
+                           mpfr_srcptr, mp_rnd_t, int, mpfr_save_expo_t *));
+
 __MPFR_DECLSPEC void mpfr_setmax _MPFR_PROTO ((mpfr_ptr, mp_exp_t));
 __MPFR_DECLSPEC void mpfr_setmin _MPFR_PROTO ((mpfr_ptr, mp_exp_t));
 
diff --git a/pow.c b/pow.c
index 2ff46a9c0..ddd06c91d 100644
--- a/pow.c
+++ b/pow.c
@@ -154,6 +154,206 @@ is_odd (mpfr_srcptr y)
   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, exp_te;                      /* 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, GMP_RNDU);            /* ln|x| */
+      mpfr_mul (t, y, t, GMP_RNDU);            /* y*ln|x| */
+      if (k_non_zero)
+        {
+          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);
+        }
+      exp_te = MPFR_GET_EXP (t);               /* FIXME: May overflow */
+      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, GMP_RNDD);            /* ln|x| */
+              mpfr_mul (t, y, t, GMP_RNDD);            /* y * ln|x| */
+              MPFR_BLOCK (flags2, mpfr_exp (t, t, GMP_RNDD));
+              /* t = 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;
+        }
+      /* estimate of the error -- see pow function in algorithms.tex.
+         The error on t is at most 1/2 + 3*2^(exp_te+1) ulps, which is
+         <= 2^(exp_te+3) for exp_te >= -1, and <= 2 ulps for exp_te <= -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 = exp_te >= -1 ? exp_te + 3 : 1;
+      if (k_non_zero)
+        {
+          if (MPFR_GET_EXP (k) > err)
+            err = MPFR_GET_EXP (k);
+          err++;
+        }
+      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:
@@ -332,6 +532,7 @@ mpfr_pow (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd_mode)
       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);
         }
@@ -355,6 +556,7 @@ mpfr_pow (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd_mode)
       if (underflow)
         {
           /* warning: mpfr_underflow rounds away from 0 for GMP_RNDN */
+          MPFR_LOG_MSG (("early underflow detection\n", 0));
           negative = MPFR_SIGN(x) < 0 && is_odd (y);
           return mpfr_underflow (z, (rnd_mode == GMP_RNDN) ? GMP_RNDZ :
                                  rnd_mode, negative ? -1 : 1);
@@ -373,6 +575,7 @@ mpfr_pow (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd_mode)
     {
       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);
@@ -391,6 +594,7 @@ mpfr_pow (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd_mode)
         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);
@@ -444,161 +648,7 @@ mpfr_pow (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd_mode)
   }
 
   /* General case */
-  {
-    /* Declaration of the intermediary variable */
-    mpfr_t t, u, k, absx;
-    int k_non_zero = 0;
-    int check_exact_case = 0;
-    /* Declaration of the size variable */
-    mp_prec_t Nz = MPFR_PREC(z);               /* target precision */
-    mp_prec_t Nt;                              /* working precision */
-    mp_exp_t err, exp_te;                      /* error */
-    MPFR_ZIV_DECL (ziv_loop);
-
-    /* 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, GMP_RNDU);            /* ln|x| */
-        mpfr_mul (t, y, t, GMP_RNDU);            /* y*ln|x| */
-        if (k_non_zero)
-          {
-            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);
-          }
-        exp_te = MPFR_GET_EXP (t);               /* FIXME: May overflow */
-        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));
-            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);
-                MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, MPFR_FLAGS_INEXACT
-                                             | MPFR_FLAGS_UNDERFLOW);
-                break;
-              }
-
-            /* Overflow. */
-            /* Note: we can probably use a low precision for this test. */
-            mpfr_log (t, absx, GMP_RNDD);            /* ln|x| */
-            mpfr_mul (t, y, t, GMP_RNDD);            /* y*ln|x| */
-            MPFR_BLOCK (flags2, mpfr_exp (t, t, GMP_RNDD));  /* 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);
-                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);
-            /* |y| < 2^Ntmin, therefore |k| < 2^Nt. */
-            continue;
-          }
-        /* estimate of the error -- see pow function in algorithms.tex.
-           The error on t is at most 1/2 + 3*2^(exp_te+1) ulps, which is
-           <= 2^(exp_te+3) for exp_te >= -1, and <= 2 ulps for exp_te <= -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 = exp_te >= -1 ? exp_te + 3 : 1;
-        if (k_non_zero)
-          {
-            if (MPFR_GET_EXP (k) > err)
-              err = MPFR_GET_EXP (k);
-            err++;
-          }
-        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;
-
-        MPFR_ASSERTN (MPFR_EMAX_MAX <= LONG_MAX);
-        mpfr_clear_flags ();
-        inex2 = mpfr_mul_2si (z, z, mpfr_get_si (k, GMP_RNDN), rnd_mode);
-        if (inex2)  /* underflow or overflow */
-          {
-            inexact = inex2;
-            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;
-    }
+  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);
diff --git a/pow_ui.c b/pow_ui.c
index 4a3d04b62..c12a3189f 100644
--- a/pow_ui.c
+++ b/pow_ui.c
@@ -24,9 +24,9 @@ MA 02110-1301, USA. */
 #define MPFR_NEED_LONGLONG_H
 #include "mpfr-impl.h"
 
-/* sets x to y^n, and return 0 if exact, non-zero otherwise */
+/* sets y to x^n, and return 0 if exact, non-zero otherwise */
 int
-mpfr_pow_ui (mpfr_ptr x, mpfr_srcptr y, unsigned long int n, mp_rnd_t rnd)
+mpfr_pow_ui (mpfr_ptr y, mpfr_srcptr x, unsigned long int n, mp_rnd_t rnd)
 {
   unsigned long m;
   mpfr_t res;
@@ -37,58 +37,61 @@ mpfr_pow_ui (mpfr_ptr x, mpfr_srcptr y, unsigned long int n, mp_rnd_t rnd)
   MPFR_ZIV_DECL (loop);
   MPFR_BLOCK_DECL (flags);
 
-  /* y^0 = 1 for any y, even a NaN */
+  MPFR_LOG_FUNC (("x[%#R]=%R n=%lu rnd=%d", x, x, n, rnd),
+                 ("y[%#R]=%R inexact=%d", y, y, inexact));
+
+  /* x^0 = 1 for any x, even a NaN */
   if (MPFR_UNLIKELY (n == 0))
-    return mpfr_set_ui (x, 1, rnd);
+    return mpfr_set_ui (y, 1, rnd);
 
-  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (y)))
+  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
     {
-      if (MPFR_IS_NAN (y))
+      if (MPFR_IS_NAN (x))
         {
-          MPFR_SET_NAN (x);
+          MPFR_SET_NAN (y);
           MPFR_RET_NAN;
         }
-      else if (MPFR_IS_INF (y))
+      else if (MPFR_IS_INF (x))
         {
           /* Inf^n = Inf, (-Inf)^n = Inf for n even, -Inf for n odd */
-          if ((MPFR_IS_NEG (y)) && ((n & 1) == 1))
-            MPFR_SET_NEG (x);
+          if (MPFR_IS_NEG (x) && (n & 1) == 1)
+            MPFR_SET_NEG (y);
           else
-            MPFR_SET_POS (x);
-          MPFR_SET_INF (x);
+            MPFR_SET_POS (y);
+          MPFR_SET_INF (y);
           MPFR_RET (0);
         }
-      else /* y is zero */
+      else /* x is zero */
         {
-          MPFR_ASSERTD (MPFR_IS_ZERO (y));
+          MPFR_ASSERTD (MPFR_IS_ZERO (x));
           /* 0^n = 0 for any n */
-          MPFR_SET_ZERO (x);
-          if (MPFR_IS_POS (y) || ((n & 1) == 0))
-            MPFR_SET_POS (x);
+          MPFR_SET_ZERO (y);
+          if (MPFR_IS_POS (x) || (n & 1) == 0)
+            MPFR_SET_POS (y);
           else
-            MPFR_SET_NEG (x);
+            MPFR_SET_NEG (y);
           MPFR_RET (0);
         }
     }
   else if (MPFR_UNLIKELY (n <= 2))
     {
       if (n < 2)
-        /* y^1 = y */
-        return mpfr_set (x, y, rnd);
+        /* x^1 = x */
+        return mpfr_set (y, x, rnd);
       else
-        /* y^2 = sqr(y) */
-        return mpfr_sqr (x, y, rnd);
+        /* x^2 = sqr(x) */
+        return mpfr_sqr (y, x, rnd);
     }
 
   /* Augment exponent range */
   MPFR_SAVE_EXPO_MARK (expo);
 
   /* setup initial precision */
-  prec = MPFR_PREC (x) + 3 + BITS_PER_MP_LIMB
-    + MPFR_INT_CEIL_LOG2 (MPFR_PREC (x));
+  prec = MPFR_PREC (y) + 3 + BITS_PER_MP_LIMB
+    + MPFR_INT_CEIL_LOG2 (MPFR_PREC (y));
   mpfr_init2 (res, prec);
 
-  rnd1 = MPFR_IS_POS (y) ? GMP_RNDU : GMP_RNDD; /* away */
+  rnd1 = MPFR_IS_POS (x) ? GMP_RNDU : GMP_RNDD; /* away */
 
   MPFR_ZIV_INIT (loop, prec);
   for (;;)
@@ -100,17 +103,17 @@ mpfr_pow_ui (mpfr_ptr x, mpfr_srcptr y, unsigned long int n, mp_rnd_t rnd)
       /* now 2^(i-1) <= n < 2^i */
       MPFR_ASSERTD (prec > (mpfr_prec_t) i);
       err = prec - 1 - (mpfr_prec_t) i;
-      /* First step: compute square from y */
+      /* First step: compute square from x */
       MPFR_BLOCK (flags,
-                  inexact = mpfr_mul (res, y, y, GMP_RNDU);
+                  inexact = mpfr_mul (res, x, x, GMP_RNDU);
                   MPFR_ASSERTD (i >= 2);
                   if (n & (1UL << (i-2)))
-                    inexact |= mpfr_mul (res, res, y, rnd1);
+                    inexact |= mpfr_mul (res, res, x, rnd1);
                   for (i -= 3; i >= 0 && !MPFR_BLOCK_EXCEP; i--)
                     {
                       inexact |= mpfr_mul (res, res, res, GMP_RNDU);
                       if (n & (1UL << i))
-                        inexact |= mpfr_mul (res, res, y, rnd1);
+                        inexact |= mpfr_mul (res, res, x, rnd1);
                     });
       /* let r(n) be the number of roundings: we have r(2)=1, r(3)=2,
          and r(2n)=2r(n)+1, r(2n+1)=2r(n)+2, thus r(n)=n-1.
@@ -122,7 +125,7 @@ mpfr_pow_ui (mpfr_ptr x, mpfr_srcptr y, unsigned long int n, mp_rnd_t rnd)
       */
       if (MPFR_LIKELY (inexact == 0
                        || MPFR_OVERFLOW (flags) || MPFR_UNDERFLOW (flags)
-                       || MPFR_CAN_ROUND (res, err, MPFR_PREC (x), rnd)))
+                       || MPFR_CAN_ROUND (res, err, MPFR_PREC (y), rnd)))
         break;
       /* Actualisation of the precision */
       MPFR_ZIV_NEXT (loop, prec);
@@ -130,26 +133,29 @@ mpfr_pow_ui (mpfr_ptr x, mpfr_srcptr y, unsigned long int n, mp_rnd_t rnd)
     }
   MPFR_ZIV_FREE (loop);
 
-  /* Check Overflow */
-  if (MPFR_OVERFLOW (flags))
-    {
-      mpfr_clear (res);
-      MPFR_SAVE_EXPO_FREE (expo);
-      return mpfr_overflow (x, rnd,
-                            (n % 2) ? MPFR_SIGN (y) : MPFR_SIGN_POS);
-    }
-  /* Check Underflow  */
-  else if (MPFR_UNDERFLOW (flags))
+  if (MPFR_UNLIKELY (MPFR_OVERFLOW (flags) || MPFR_UNDERFLOW (flags)))
     {
+      mpz_t z;
+
+      /* Internal overflow or underflow. However the approximation error has
+       * not been taken into account. So, let's solve this problem by using
+       * mpfr_pow_z, which can handle it. This case could be improved in the
+       * future, without having to use mpfr_pow_z.
+       */
+      MPFR_LOG_MSG (("Internal overflow or underflow,"
+                     " let's use mpfr_pow_z.\n", 0));
       mpfr_clear (res);
       MPFR_SAVE_EXPO_FREE (expo);
-      return mpfr_underflow (x, rnd == GMP_RNDN ? GMP_RNDZ : rnd,
-                             (n % 2) ? MPFR_SIGN(y) : MPFR_SIGN_POS);
+      mpz_init (z);
+      mpz_set_ui (z, n);
+      inexact = mpfr_pow_z (y, x, z, rnd);
+      mpz_clear (z);
+      return inexact;
     }
 
-  inexact = mpfr_set (x, res, rnd);
+  inexact = mpfr_set (y, res, rnd);
   mpfr_clear (res);
 
   MPFR_SAVE_EXPO_FREE (expo);
-  return mpfr_check_range (x, inexact, rnd);
+  return mpfr_check_range (y, inexact, rnd);
 }
diff --git a/pow_z.c b/pow_z.c
index b2b4ec6a4..4a364e400 100644
--- a/pow_z.c
+++ b/pow_z.c
@@ -39,6 +39,9 @@ mpfr_pow_pos_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mp_rnd_t rnd, int cr)
   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))
@@ -97,13 +100,46 @@ mpfr_pow_pos_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mp_rnd_t rnd, int cr)
 
   /* Check Overflow */
   if (MPFR_OVERFLOW (flags))
-    inexact = mpfr_overflow (y, rnd, mpz_odd_p (absz) ?
-                             MPFR_SIGN (x) : MPFR_SIGN_POS);
+    {
+      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))
-    inexact = mpfr_underflow (y, rnd == GMP_RNDN ? GMP_RNDZ : rnd,
-                              mpz_odd_p (absz) ? MPFR_SIGN (x) :
-                              MPFR_SIGN_POS);
+    {
+      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
+        {
+          MPFR_LOG_MSG (("underflow\n", 0));
+          inexact = mpfr_underflow (y, rnd, mpz_odd_p (absz) ?
+                                    MPFR_SIGN (x) : MPFR_SIGN_POS);
+        }
+    }
   else
     inexact = mpfr_set (y, res, rnd);
 
@@ -130,6 +166,9 @@ mpfr_pow_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mp_rnd_t rnd)
   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);
@@ -173,35 +212,36 @@ mpfr_pow_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mp_rnd_t rnd)
     }
 
   /* detect exact powers: x^-n is exact iff x is a power of 2
-     Do it if n > 0 too (faster). */
+     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);
+      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))
         {
-          /* The following test is necessary because in the rounding to the
-           * nearest mode, mpfr_underflow always rounds away from 0. In
-           * this rounding mode, we need to round to 0 if:
-           *   _ |y| < 2^(emin-2), or
-           *   _ |y| = 2^(emin-2) and the absolute value of the exact
-           *     result is <= 2^(emin-2).
-           * NOTE: y is a power of 2 and inexact = 0!
-           */
-          if (rnd == GMP_RNDN && mpz_cmp_si (tmp, __gmpfr_emin-1) < 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))
-        inexact = mpfr_overflow (y, rnd, MPFR_SIGN (y));
+        {
+          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);
diff --git a/tests/tpow_all.c b/tests/tpow_all.c
index 792334bf1..d1a6ba153 100644
--- a/tests/tpow_all.c
+++ b/tests/tpow_all.c
@@ -58,15 +58,37 @@ err (const char *s, int i, int j, int rnd, mpfr_srcptr z, int inex)
   exit (1);
 }
 
+/* Arguments:
+ *   spx: non-zero if px is a stringm zero if px is a MPFR number.
+ *   px: value of x (string or MPFR number).
+ *   sy: value of y (string).
+ *   rnd: rounding mode.
+ *   z1: expected result (null pointer if unknown pure FP value).
+ *   inex1: expected ternary value (if z1 is not a null pointer).
+ *   z2: computed result.
+ *   inex2: computed ternary value.
+ *   flags1: expected flags (computed flags in __gmpfr_flags).
+ *   s: string about the context.
+ */
 static void
 cmpres (int spx, const void *px, const char *sy, mp_rnd_t rnd,
         mpfr_srcptr z1, int inex1, mpfr_srcptr z2, int inex2,
-        const char *s)
+        unsigned int flags1, const char *s)
 {
-  if (MPFR_IS_NAN (z1) && MPFR_IS_NAN (z2))
-    return;
-  if (mpfr_equal_p (z1, z2) && SAME_SIGN (inex1, inex2))
-    return;
+  unsigned int flags2 = __gmpfr_flags;
+
+  if (flags1 == flags2)
+    {
+      if (z1 == NULL)
+        {
+          if (MPFR_IS_PURE_FP (z2))
+            return;
+        }
+      else if (SAME_SIGN (inex1, inex2) &&
+               ((MPFR_IS_NAN (z1) && MPFR_IS_NAN (z2)) ||
+                mpfr_equal_p (z1, z2)))
+        return;
+    }
 
   printf ("Error with %s\nx = ", s);
   if (spx)
@@ -78,11 +100,18 @@ cmpres (int spx, const void *px, const char *sy, mp_rnd_t rnd,
     }
   printf ("y = %s, %s\n", sy, mpfr_print_rnd_mode (rnd));
   printf ("Expected ");
-  mpfr_out_str (stdout, 16, 0, z1, GMP_RNDN);
-  printf (", inex = %d\n", SIGN (inex1));
+  if (z1 == NULL)
+    {
+      printf ("pure FP value, flags = %u\n", flags1);
+    }
+  else
+    {
+      mpfr_out_str (stdout, 16, 0, z1, GMP_RNDN);
+      printf (", inex = %d, flags = %u\n", SIGN (inex1), flags1);
+    }
   printf ("Got      ");
   mpfr_out_str (stdout, 16, 0, z2, GMP_RNDN);
-  printf (", inex = %d\n", SIGN (inex2));
+  printf (", inex = %d, flags = %u\n", SIGN (inex2), flags2);
   if (all_cmpres_errors != 0)
     all_cmpres_errors = -1;
   else
@@ -103,7 +132,8 @@ is_odd (mpfr_srcptr x)
    rounding mode. */
 static void
 test_others (const void *sx, const char *sy, mp_rnd_t rnd,
-             mpfr_srcptr x, mpfr_srcptr y, mpfr_srcptr z1, int inex1)
+             mpfr_srcptr x, mpfr_srcptr y, mpfr_srcptr z1,
+             int inex1, unsigned int flags)
 {
   mpfr_t z2;
   int inex2;
@@ -116,7 +146,8 @@ test_others (const void *sx, const char *sy, mp_rnd_t rnd,
 
   __gmpfr_flags = MPFR_FLAGS_ALL;
   inex2 = mpfr_pow (z2, x, y, rnd);
-  cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, "mpfr_pow, flags set");
+  cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL,
+          "mpfr_pow, flags set");
 
   /* If y is an integer that fits in an unsigned long and is not -0,
      we can test mpfr_pow_ui. */
@@ -127,11 +158,11 @@ test_others (const void *sx, const char *sy, mp_rnd_t rnd,
 
       mpfr_clear_flags ();
       inex2 = mpfr_pow_ui (z2, x, yy, rnd);
-      cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2,
+      cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags,
               "mpfr_pow_ui, flags cleared");
       __gmpfr_flags = MPFR_FLAGS_ALL;
       inex2 = mpfr_pow_ui (z2, x, yy, rnd);
-      cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2,
+      cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL,
               "mpfr_pow_ui, flags set");
 
       /* If x is an integer that fits in an unsigned long and is not -0,
@@ -143,11 +174,11 @@ test_others (const void *sx, const char *sy, mp_rnd_t rnd,
 
           mpfr_clear_flags ();
           inex2 = mpfr_ui_pow_ui (z2, xx, yy, rnd);
-          cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2,
+          cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags,
                   "mpfr_ui_pow_ui, flags cleared");
           __gmpfr_flags = MPFR_FLAGS_ALL;
           inex2 = mpfr_ui_pow_ui (z2, xx, yy, rnd);
-          cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2,
+          cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL,
                   "mpfr_ui_pow_ui, flags set");
         }
     }
@@ -165,11 +196,11 @@ test_others (const void *sx, const char *sy, mp_rnd_t rnd,
 
           mpfr_clear_flags ();
           inex2 = mpfr_pow_si (z2, x, yy, rnd);
-          cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2,
+          cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags,
                   "mpfr_pow_si, flags cleared");
           __gmpfr_flags = MPFR_FLAGS_ALL;
           inex2 = mpfr_pow_si (z2, x, yy, rnd);
-          cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2,
+          cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL,
                   "mpfr_pow_si, flags set");
         }
 
@@ -178,11 +209,11 @@ test_others (const void *sx, const char *sy, mp_rnd_t rnd,
       mpfr_get_z (yyy, y, GMP_RNDN);
       mpfr_clear_flags ();
       inex2 = mpfr_pow_z (z2, x, yyy, rnd);
-      cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2,
+      cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags,
               "mpfr_pow_z, flags cleared");
       __gmpfr_flags = MPFR_FLAGS_ALL;
       inex2 = mpfr_pow_z (z2, x, yyy, rnd);
-      cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2,
+      cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL,
               "mpfr_pow_z, flags set");
       mpz_clear (yyy);
     }
@@ -196,11 +227,11 @@ test_others (const void *sx, const char *sy, mp_rnd_t rnd,
 
       mpfr_clear_flags ();
       inex2 = mpfr_ui_pow (z2, xx, y, rnd);
-      cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2,
+      cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags,
               "mpfr_ui_pow, flags cleared");
       __gmpfr_flags = MPFR_FLAGS_ALL;
       inex2 = mpfr_ui_pow (z2, xx, y, rnd);
-      cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2,
+      cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL,
               "mpfr_ui_pow, flags set");
 
       /* If x = 2, we can test mpfr_exp2. */
@@ -208,11 +239,11 @@ test_others (const void *sx, const char *sy, mp_rnd_t rnd,
         {
           mpfr_clear_flags ();
           inex2 = mpfr_exp2 (z2, y, rnd);
-          cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2,
+          cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags,
                   "mpfr_exp2, flags cleared");
           __gmpfr_flags = MPFR_FLAGS_ALL;
           inex2 = mpfr_exp2 (z2, y, rnd);
-          cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2,
+          cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL,
                   "mpfr_exp2, flags set");
         }
 
@@ -221,11 +252,11 @@ test_others (const void *sx, const char *sy, mp_rnd_t rnd,
         {
           mpfr_clear_flags ();
           inex2 = mpfr_exp10 (z2, y, rnd);
-          cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2,
+          cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags,
                   "mpfr_exp10, flags cleared");
           __gmpfr_flags = MPFR_FLAGS_ALL;
           inex2 = mpfr_exp10 (z2, y, rnd);
-          cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2,
+          cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL,
                   "mpfr_exp10, flags set");
         }
     }
@@ -248,6 +279,7 @@ tst (void)
       RND_LOOP (rnd)
         {
           int exact, inex;
+          unsigned int flags;
 
           if (mpfr_set_str (x, val[i], 10, GMP_RNDN) ||
               mpfr_set_str (y, val[j], 10, GMP_RNDN))
@@ -257,6 +289,7 @@ tst (void)
             }
           mpfr_clear_flags ();
           inex = mpfr_pow (z, x, y, (mp_rnd_t) rnd);
+          flags = __gmpfr_flags;
           if (mpfr_underflow_p ())
             err ("got underflow", i, j, rnd, z, inex);
           if (mpfr_overflow_p ())
@@ -328,13 +361,99 @@ tst (void)
               if ((MPFR_IS_NEG (x) && is_odd (y)) ^ MPFR_IS_NEG (z))
                 err ("wrong sign", i, j, rnd, z, inex);
             }
-          test_others (val[i], val[j], (mp_rnd_t) rnd, x, y, z, inex);
+          test_others (val[i], val[j], (mp_rnd_t) rnd, x, y, z, inex, flags);
         }
   mpfr_clears (x, y, z, tmp, (mpfr_ptr) 0);
 }
 
 static void
-underflow_up (int extended_emin)
+underflow_up1 (int extended_emin)
+{
+  mpfr_t delta, x, y, z, z0;
+  mp_exp_t n;
+  int inex;
+  int rnd;
+  int i;
+
+  n = mpfr_get_emin ();
+  if (n < LONG_MIN)
+    return;
+
+  mpfr_init2 (delta, 2);
+  inex = mpfr_set_ui_2exp (delta, 1, -2, GMP_RNDN);
+  MPFR_ASSERTN (inex == 0);
+
+  mpfr_init2 (x, 8);
+  inex = mpfr_set_ui (x, 2, GMP_RNDN);
+  MPFR_ASSERTN (inex == 0);
+
+  mpfr_init2 (y, sizeof (long) * CHAR_BIT + 4);
+  inex = mpfr_set_si (y, n, GMP_RNDN);
+  MPFR_ASSERTN (inex == 0);
+
+  mpfr_init2 (z0, 2);
+  mpfr_set_ui (z0, 0, GMP_RNDN);
+
+  mpfr_init2 (z, 32);
+
+  for (i = 0; i <= 12; i++)
+    {
+      unsigned int flags = 0;
+      char sy[16];
+
+      /* Test 2^(emin - i/4).
+       * --> Underflow iff i > 4.
+       * --> Zero in GMP_RNDN iff i >= 8.
+       */
+
+      if (i != 0 && i != 4)
+        flags |= MPFR_FLAGS_INEXACT;
+      if (i > 4)
+        flags |= MPFR_FLAGS_UNDERFLOW;
+
+      sprintf (sy, "emin - %d/4", i);
+
+      RND_LOOP (rnd)
+        {
+          int zero;
+
+          zero = (i > 4 && (rnd == GMP_RNDZ || rnd == GMP_RNDD)) ||
+            (i >= 8 && rnd == GMP_RNDN);
+
+          mpfr_clear_flags ();
+          inex = mpfr_pow (z, x, y, (mp_rnd_t) rnd);
+          cmpres (1, "2", sy, (mp_rnd_t) rnd,
+                  zero ? z0 : (mpfr_ptr) NULL, -1, z, inex, flags,
+                  extended_emin ? "underflow_up1 and extended emin" :
+                  "underflow_up1");
+          test_others ("2", sy, (mp_rnd_t) rnd, x, y, z, inex, flags);
+        }
+
+      inex = mpfr_sub (y, y, delta, GMP_RNDN);
+      MPFR_ASSERTN (inex == 0);
+    }
+
+  mpfr_clears (delta, x, y, z, z0, (mpfr_ptr) 0);
+}
+
+/* With pow.c r5497, the following test fails on a 64-bit Linux machine
+ * due to a double-rounding problem when rescaling the result:
+ *   Error with underflow_up2 and extended emin
+ *   x = 7.fffffffffffffff0@-1,
+ *   y = 4611686018427387904, GMP_RNDN
+ *   Expected 1.0000000000000000@-1152921504606846976, inex = 1, flags = 9
+ *   Got      0, inex = -1, flags = 9
+ * With pow_ui.c r5423, the following test fails on a 64-bit Linux machine
+ * as underflows and overflows are not handled correctly (the approximation
+ * error is ignored):
+ *   Error with mpfr_pow_ui, flags cleared
+ *   x = 7.fffffffffffffff0@-1,
+ *   y = 4611686018427387904, GMP_RNDN
+ *   Expected 1.0000000000000000@-1152921504606846976, inex = 1, flags = 9
+ *   Got      0, inex = -1, flags = 9
+ */
+static void
+underflow_up2 (int extended_emin)
 {
   mpfr_t x, y, z, z0, eps;
   mp_exp_t n;
@@ -369,32 +488,30 @@ underflow_up (int extended_emin)
       int expected_inex;
       char sy[256];
 
-      mpfr_clear_flags ();
-      inex = mpfr_pow (z, x, y, (mp_rnd_t) rnd);
-      if (__gmpfr_flags != ufinex)
-        {
-          printf ("Error in underflow_up for %s",
-                  mpfr_print_rnd_mode ((mp_rnd_t) rnd));
-          if (extended_emin)
-            printf (" and extended emin");
-          printf ("\n");
-          printf ("got %u instead of %u\n", __gmpfr_flags, ufinex);
-          exit (1);
-        }
       mpfr_set_ui (z0, 0, GMP_RNDN);
       expected_inex = rnd == GMP_RNDN || rnd == GMP_RNDU ?
         (mpfr_nextabove (z0), 1) : -1;
       sprintf (sy, "%lu", (unsigned long) n);
-      cmpres (0, x, sy, (mp_rnd_t) rnd, z0, expected_inex, z, inex,
-              extended_emin ? "underflow_up and extended emin" :
-              "underflow_up");
-      test_others (NULL, sy, (mp_rnd_t) rnd, x, y, z, inex);
+
+      mpfr_clear_flags ();
+      inex = mpfr_pow (z, x, y, (mp_rnd_t) rnd);
+      cmpres (0, x, sy, (mp_rnd_t) rnd, z0, expected_inex, z, inex, ufinex,
+              extended_emin ? "underflow_up2 and extended emin" :
+              "underflow_up2");
+      test_others (NULL, sy, (mp_rnd_t) rnd, x, y, z, inex, ufinex);
     }
 
   mpfr_clears (x, y, z, z0, eps, (mpfr_ptr) 0);
 }
 
 static void
+underflow_up (int extended_emin)
+{
+  underflow_up1 (extended_emin);
+  underflow_up2 (extended_emin);
+}
+
+static void
 underflow (void)
 {
   mp_exp_t emin;