Fixed various bugs in functions mpfr_pow, mpfr_pow_ui, mpfr_pow_si and

mpfr_pow_z (overflow and underflow handling, and mpfr_pow(z, x, y, rnd)
with integer |y| >= 2^256) and added testcases for these bugs (and more)
by merging
  -r5313:5318 -r5415:5417 -r5426:5434 -c5445 -c5455 -c5505
  -r5529:5533 -r5534:5544 -r5546:5557
from the trunk. Details:
* pow.c:
  - Moved the general case from mpfr_pow() to a new internal function
    mpfr_pow_general().
  - Fixed overflow/underflow handling.
  - Fixed a bug in mpfr_pow(z, x, y, rnd) that occurred with integer
    |y| >= 2^256 (unless overflow or underflow had been detected).
  - Added log messages.
* mpfr-impl.h: added mpfr_pow_general prototype.
* pow_z.c:
  - Fixed overflow/underflow handling.
    + The underflow cases in rounding to nearest are 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!).
    + If z < 0, internally evaluate X^Z rounded towards 1 (or -1) to
      avoid a spurious overflow/underflow, where X = 1/x and Z = -z.
  - Added log messages.
* pow_si.c:
  - Fixed overflow/underflow handling for n < 0.
  - 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.
  - Added log messages.
* tests/tpow.c:
  - Added specific testcase for the bug with integer |y| >= 2^256.
  - Added specific testcase for a bug from r4940 (the rounding modes
    to compute an upper bound on exp(y*ln|x|) are incorrect); this
    bug could affect only underflow cases and possibly cases near
    overflow.
* tests/tpow_all.c:
  - Added several tests of all power functions in MPFR or equivalent
    (including mpfr_ui_div, mpfr_exp2 and mpfr_exp10), in particular
    in underflow and overflow cases.
* tests/tpow_z.c:
  - Added specific testcase for the mpfr_pow_z bug with z < 0.

git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/branches/2.3@5558 280ebfd0-de03-0410-8827-d642c229c3f4

8 files changed, 1270 insertions, 431 deletions

diff --git a/mpfr-impl.h b/mpfr-impl.h
index 943df1b8c..63bc2542a 100644
--- a/mpfr-impl.h
+++ b/mpfr-impl.h
@@ -567,8 +567,18 @@ typedef unsigned long int  mpfr_uexp_t;
 #ifndef mp_exp_unsigned_t
 # define mp_exp_unsigned_t mpfr_uexp_t
 #endif
-#define mpfr_get_exp_t(x,r) mpfr_get_si((x),(r))
-#define mpfr_set_exp_t(x,e,r) mpfr_set_si((x),(e),(r))
+
+#if MPFR_EXP_MIN >= LONG_MIN && MPFR_EXP_MAX <= LONG_MAX
+typedef long int mpfr_eexp_t;
+# define mpfr_get_exp_t(x,r) mpfr_get_si((x),(r))
+# define mpfr_set_exp_t(x,e,r) mpfr_set_si((x),(e),(r))
+#elif defined (_MPFR_H_HAVE_INTMAX_T)
+typedef intmax_t mpfr_eexp_t;
+# define mpfr_get_exp_t(x,r) mpfr_get_sj((x),(r))
+# define mpfr_set_exp_t(x,e,r) mpfr_set_sj((x),(e),(r))
+#else
+# error "Cannot define mpfr_get_exp_t and mpfr_set_exp_t"
+#endif
 
 /* Invalid exponent value (to track bugs...) */
 #define MPFR_EXP_INVALID \
@@ -1547,6 +1557,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 53c33bd45..b910b4352 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, 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_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, 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;
+        }
+      /* 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,32 +532,51 @@ 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);
         }
     }
 
-  /* detect underflows: for x > 0, y < 0, |x^y| = |(1/x)^(-y)|
-                        <= 2^((1-EXP(x))*(-y)) */
-  if (MPFR_IS_NEG(y) && MPFR_EXP(x) > 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;
-      int negative, underflow;
+      mpfr_eexp_t ebound;
+      int inex2;
 
-      /* We must restore the flags if no underflow. */
+      /* We must restore the flags. */
       MPFR_SAVE_EXPO_MARK (expo);
-      mpfr_init2 (tmp, 53);
-      mpfr_neg (tmp, y, GMP_RNDZ);
-      mpfr_mul_si (tmp, tmp, 1 - MPFR_EXP(x), GMP_RNDZ);
-      underflow = mpfr_cmp_si (tmp, __gmpfr_emin - 2) <= 0;
+      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 (underflow)
+      if (MPFR_UNLIKELY (ebound <=
+                         __gmpfr_emin - (rnd_mode == GMP_RNDN ? 2 : 1)))
         {
           /* warning: mpfr_underflow rounds away from 0 for GMP_RNDN */
-          negative = MPFR_SIGN(x) < 0 && is_odd (y);
-          return mpfr_underflow (z, (rnd_mode == GMP_RNDN) ? GMP_RNDZ :
-                                 rnd_mode, negative ? -1 : 1);
+          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);
         }
     }
 
@@ -365,11 +584,15 @@ mpfr_pow (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd_mode)
      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). */
+     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);
@@ -388,6 +611,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);
@@ -395,9 +619,9 @@ mpfr_pow (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd_mode)
         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 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). */
+           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);
@@ -441,155 +665,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));
-
-    /* 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 */
-        if (check_exact_case == 0)
-          {
-            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_si.c b/pow_si.c
index f4e4c1e40..11a26ae8f 100644
--- a/pow_si.c
+++ b/pow_si.c
@@ -31,6 +31,9 @@ MA 02110-1301, USA. */
 int
 mpfr_pow_si (mpfr_ptr y, mpfr_srcptr x, long int n, mp_rnd_t rnd)
 {
+  MPFR_LOG_FUNC (("x[%#R]=%R n=%ld rnd=%d", x, x, n, rnd),
+                 ("y[%#R]=%R", y, y));
+
   if (n >= 0)
     return mpfr_pow_ui (y, x, n, rnd);
   else
@@ -131,52 +134,98 @@ mpfr_pow_si (mpfr_ptr y, mpfr_srcptr x, long int n, mp_rnd_t rnd)
         /* Declaration of the intermediary variable */
         mpfr_t t;
         /* Declaration of the size variable */
-        mp_prec_t Ny = MPFR_PREC (y);               /* target precision */
+        mp_prec_t Ny;                              /* target precision */
         mp_prec_t Nt;                              /* working precision */
-        mp_exp_t  err;                             /* error */
+        mp_rnd_t rnd1;
+        int size_n;
         int inexact;
         unsigned long abs_n;
         MPFR_SAVE_EXPO_DECL (expo);
         MPFR_ZIV_DECL (loop);
 
         abs_n = - (unsigned long) n;
+        count_leading_zeros (size_n, (mp_limb_t) abs_n);
+        size_n = BITS_PER_MP_LIMB - size_n;
 
-        /* compute the precision of intermediary variable */
-        /* the optimal number of bits : see algorithms.tex */
-        Nt = Ny + 3 + MPFR_INT_CEIL_LOG2 (Ny);
+        /* initial working precision */
+        Ny = MPFR_PREC (y);
+        Nt = Ny + size_n + 3 + MPFR_INT_CEIL_LOG2 (Ny);
 
         MPFR_SAVE_EXPO_MARK (expo);
 
         /* initialise of intermediary   variable */
         mpfr_init2 (t, Nt);
 
+        /* We will compute rnd(rnd1(1/x) ^ |n|), where rnd1 is the rounding
+           toward sign(x), to avoid spurious overflow or underflow, as in
+           mpfr_pow_z. */
+        rnd1 = MPFR_EXP (x) < 1 ? GMP_RNDZ :
+          (MPFR_SIGN (x) > 0 ? GMP_RNDU : GMP_RNDD);
+
         MPFR_ZIV_INIT (loop, Nt);
         for (;;)
           {
-            /* compute 1/(x^n), with n > 0 */
-            mpfr_pow_ui (t, x, abs_n, GMP_RNDN);
-            mpfr_ui_div (t, 1, t, GMP_RNDN);
-            /* FIXME: old code improved, but I think this is still incorrect. */
-            if (MPFR_UNLIKELY (MPFR_IS_ZERO (t)))
+            MPFR_BLOCK_DECL (flags);
+
+            /* compute (1/x)^|n| */
+            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_ui (t, t, abs_n, rnd));
+            /* t = (1/x)^|n|*(1+theta')^(|n|+1) where |theta'| <= 2^(-Nt).
+               If (|n|+1)*2^(-Nt) <= 1/2, which is satisfied as soon as
+               Nt >= bits(n)+2, then we can use Lemma \ref{lemma_graillat}
+               from algorithms.tex, which yields x^n*(1+theta) with
+               |theta| <= 2(|n|+1)*2^(-Nt), thus the error is bounded by
+               2(|n|+1) ulps <= 2^(bits(n)+2) ulps. */
+            if (MPFR_UNLIKELY (MPFR_OVERFLOW (flags)))
               {
+              overflow:
                 MPFR_ZIV_FREE (loop);
                 mpfr_clear (t);
                 MPFR_SAVE_EXPO_FREE (expo);
-                return mpfr_underflow (y, rnd == GMP_RNDN ? GMP_RNDZ : rnd,
-                                       abs_n & 1 ? MPFR_SIGN (x) :
-                                       MPFR_SIGN_POS);
+                MPFR_LOG_MSG (("overflow\n", 0));
+                return mpfr_overflow (y, rnd, abs_n & 1 ?
+                                      MPFR_SIGN (x) : MPFR_SIGN_POS);
               }
-            if (MPFR_UNLIKELY (MPFR_IS_INF (t)))
+            if (MPFR_UNLIKELY (MPFR_UNDERFLOW (flags)))
               {
                 MPFR_ZIV_FREE (loop);
                 mpfr_clear (t);
-                MPFR_SAVE_EXPO_FREE (expo);
-                return mpfr_overflow (y, rnd, abs_n & 1 ? MPFR_SIGN (x) :
-                                      MPFR_SIGN_POS);
+                MPFR_LOG_MSG (("underflow\n", 0));
+                if (rnd == GMP_RNDN)
+                  {
+                    mpfr_t y2, nn;
+
+                    /* 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 (nn, sizeof (long) * CHAR_BIT);
+                    inexact = mpfr_set_si (nn, n, GMP_RNDN);
+                    MPFR_ASSERTN (inexact == 0);
+                    inexact = mpfr_pow_general (y2, x, nn, rnd, 1,
+                                                (mpfr_save_expo_t *) NULL);
+                    mpfr_clear (nn);
+                    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, abs_n & 1 ?
+                                           MPFR_SIGN (x) : MPFR_SIGN_POS);
+                  }
               }
             /* error estimate -- see pow function in algorithms.ps */
-            err = Nt - 3;
-            if (MPFR_LIKELY (MPFR_CAN_ROUND (t, err, Ny, rnd)))
+            if (MPFR_LIKELY (MPFR_CAN_ROUND (t, Nt - size_n - 2, Ny, rnd)))
               break;
 
             /* actualisation of the precision */
@@ -187,6 +236,8 @@ mpfr_pow_si (mpfr_ptr y, mpfr_srcptr x, long int n, mp_rnd_t rnd)
 
         inexact = mpfr_set (y, t, rnd);
         mpfr_clear (t);
+
+      end:
         MPFR_SAVE_EXPO_FREE (expo);
         return mpfr_check_range (y, inexact, rnd);
       }
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 2c8b545f2..5a035b105 100644
--- a/pow_z.c
+++ b/pow_z.c
@@ -23,59 +23,77 @@ 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)
+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;
+  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);
 
-  rnd1 = MPFR_IS_POS (x) ? GMP_RNDU : GMP_RNDD; /* away */
   absz[0] = z[0];
   SIZ (absz) = ABS(SIZ(absz)); /* Hack to get abs(z) */
   MPFR_MPZ_SIZEINBASE2 (size_z, z);
 
-  prec = MPFR_PREC (y) + 3 + size_z + MPFR_INT_CEIL_LOG2 (MPFR_PREC (y));
+  /* 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 pow_ui.c for the error analusis, which is identical */
+      /* 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;
 
-      /* First step: compute square from y */
       MPFR_BLOCK (flags,
-                  inexact = mpfr_mul (res, x, x, GMP_RNDU);
+                  inexmul = mpfr_mul (res, x, x, rnd2);
                   MPFR_ASSERTD (i >= 2);
                   if (mpz_tstbit (absz, i - 2))
-                    inexact |= mpfr_mul (res, res, x, rnd1);
+                    inexmul |= mpfr_mul (res, res, x, rnd1);
                   for (i -= 3; i >= 0 && !MPFR_BLOCK_EXCEP; i--)
                     {
-                      inexact |= mpfr_mul (res, res, res, GMP_RNDU);
+                      inexmul |= mpfr_mul (res, res, res, rnd2);
                       if (mpz_tstbit (absz, i))
-                        inexact |= mpfr_mul (res, res, x, rnd1);
+                        inexmul |= mpfr_mul (res, res, x, rnd1);
                     });
-      /*    printf ("pow_z ");
-            mpfr_dump_mant (MPFR_MANT (res), prec, MPFR_PREC (x), err); */
-      if (MPFR_LIKELY (inexact == 0
+      if (MPFR_LIKELY (inexmul == 0 || cr == 0
                        || MPFR_OVERFLOW (flags) || MPFR_UNDERFLOW (flags)
                        || MPFR_CAN_ROUND (res, err, MPFR_PREC (y), rnd)))
         break;
-      /* Actualisation of the precision */
+      /* Can't decide correct rounding, increase the precision */
       MPFR_ZIV_NEXT (loop, prec);
       mpfr_set_prec (res, prec);
     }
@@ -83,13 +101,46 @@ mpfr_pow_pos_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mp_rnd_t rnd)
 
   /* 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);
+    {
+      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);
 
@@ -97,10 +148,16 @@ mpfr_pow_pos_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mp_rnd_t rnd)
   return inexact;
 }
 
-/* The computation of y=pow(x,z) is done by
- *    y=pow_ui(x,z) if z>0
- * else
- *    y=1/pow_ui(x,z) if z<0
+/* 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
@@ -110,6 +167,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);
@@ -153,35 +213,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);
@@ -191,49 +252,102 @@ mpfr_pow_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mp_rnd_t rnd)
   MPFR_SAVE_EXPO_MARK (expo);
 
   if (mpz_sgn (z) > 0)
-    inexact = mpfr_pow_pos_z (y, x, z, rnd);
+    {
+      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);
 
-      /* compute the precision of intermediary variable */
-      Nt = MAX (MPFR_PREC (x), MPFR_PREC (y));
+      MPFR_MPZ_SIZEINBASE2 (size_z, z);
 
-      /* the optimal number of bits : see algorithms.ps */
-      Nt = Nt + 3 + MPFR_INT_CEIL_LOG2 (Nt);
+      /* 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 */
+      /* 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 (;;)
         {
-          /* compute 1/(x^n) n>0 */
-          mpfr_pow_pos_z (t, x, z, GMP_RNDN);
-          mpfr_ui_div (t, 1, t, GMP_RNDN);
-          /* FIXME: old code improved, but I think this is still incorrect. */
-          if (MPFR_UNLIKELY (MPFR_IS_ZERO (t)))
+          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);
-              return mpfr_underflow (y, rnd == GMP_RNDN ? GMP_RNDZ : rnd,
-                                     mpz_odd_p (z) ? MPFR_SIGN (x) :
-                                     MPFR_SIGN_POS);
+              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_IS_INF (t)))
+          if (MPFR_UNLIKELY (MPFR_UNDERFLOW (flags)))
             {
               MPFR_ZIV_FREE (loop);
               mpfr_clear (t);
-              MPFR_SAVE_EXPO_FREE (expo);
-              return mpfr_overflow (y, rnd,
-                                    mpz_odd_p (z) ? MPFR_SIGN (x) :
-                                    MPFR_SIGN_POS);
+              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-3, MPFR_PREC (y), rnd)))
+          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);
@@ -245,7 +359,7 @@ mpfr_pow_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mp_rnd_t rnd)
       mpfr_clear (t);
     }
 
-  MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);
+ end:
   MPFR_SAVE_EXPO_FREE (expo);
   return mpfr_check_range (y, inexact, rnd);
 }
diff --git a/tests/tpow.c b/tests/tpow.c
index 470f21cf6..85dd715f6 100644
--- a/tests/tpow.c
+++ b/tests/tpow.c
@@ -1243,6 +1243,97 @@ bug20071218 (void)
   mpfr_clears (x, y, z, t, (mpfr_ptr) 0);
 }
 
+/* With revision 5429, this gives:
+ *   pow.c:43:  assertion failed: !mpfr_integer_p (y)
+ * This is fixed in revision 5432.
+ */
+static void
+bug20080721 (void)
+{
+  mpfr_t x, y, z, t[2];
+  int inex;
+  int rnd;
+  int err = 0;
+
+  /* Note: input values have been chosen in a way to select the
+   * general case. If mpfr_pow is modified, in particular line
+   *     if (y_is_integer && (MPFR_GET_EXP (y) <= 256))
+   * make sure that this test still does what we want.
+   */
+  mpfr_inits2 (4913, x, y, (mpfr_ptr) 0);
+  mpfr_inits2 (8, z, t[0], t[1], (mpfr_ptr) 0);
+  mpfr_set_si (x, -1, GMP_RNDN);
+  mpfr_nextbelow (x);
+  mpfr_set_ui_2exp (y, 1, mpfr_get_prec (y) - 1, GMP_RNDN);
+  inex = mpfr_add_ui (y, y, 1, GMP_RNDN);
+  MPFR_ASSERTN (inex == 0);
+  mpfr_set_str_binary (t[0], "-0.10101101e2");
+  mpfr_set_str_binary (t[1], "-0.10101110e2");
+  RND_LOOP (rnd)
+    {
+      int i, inex0;
+
+      i = (rnd == GMP_RNDN || rnd == GMP_RNDD);
+      inex0 = i ? -1 : 1;
+      mpfr_clear_flags ();
+      inex = mpfr_pow (z, x, y, (mp_rnd_t) rnd);
+      if (__gmpfr_flags != MPFR_FLAGS_INEXACT || ! SAME_SIGN (inex, inex0)
+          || MPFR_IS_NAN (z) || mpfr_cmp (z, t[i]) != 0)
+        {
+          unsigned int flags = __gmpfr_flags;
+
+          printf ("Error in bug20080721 with %s\n",
+                  mpfr_print_rnd_mode ((mp_rnd_t) rnd));
+          printf ("expected ");
+          mpfr_out_str (stdout, 2, 0, t[i], GMP_RNDN);
+          printf (", inex = %d, flags = %u\n", inex0, MPFR_FLAGS_INEXACT);
+          printf ("got      ");
+          mpfr_out_str (stdout, 2, 0, z, GMP_RNDN);
+          printf (", inex = %d, flags = %u\n", inex, flags);
+          err = 1;
+        }
+    }
+  mpfr_clears (x, y, z, t[0], t[1], (mpfr_ptr) 0);
+  if (err)
+    exit (1);
+}
+
+/* The following test fails in r5552 (32-bit and 64-bit). This is due to:
+ *   mpfr_log (t, absx, GMP_RNDU);
+ *   mpfr_mul (t, y, t, GMP_RNDU);
+ * in pow.c, that is supposed to compute an upper bound on exp(y*ln|x|),
+ * but this is incorrect if y is negative.
+ */
+static void
+bug20080820 (void)
+{
+  mp_exp_t emin;
+  mpfr_t x, y, z1, z2;
+
+  emin = mpfr_get_emin ();
+  mpfr_set_emin (MPFR_EMIN_MIN);
+  mpfr_init2 (x, 80);
+  mpfr_init2 (y, sizeof (mp_exp_t) * CHAR_BIT + 32);
+  mpfr_init2 (z1, 2);
+  mpfr_init2 (z2, 80);
+  mpfr_set_ui (x, 2, GMP_RNDN);
+  mpfr_nextbelow (x);
+  mpfr_set_exp_t (y, mpfr_get_emin () - 2, GMP_RNDN);
+  mpfr_nextabove (y);
+  mpfr_pow (z1, x, y, GMP_RNDN);
+  mpfr_pow (z2, x, y, GMP_RNDN);
+  /* As x > 0, the rounded value of x^y to nearest in precision p is equal
+     to 0 iff x^y <= 2^(emin - 2). In particular, this does not depend on
+     the precision p. Hence the following test. */
+  if (MPFR_IS_ZERO (z1) && MPFR_NOTZERO (z2))
+    {
+      printf ("Error in bug20080820\n");
+      exit (1);
+    }
+  mpfr_clears (x, y, z1, z2, (mpfr_ptr) 0);
+  set_emin (emin);
+}
+
 int
 main (int argc, char **argv)
 {
@@ -1268,6 +1359,8 @@ main (int argc, char **argv)
   bug20071104 ();
   bug20071128 ();
   bug20071218 ();
+  bug20080721 ();
+  bug20080820 ();
 
   test_generic (2, 100, 100);
   test_generic_ui (2, 100, 100);
diff --git a/tests/tpow_all.c b/tests/tpow_all.c
index ddd59de28..5e25cb643 100644
--- a/tests/tpow_all.c
+++ b/tests/tpow_all.c
@@ -26,18 +26,35 @@ MA 02110-1301, USA. */
    expected result, in particular because such special cases are handled
    in different ways in each function. */
 
+/* Execute with at least an argument to report all the errors found by
+   comparisons. */
+
+#include <stdio.h>
 #include <stdlib.h>
+#include <limits.h>
 
 #include "mpfr-test.h"
 
+/* Behavior of cmpres (called by test_others):
+ *   0: stop as soon as an error is found.
+ *   1: report all errors found by test_others.
+ *  -1: the 1 is changed to this value as soon as an error has been found.
+ */
+static int all_cmpres_errors;
+
+/* Non-zero when extended exponent range */
+static int ext = 0;
+
 static char *val[] =
-  { "@NaN@", "-@Inf@", "-4", "-3", "-2", "-1.5", "-1", "-0.5", "-0",
-    "0", "0.5", "1", "1.5", "2", "3", "4", "@Inf@" };
+  { "min", "min+", "max", "@NaN@", "-@Inf@", "-4", "-3", "-2", "-1.5",
+    "-1", "-0.5", "-0", "0", "0.5", "1", "1.5", "2", "3", "4", "@Inf@" };
 
 static void
 err (const char *s, int i, int j, int rnd, mpfr_srcptr z, int inex)
 {
   puts (s);
+  if (ext)
+    puts ("extended exponent range");
   printf ("x = %s, y = %s, %s\n", val[i], val[j],
           mpfr_print_rnd_mode ((mp_rnd_t) rnd));
   printf ("z = ");
@@ -46,32 +63,252 @@ 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).
+ *   s1, s2: strings about the context.
+ */
 static void
-cmpres (int i, int j, int rnd, mpfr_srcptr z1, int inex1,
-        mpfr_srcptr z2, int inex2, const char *s)
+cmpres (int spx, const void *px, const char *sy, mp_rnd_t rnd,
+        mpfr_srcptr z1, int inex1, mpfr_srcptr z2, int inex2,
+        unsigned int flags1, const char *s1, const char *s2)
 {
-  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;
 
-  printf ("Error with %s\n", s);
-  printf ("x = %s, y = %s, %s\n", val[i], val[j],
-          mpfr_print_rnd_mode ((mp_rnd_t) rnd));
+  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 in %s\nwith %s%s\nx = ", s1, s2,
+          ext ? ", extended exponent range" : "");
+  if (spx)
+    printf ("%s, ", (char *) px);
+  else
+    {
+      mpfr_out_str (stdout, 16, 0, (mpfr_ptr) px, GMP_RNDN);
+      puts (",");
+    }
+  printf ("y = %s, %s\n", sy, mpfr_print_rnd_mode (rnd));
   printf ("Expected ");
-  mpfr_out_str (stdout, 10, 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, 10, 0, z2, GMP_RNDN);
-  printf (", inex = %d\n", SIGN (inex2));
-  exit (1);
+  mpfr_out_str (stdout, 16, 0, z2, GMP_RNDN);
+  printf (", inex = %d, flags = %u\n", SIGN (inex2), flags2);
+  if (all_cmpres_errors != 0)
+    all_cmpres_errors = -1;
+  else
+    exit (1);
 }
 
 static int
 is_odd (mpfr_srcptr x)
 {
-  /* does not work for large integers */
-  return mpfr_integer_p (x) && (mpfr_get_si (x, GMP_RNDN) & 1);
+  /* works only with the values from val[] */
+  return mpfr_integer_p (x) && mpfr_fits_slong_p (x, GMP_RNDN) &&
+    (mpfr_get_si (x, GMP_RNDN) & 1);
+}
+
+/* Compare the result (z1,inex1) of mpfr_pow with all flags cleared
+   with those of mpfr_pow with all flags set and of the other power
+   functions. Arguments x and y are the input values; sx and sy are
+   their string representations (sx may be null); rnd contains the
+   rounding mode; s is a string containing the function that called
+   test_others. */
+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, unsigned int flags, const char *s)
+{
+  mpfr_t z2;
+  int inex2;
+  int spx = sx != NULL;
+
+  if (!spx)
+    sx = x;
+
+  mpfr_init2 (z2, mpfr_get_prec (z1));
+
+  __gmpfr_flags = MPFR_FLAGS_ALL;
+  inex2 = mpfr_pow (z2, x, y, rnd);
+  cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL,
+          s, "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. */
+  if (MPFR_IS_POS (y) && mpfr_integer_p (y) &&
+      mpfr_fits_ulong_p (y, GMP_RNDN))
+    {
+      unsigned long yy = mpfr_get_ui (y, GMP_RNDN);
+
+      mpfr_clear_flags ();
+      inex2 = mpfr_pow_ui (z2, x, yy, rnd);
+      cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags,
+              s, "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, MPFR_FLAGS_ALL,
+              s, "mpfr_pow_ui, flags set");
+
+      /* If x is an integer that fits in an unsigned long and is not -0,
+         we can also test mpfr_ui_pow_ui. */
+      if (MPFR_IS_POS (x) && mpfr_integer_p (x) &&
+          mpfr_fits_ulong_p (x, GMP_RNDN))
+        {
+          unsigned long xx = mpfr_get_ui (x, GMP_RNDN);
+
+          mpfr_clear_flags ();
+          inex2 = mpfr_ui_pow_ui (z2, xx, yy, rnd);
+          cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags,
+                  s, "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, MPFR_FLAGS_ALL,
+                  s, "mpfr_ui_pow_ui, flags set");
+        }
+    }
+
+  /* If y is an integer but not -0 and not huge, we can test mpfr_pow_z,
+     and possibly mpfr_pow_si (and possibly mpfr_ui_div). */
+  if (MPFR_IS_ZERO (y) ? MPFR_IS_POS (y) :
+      (mpfr_integer_p (y) && MPFR_GET_EXP (y) < 256))
+    {
+      mpz_t yyy;
+
+      /* If y fits in a long, we can test mpfr_pow_si. */
+      if (mpfr_fits_slong_p (y, GMP_RNDN))
+        {
+          long yy = mpfr_get_si (y, GMP_RNDN);
+
+          mpfr_clear_flags ();
+          inex2 = mpfr_pow_si (z2, x, yy, rnd);
+          cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags,
+                  s, "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, MPFR_FLAGS_ALL,
+                  s, "mpfr_pow_si, flags set");
+
+          /* If y is -1, we can test mpfr_ui_div. */
+          if (yy == -1)
+            {
+              mpfr_clear_flags ();
+              inex2 = mpfr_ui_div (z2, 1, x, rnd);
+              cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags,
+                      s, "mpfr_ui_div, flags cleared");
+              __gmpfr_flags = MPFR_FLAGS_ALL;
+              inex2 = mpfr_ui_div (z2, 1, x, rnd);
+              cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL,
+                      s, "mpfr_ui_div, flags set");
+            }
+        }
+
+      /* Test mpfr_pow_z. */
+      mpz_init (yyy);
+      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, flags,
+              s, "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, MPFR_FLAGS_ALL,
+              s, "mpfr_pow_z, flags set");
+      mpz_clear (yyy);
+    }
+
+  /* If x is an integer that fits in an unsigned long and is not -0,
+     we can test mpfr_ui_pow. */
+  if (MPFR_IS_POS (x) && mpfr_integer_p (x) &&
+      mpfr_fits_ulong_p (x, GMP_RNDN))
+    {
+      unsigned long xx = mpfr_get_ui (x, GMP_RNDN);
+
+      mpfr_clear_flags ();
+      inex2 = mpfr_ui_pow (z2, xx, y, rnd);
+      cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags,
+              s, "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, MPFR_FLAGS_ALL,
+              s, "mpfr_ui_pow, flags set");
+
+      /* If x = 2, we can test mpfr_exp2. */
+      if (xx == 2)
+        {
+          mpfr_clear_flags ();
+          inex2 = mpfr_exp2 (z2, y, rnd);
+          cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags,
+                  s, "mpfr_exp2, flags cleared");
+          __gmpfr_flags = MPFR_FLAGS_ALL;
+          inex2 = mpfr_exp2 (z2, y, rnd);
+          cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL,
+                  s, "mpfr_exp2, flags set");
+        }
+
+      /* If x = 10, we can test mpfr_exp10. */
+      if (xx == 10)
+        {
+          mpfr_clear_flags ();
+          inex2 = mpfr_exp10 (z2, y, rnd);
+          cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, flags,
+                  s, "mpfr_exp10, flags cleared");
+          __gmpfr_flags = MPFR_FLAGS_ALL;
+          inex2 = mpfr_exp10 (z2, y, rnd);
+          cmpres (spx, sx, sy, rnd, z1, inex1, z2, inex2, MPFR_FLAGS_ALL,
+                  s, "mpfr_exp10, flags set");
+        }
+    }
+
+  mpfr_clear (z2);
+}
+
+static int
+my_setstr (mpfr_ptr t, const char *s)
+{
+  if (strcmp (s, "min") == 0)
+    {
+      mpfr_setmin (t, mpfr_get_emin ());
+      MPFR_SET_POS (t);
+      return 0;
+    }
+  if (strcmp (s, "min+") == 0)
+    {
+      mpfr_setmin (t, mpfr_get_emin ());
+      MPFR_SET_POS (t);
+      mpfr_nextabove (t);
+      return 0;
+    }
+  if (strcmp (s, "max") == 0)
+    {
+      mpfr_setmax (t, mpfr_get_emax ());
+      MPFR_SET_POS (t);
+      return 0;
+    }
+  return mpfr_set_str (t, s, 10, GMP_RNDN);
 }
 
 static void
@@ -80,184 +317,409 @@ tst (void)
   int sv = sizeof (val) / sizeof (*val);
   int i, j;
   int rnd;
-  mpfr_t x, y, z1, z2;
+  mpfr_t x, y, z, tmp;
 
-  mpfr_inits2 (53, x, y, z1, z2, (mpfr_ptr) 0);
+  mpfr_inits2 (53, x, y, z, tmp, (mpfr_ptr) 0);
 
   for (i = 0; i < sv; i++)
     for (j = 0; j < sv; j++)
       RND_LOOP (rnd)
         {
-          int exact, inex1, inex2;
+          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))
+          if (my_setstr (x, val[i]) || my_setstr (y, val[j]))
             {
               printf ("internal error for (%d,%d,%d)\n", i, j, rnd);
               exit (1);
             }
           mpfr_clear_flags ();
-          inex1 = mpfr_pow (z1, x, y, (mp_rnd_t) rnd);
-          if (mpfr_underflow_p ())
-            err ("got underflow", i, j, rnd, z1, inex1);
-          if (mpfr_overflow_p ())
-            err ("got overflow", i, j, rnd, z1, inex1);
-          if (! MPFR_IS_NAN (z1) && mpfr_nanflag_p ())
-            err ("got NaN flag without NaN value", i, j, rnd, z1, inex1);
-          if (MPFR_IS_NAN (z1) && ! mpfr_nanflag_p ())
-            err ("got NaN value without NaN flag", i, j, rnd, z1, inex1);
-          if (inex1 != 0 && ! mpfr_inexflag_p ())
+          inex = mpfr_pow (z, x, y, (mp_rnd_t) rnd);
+          flags = __gmpfr_flags;
+          if (! MPFR_IS_NAN (z) && mpfr_nanflag_p ())
+            err ("got NaN flag without NaN value", i, j, rnd, z, inex);
+          if (MPFR_IS_NAN (z) && ! mpfr_nanflag_p ())
+            err ("got NaN value without NaN flag", i, j, rnd, z, inex);
+          if (inex != 0 && ! mpfr_inexflag_p ())
             err ("got non-zero ternary value without inexact flag",
-                 i, j, rnd, z1, inex1);
-          if (inex1 == 0 && mpfr_inexflag_p ())
+                 i, j, rnd, z, inex);
+          if (inex == 0 && mpfr_inexflag_p ())
             err ("got null ternary value with inexact flag",
-                 i, j, rnd, z1, inex1);
-          exact = MPFR_IS_SINGULAR (z1) ||
-            (mpfr_mul_2ui (z2, z1, 16, GMP_RNDN), mpfr_integer_p (z2));
-          if (exact && inex1 != 0)
-            err ("got exact value with ternary flag different from 0",
-                 i, j, rnd, z1, inex1);
-          if (! exact && inex1 == 0)
-            err ("got inexact value with ternary flag equal to 0",
-                 i, j, rnd, z1, inex1);
+                 i, j, rnd, z, inex);
+          if (i >= 3 && j >= 3)
+            {
+              if (mpfr_underflow_p ())
+                err ("got underflow", i, j, rnd, z, inex);
+              if (mpfr_overflow_p ())
+                err ("got overflow", i, j, rnd, z, inex);
+              exact = MPFR_IS_SINGULAR (z) ||
+                (mpfr_mul_2ui (tmp, z, 16, GMP_RNDN), mpfr_integer_p (tmp));
+              if (exact && inex != 0)
+                err ("got exact value with ternary flag different from 0",
+                     i, j, rnd, z, inex);
+              if (! exact && inex == 0)
+                err ("got inexact value with ternary flag equal to 0",
+                     i, j, rnd, z, inex);
+            }
           if (MPFR_IS_ZERO (x) && ! MPFR_IS_NAN (y) && MPFR_NOTZERO (y))
             {
-              if (MPFR_IS_NEG (y) && ! MPFR_IS_INF (z1))
-                err ("expected an infinity", i, j, rnd, z1, inex1);
-              if (MPFR_IS_POS (y) && ! MPFR_IS_ZERO (z1))
-                err ("expected a zero", i, j, rnd, z1, inex1);
-              if ((MPFR_IS_NEG (x) && is_odd (y)) ^ MPFR_IS_NEG (z1))
-                err ("wrong sign", i, j, rnd, z1, inex1);
+              if (MPFR_IS_NEG (y) && ! MPFR_IS_INF (z))
+                err ("expected an infinity", i, j, rnd, z, inex);
+              if (MPFR_IS_POS (y) && ! MPFR_IS_ZERO (z))
+                err ("expected a zero", i, j, rnd, z, inex);
+              if ((MPFR_IS_NEG (x) && is_odd (y)) ^ MPFR_IS_NEG (z))
+                err ("wrong sign", i, j, rnd, z, inex);
             }
           if (! MPFR_IS_NAN (x) && mpfr_cmp_si (x, -1) == 0)
             {
               /* x = -1 */
               if (! (MPFR_IS_INF (y) || mpfr_integer_p (y)) &&
-                  ! MPFR_IS_NAN (z1))
-                err ("expected NaN", i, j, rnd, z1, inex1);
+                  ! MPFR_IS_NAN (z))
+                err ("expected NaN", i, j, rnd, z, inex);
               if ((MPFR_IS_INF (y) || (mpfr_integer_p (y) && ! is_odd (y)))
-                  && ! mpfr_equal_p (z1, __gmpfr_one))
-                err ("expected 1", i, j, rnd, z1, inex1);
+                  && ! mpfr_equal_p (z, __gmpfr_one))
+                err ("expected 1", i, j, rnd, z, inex);
               if (is_odd (y) &&
-                  (MPFR_IS_NAN (z1) || mpfr_cmp_si (z1, -1) != 0))
-                err ("expected -1", i, j, rnd, z1, inex1);
+                  (MPFR_IS_NAN (z) || mpfr_cmp_si (z, -1) != 0))
+                err ("expected -1", i, j, rnd, z, inex);
             }
           if ((mpfr_equal_p (x, __gmpfr_one) || MPFR_IS_ZERO (y)) &&
-              ! mpfr_equal_p (z1, __gmpfr_one))
-            err ("expected 1", i, j, rnd, z1, inex1);
+              ! mpfr_equal_p (z, __gmpfr_one))
+            err ("expected 1", i, j, rnd, z, inex);
           if (MPFR_IS_PURE_FP (x) && MPFR_IS_NEG (x) &&
               MPFR_IS_FP (y) && ! mpfr_integer_p (y) &&
-              ! MPFR_IS_NAN (z1))
-            err ("expected NaN", i, j, rnd, z1, inex1);
+              ! MPFR_IS_NAN (z))
+            err ("expected NaN", i, j, rnd, z, inex);
           if (MPFR_IS_INF (y) && MPFR_NOTZERO (x))
             {
               int cmpabs1 = mpfr_cmpabs (x, __gmpfr_one);
 
               if ((MPFR_IS_NEG (y) ? (cmpabs1 < 0) : (cmpabs1 > 0)) &&
-                  ! (MPFR_IS_POS (z1) && MPFR_IS_INF (z1)))
-                err ("expected +Inf", i, j, rnd, z1, inex1);
+                  ! (MPFR_IS_POS (z) && MPFR_IS_INF (z)))
+                err ("expected +Inf", i, j, rnd, z, inex);
               if ((MPFR_IS_NEG (y) ? (cmpabs1 > 0) : (cmpabs1 < 0)) &&
-                  ! (MPFR_IS_POS (z1) && MPFR_IS_ZERO (z1)))
-                err ("expected +0", i, j, rnd, z1, inex1);
+                  ! (MPFR_IS_POS (z) && MPFR_IS_ZERO (z)))
+                err ("expected +0", i, j, rnd, z, inex);
             }
           if (MPFR_IS_INF (x) && ! MPFR_IS_NAN (y) && MPFR_NOTZERO (y))
             {
-              if (MPFR_IS_POS (y) && ! MPFR_IS_INF (z1))
-                err ("expected an infinity", i, j, rnd, z1, inex1);
-              if (MPFR_IS_NEG (y) && ! MPFR_IS_ZERO (z1))
-                err ("expected a zero", i, j, rnd, z1, inex1);
-              if ((MPFR_IS_NEG (x) && is_odd (y)) ^ MPFR_IS_NEG (z1))
-                err ("wrong sign", i, j, rnd, z1, inex1);
+              if (MPFR_IS_POS (y) && ! MPFR_IS_INF (z))
+                err ("expected an infinity", i, j, rnd, z, inex);
+              if (MPFR_IS_NEG (y) && ! MPFR_IS_ZERO (z))
+                err ("expected a zero", i, j, rnd, z, inex);
+              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, flags,
+                       "tst");
+        }
+  mpfr_clears (x, y, z, tmp, (mpfr_ptr) 0);
+}
 
-          __gmpfr_flags = MPFR_FLAGS_ALL;
-          inex2 = mpfr_pow (z2, x, y, (mp_rnd_t) rnd);
-          cmpres (i, j, rnd, z1, inex1, z2, inex2, "mpfr_pow, flags set");
+static void
+underflow_up1 (void)
+{
+  mpfr_t delta, x, y, z, z0;
+  mp_exp_t n;
+  int inex;
+  int rnd;
+  int i;
 
-          if (mpfr_integer_p (y))
-            {
-              long yy = mpfr_get_si (y, GMP_RNDN);
+  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);
 
-              /* If y >= 0 and y is not -0, we can test mpfr_pow_ui. */
-              if (MPFR_IS_POS (y))
+  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, "underflow_up1", "mpfr_pow");
+          test_others ("2", sy, (mp_rnd_t) rnd, x, y, z, inex, flags,
+                       "underflow_up1");
+        }
+
+      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 exponent range
+ *   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 (void)
+{
+  mpfr_t x, y, z, z0, eps;
+  mp_exp_t n;
+  int inex;
+  int rnd;
+
+  n = 1 - mpfr_get_emin ();
+  MPFR_ASSERTN (n > 1);
+  if (n > ULONG_MAX)
+    return;
+
+  mpfr_init2 (eps, 2);
+  mpfr_set_ui_2exp (eps, 1, -1, GMP_RNDN);  /* 1/2 */
+  mpfr_div_ui (eps, eps, n, GMP_RNDZ);      /* 1/(2n) rounded toward zero */
+
+  mpfr_init2 (x, sizeof (unsigned long) * CHAR_BIT + 1);
+  inex = mpfr_ui_sub (x, 1, eps, GMP_RNDN);
+  MPFR_ASSERTN (inex == 0);  /* since n < 2^(size_of_long_in_bits) */
+  inex = mpfr_div_2ui (x, x, 1, GMP_RNDN);  /* 1/2 - eps/2 exactly */
+  MPFR_ASSERTN (inex == 0);
+
+  mpfr_init2 (y, sizeof (unsigned long) * CHAR_BIT);
+  inex = mpfr_set_ui (y, n, GMP_RNDN);
+  MPFR_ASSERTN (inex == 0);
+
+  /* 0 < eps < 1 / (2n), thus (1 - eps)^n > 1/2,
+     and 1/2 (1/2)^n < (1/2 - eps/2)^n < (1/2)^n. */
+  mpfr_inits2 (64, z, z0, (mpfr_ptr) 0);
+  RND_LOOP (rnd)
+    {
+      unsigned int ufinex = MPFR_FLAGS_UNDERFLOW | MPFR_FLAGS_INEXACT;
+      int expected_inex;
+      char sy[256];
+
+      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);
+
+      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,
+              "underflow_up2", "mpfr_pow");
+      test_others (NULL, sy, (mp_rnd_t) rnd, x, y, z, inex, ufinex,
+                   "underflow_up2");
+    }
+
+  mpfr_clears (x, y, z, z0, eps, (mpfr_ptr) 0);
+}
+
+static void
+underflow_up3 (void)
+{
+  mpfr_t x, y, z, z0;
+  int inex;
+  int rnd;
+  int i;
+
+  mpfr_init2 (x, 64);
+  mpfr_init2 (y, sizeof (mp_exp_t) * CHAR_BIT);
+  mpfr_init2 (z, 32);
+  mpfr_init2 (z0, 2);
+
+  inex = mpfr_set_exp_t (y, mpfr_get_emin () - 2, GMP_RNDN);
+  MPFR_ASSERTN (inex == 0);
+  for (i = -1; i <= 1; i++)
+    RND_LOOP (rnd)
+      {
+        unsigned int ufinex = MPFR_FLAGS_UNDERFLOW | MPFR_FLAGS_INEXACT;
+        int expected_inex;
+
+        mpfr_set_ui (x, 2, GMP_RNDN);
+        if (i < 0)
+          mpfr_nextbelow (x);
+        if (i > 0)
+          mpfr_nextabove (x);
+        /* x = 2 + i * eps, y = emin - 2, x^y ~= 2^(emin - 2) */
+
+        expected_inex = rnd == GMP_RNDU || (rnd == GMP_RNDN && i < 0) ?
+          1 : -1;
+
+        mpfr_set_ui (z0, 0, GMP_RNDN);
+        if (expected_inex > 0)
+          mpfr_nextabove (z0);
+
+        mpfr_clear_flags ();
+        inex = mpfr_pow (z, x, y, (mp_rnd_t) rnd);
+        cmpres (0, x, "emin - 2", (mp_rnd_t) rnd, z0, expected_inex, z, inex,
+                ufinex, "underflow_up3", "mpfr_pow");
+        test_others (NULL, "emin - 2", (mp_rnd_t) rnd, x, y, z, inex, ufinex,
+                     "underflow_up3");
+      }
+
+  mpfr_clears (x, y, z, z0, (mpfr_ptr) 0);
+}
+
+static void
+underflow_up (void)
+{
+  underflow_up1 ();
+  underflow_up2 ();
+  underflow_up3 ();
+}
+
+static void
+overflow_inv (void)
+{
+  mpfr_t x, y, z;
+  int precx;
+  int s, t;
+  int inex;
+  int rnd;
+
+  mpfr_init2 (y, 2);
+  mpfr_init2 (z, 8);
+
+  mpfr_set_si (y, -1, GMP_RNDN);
+  for (precx = 10; precx <= 100; precx += 90)
+    {
+      char *sp = precx == 10 ?
+        "overflow_inv (precx = 10)" : "overflow_inv (precx = 100)";
+
+      mpfr_init2 (x, precx);
+      for (s = -1; s <= 1; s += 2)
+        {
+          inex = mpfr_set_si_2exp (x, s, - mpfr_get_emax (), GMP_RNDN);
+          MPFR_ASSERTN (inex == 0);
+          for (t = 0; t <= 5; t++)
+            {
+              /* If precx = 10:
+               * x = s * 2^(-emax) * (1 + t * 2^(-9)), so that
+               * 1/x = s * 2^emax * (1 - t * 2^(-9) + eps) with eps > 0.
+               * Values of (1/x) / 2^emax and overflow condition for x > 0:
+               * t = 0: 1                           o: always
+               * t = 1: 0.11111111 100000000011...  o: GMP_RNDN and GMP_RNDU
+               * t = 2: 0.11111111 000000001111...  o: GMP_RNDU
+               * t = 3: 0.11111110 100000100011...  o: never
+               *
+               * If precx = 100:
+               * t = 0: always overflow
+               * t > 0: overflow for GMP_RNDN and GMP_RNDU.
+               */
+              RND_LOOP (rnd)
                 {
-                  MPFR_ASSERTN (yy >= 0);
+                  int inf, overflow;
+
+                  overflow = t == 0 ||
+                    ((mp_rnd_t) rnd == GMP_RNDN && (precx > 10 || t == 1)) ||
+                    ((mp_rnd_t) rnd == (s < 0 ? GMP_RNDD : GMP_RNDU) &&
+                     (precx > 10 || t <= 2));
+                  inf = overflow &&
+                    ((mp_rnd_t) rnd == GMP_RNDN ||
+                     (mp_rnd_t) rnd == (s < 0 ? GMP_RNDD : GMP_RNDU));
                   mpfr_clear_flags ();
-                  inex2 = mpfr_pow_ui (z2, x, yy, (mp_rnd_t) rnd);
-                  cmpres (i, j, rnd, z1, inex1, z2, inex2,
-                          "mpfr_pow_ui, flags cleared");
-                  __gmpfr_flags = MPFR_FLAGS_ALL;
-                  inex2 = mpfr_pow_ui (z2, x, yy, (mp_rnd_t) rnd);
-                  cmpres (i, j, rnd, z1, inex1, z2, inex2,
-                          "mpfr_pow_ui, flags set");
-
-                  /* If x >= 0 and x is not -0, we can test mpfr_ui_pow_ui. */
-                  if (mpfr_integer_p (x) && MPFR_IS_POS (x))
+                  inex = mpfr_pow (z, x, y, (mp_rnd_t) rnd);
+                  if (overflow ^ !! mpfr_overflow_p ())
+                    {
+                      printf ("Bad overflow flag in %s\nfor mpfr_pow%s\n"
+                              "s = %d, t = %d, %s\n", sp,
+                              ext ? ", extended exponent range" : "",
+                              s, t, mpfr_print_rnd_mode ((mp_rnd_t) rnd));
+                      exit (1);
+                    }
+                  if (overflow && (inf ^ !! MPFR_IS_INF (z)))
                     {
-                      unsigned long xx = mpfr_get_ui (x, GMP_RNDN);
-
-                      mpfr_clear_flags ();
-                      inex2 = mpfr_ui_pow_ui (z2, xx, yy, (mp_rnd_t) rnd);
-                      cmpres (i, j, rnd, z1, inex1, z2, inex2,
-                              "mpfr_ui_pow_ui, flags cleared");
-                      __gmpfr_flags = MPFR_FLAGS_ALL;
-                      inex2 = mpfr_ui_pow_ui (z2, xx, yy, (mp_rnd_t) rnd);
-                      cmpres (i, j, rnd, z1, inex1, z2, inex2,
-                              "mpfr_ui_pow_ui, flags set");
+                      printf ("Bad value in %s\nfor mpfr_pow%s\n"
+                              "s = %d, t = %d, %s\nGot ", sp,
+                              ext ? ", extended exponent range" : "",
+                              s, t, mpfr_print_rnd_mode ((mp_rnd_t) rnd));
+                      mpfr_out_str (stdout, 16, 0, z, GMP_RNDN);
+                      printf (" instead of %s value.\n",
+                              inf ? "infinite" : "finite");
+                      exit (1);
                     }
-                }
+                  test_others (NULL, "-1", (mp_rnd_t) rnd, x, y, z,
+                               inex, __gmpfr_flags, sp);
+                }  /* RND_LOOP */
+              mpfr_nexttoinf (x);
+            }  /* for (t = ...) */
+        }  /* for (s = ...) */
+      mpfr_clear (x);
+    }  /* for (precx = ...) */
 
-              /* We can test mpfr_pow_si and mpfr_pow_z when y is not -0. */
-              if (MPFR_IS_POS (y) || MPFR_NOTZERO (y))
-                {
-                  mpz_t yyy;
+  mpfr_clears (y, z, (mpfr_ptr) 0);
+}
 
-                  mpfr_clear_flags ();
-                  inex2 = mpfr_pow_si (z2, x, yy, (mp_rnd_t) rnd);
-                  cmpres (i, j, rnd, z1, inex1, z2, inex2,
-                          "mpfr_pow_si, flags cleared");
-                  __gmpfr_flags = MPFR_FLAGS_ALL;
-                  inex2 = mpfr_pow_si (z2, x, yy, (mp_rnd_t) rnd);
-                  cmpres (i, j, rnd, z1, inex1, z2, inex2,
-                          "mpfr_pow_si, flags set");
-
-                  mpz_init (yyy);
-                  mpfr_get_z (yyy, y, GMP_RNDN);
-                  mpfr_clear_flags ();
-                  inex2 = mpfr_pow_z (z2, x, yyy, (mp_rnd_t) rnd);
-                  cmpres (i, j, rnd, z1, inex1, z2, inex2,
-                          "mpfr_pow_z, flags cleared");
-                  __gmpfr_flags = MPFR_FLAGS_ALL;
-                  inex2 = mpfr_pow_z (z2, x, yyy, (mp_rnd_t) rnd);
-                  cmpres (i, j, rnd, z1, inex1, z2, inex2,
-                          "mpfr_pow_z, flags set");
-                  mpz_clear (yyy);
-                }
-            }
+static void
+alltst (void)
+{
+  mp_exp_t emin, emax;
 
-          /* If x >= 0 and x is not -0, we can test mpfr_ui_pow. */
-          if (mpfr_integer_p (x) && MPFR_IS_POS (x))
-            {
-              unsigned long xx = mpfr_get_ui (x, GMP_RNDN);
+  ext = 0;
+  tst ();
+  underflow_up ();
+  overflow_inv ();
 
-              mpfr_clear_flags ();
-              inex2 = mpfr_ui_pow (z2, xx, y, (mp_rnd_t) rnd);
-              cmpres (i, j, rnd, z1, inex1, z2, inex2,
-                      "mpfr_ui_pow, flags cleared");
-              __gmpfr_flags = MPFR_FLAGS_ALL;
-              inex2 = mpfr_ui_pow (z2, xx, y, (mp_rnd_t) rnd);
-              cmpres (i, j, rnd, z1, inex1, z2, inex2,
-                      "mpfr_ui_pow, flags set");
-            }
-        }
-  mpfr_clears (x, y, z1, z2, (mpfr_ptr) 0);
+  emin = mpfr_get_emin ();
+  emax = mpfr_get_emax ();
+  set_emin (MPFR_EMIN_MIN);
+  set_emax (MPFR_EMAX_MAX);
+  if (mpfr_get_emin () != emin || mpfr_get_emax () != emax)
+    {
+      ext = 1;
+      tst ();
+      underflow_up ();
+      overflow_inv ();
+      set_emin (emin);
+      set_emax (emax);
+    }
 }
 
 int
-main (void)
+main (int argc, char *argv[])
 {
   tests_start_mpfr ();
-  tst ();
+  all_cmpres_errors = argc > 1;
+  alltst ();
   tests_end_mpfr ();
-  return 0;
+  return all_cmpres_errors < 0;
 }
diff --git a/tests/tpow_z.c b/tests/tpow_z.c
index 7ece9e8a0..08f7b455a 100644
--- a/tests/tpow_z.c
+++ b/tests/tpow_z.c
@@ -295,6 +295,29 @@ check_overflow (void)
   mpz_clear (z);
 }
 
+/* bug reported by Carl Witty (32-bit architecture) */
+static void
+bug20080223 (void)
+{
+  mpfr_t a, exp, answer;
+
+  mpfr_init2 (a, 53);
+  mpfr_init2 (exp, 53);
+  mpfr_init2 (answer, 53);
+
+  mpfr_set_si (exp, -1073741824, GMP_RNDN);
+
+  mpfr_set_str (a, "1.999999999", 10, GMP_RNDN);
+  /* a = 562949953139837/2^48 */
+  mpfr_pow (answer, a, exp, GMP_RNDN);
+  mpfr_set_str_binary (a, "0.110110101111011001110000111111100011101000111011101E-1073741823");
+  MPFR_ASSERTN(mpfr_cmp0 (answer, a) == 0);
+
+  mpfr_clear (a);
+  mpfr_clear (exp);
+  mpfr_clear (answer);
+}
+
 int
 main (void)
 {
@@ -304,6 +327,7 @@ main (void)
   check_integer (2, 163, 100);
   check_regression ();
   bug20071104 ();
+  bug20080223 ();
   check_overflow ();
 
   tests_end_mpfr ();