--- pow.c	Mon Apr 15 14:51:05 2002
+++ pow.c.new	Fri Oct 10 11:10:22 2003
@@ -21,12 +21,110 @@
 
 #include "gmp.h"
 #include "gmp-impl.h"
+#include "longlong.h"
 #include "mpfr.h"
 #include "mpfr-impl.h"
 
+static int mpfr_pow_is_exact _PROTO((mpfr_srcptr, mpfr_srcptr));
+
+/* return non zero iff x^y is exact.
+   Assumes x and y are ordinary numbers (neither NaN nor Inf),
+   and y is not zero.
+*/
+int
+mpfr_pow_is_exact (mpfr_srcptr x, mpfr_srcptr y)
+{
+  mp_exp_t d;
+  unsigned long i, c;
+  mp_limb_t *yp;
+  
+  if ((mpfr_sgn (x) < 0) && (mpfr_isinteger (y) == 0))
+      return 0;
+
+  if (mpfr_sgn (y) < 0)
+    return mpfr_cmp_si_2exp (x, MPFR_SIGN(x), MPFR_EXP(x) - 1) == 0;
+
+  /* compute d such that y = c*2^d with c odd integer */
+  d = MPFR_EXP(y) - MPFR_PREC(y);
+  /* since y is not zero, necessarily one of the mantissa limbs is not zero,
+     thus we can simply loop until we find a non zero limb */
+  yp = MPFR_MANT(y);
+  for (i = 0; yp[i] == 0; i++, d += BITS_PER_MP_LIMB);
+  /* now yp[i] is not zero */
+  count_trailing_zeros (c, yp[i]);
+  d += c;
+  
+  if (d < 0)
+    {
+      mpz_t a;
+      mp_exp_t b;
+
+      mpz_init (a);
+      b = mpfr_get_z_exp (a, x); /* x = a * 2^b */
+      c = mpz_scan1 (a, 0);
+      mpz_div_2exp (a, a, c);
+      b += c;
+      /* now a is odd */
+      while (d != 0)
+        {
+          if (mpz_perfect_square_p (a))
+            {
+              d++;
+              mpz_sqrt (a, a);
+            }
+          else
+            {
+              mpz_clear (a);
+              return 0;
+            }
+        }
+      mpz_clear (a);
+    }
+
+    return 1;
+}
+
+/* Return 1 if y is an odd integer, 0 otherwise. */
+static int
+is_odd (mpfr_srcptr y)
+{
+  mp_exp_t expo;
+  mp_prec_t prec;
+  mp_size_t yn;
+  mp_limb_t *yp;
+
+  MPFR_ASSERTN(MPFR_IS_FP(y));
+
+  if (MPFR_IS_ZERO(y))
+    return 0;
+
+  expo = MPFR_EXP (y);
+  if (expo <= 0)
+    return 0;  /* |y| < 1 and not 0 */
+
+  prec = MPFR_PREC(y);
+  if (expo > prec)
+    return 0;  /* y is a multiple of 2^(expo-prec), thus not odd */
+
+  /* 0 < expo <= prec */
+
+  yn = (prec - 1) / BITS_PER_MP_LIMB;  /* index of last limb */
+  yn -= (mp_size_t) (expo / BITS_PER_MP_LIMB);
+  MPFR_ASSERTN(yn >= 0);
+  /* now the index of the last limb containing bits of the fractional part */
+
+  yp = MPFR_MANT(y);
+  if (expo % BITS_PER_MP_LIMB == 0 ? (yp[yn+1] & 1) == 0 || yp[yn] != 0
+      : yp[yn] << ((expo % BITS_PER_MP_LIMB) - 1) != GMP_LIMB_HIGHBIT)
+    return 0;
+  while (--yn >= 0)
+    if (yp[yn] != 0)
+      return 0;
+  return 1;
+}
+
 /* The computation of z = pow(x,y) is done by
    z = exp(y * log(x)) = x^y */
-
 int
 mpfr_pow (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd_mode)
 {
@@ -129,15 +227,30 @@
 
   if (MPFR_IS_ZERO(y))
     {
-      return mpfr_set_ui(z,1,GMP_RNDN);
+      return mpfr_set_ui (z, 1, GMP_RNDN);
     }
 
-  if(mpfr_isinteger(y))
+  /* detect overflows: |x^y| >= 2^EMAX when (EXP(x)-1) * y >= EMAX for y > 0,
+                                       or   EXP(x) * y     >= EMAX for y < 0 */
+  {
+    double exy;
+    int negative;
+
+    exy = (double) (mpfr_sgn (y) > 0) ? MPFR_EXP(x) - 1 : MPFR_EXP(x);
+    exy *= mpfr_get_d (y, GMP_RNDZ);
+    if (exy >= (double) __mpfr_emax)
+      {
+        negative = MPFR_SIGN(x) < 0 && is_odd (y);
+        return mpfr_set_overflow (z, rnd_mode, negative ? -1 : 1);
+      }
+  }  
+
+  if (mpfr_isinteger (y))
     {
       mpz_t zi;
       long int zii;
       int exptol;
-    
+
       mpz_init(zi);  
       exptol = mpfr_get_z_exp (zi, y);
         
@@ -149,7 +262,7 @@
       zii=mpz_get_ui(zi);
         
       mpz_clear(zi);
-      return mpfr_pow_si(z,x,zii,rnd_mode); 
+      return mpfr_pow_si (z, x, zii, rnd_mode); 
     }
 
   if (MPFR_IS_INF(x))
@@ -191,6 +304,7 @@
   {
     /* Declaration of the intermediary variable */
       mpfr_t t, te, ti;
+      int loop = 0, ok;
 
       /* Declaration of the size variable */
       mp_prec_t Nx = MPFR_PREC(x);   /* Precision of input variable */
@@ -209,29 +323,41 @@
       mpfr_init(ti);
       mpfr_init(te);             
 
-      do {
+      do
+        {
+
+          loop ++;
 
-	/* reactualisation of the precision */
-	mpfr_set_prec(t,Nt);                    
-	mpfr_set_prec(ti,Nt);                    
-   	mpfr_set_prec(te,Nt);             
-
-	/* compute   exp(y*ln(x))*/
-        mpfr_log(ti,x,GMP_RNDU);         /* ln(n) */
-        mpfr_mul(te,y,ti,GMP_RNDU);       /* y*ln(n) */
-        mpfr_exp(t,te,GMP_RNDN);         /* exp(x*ln(n))*/
+          /* reactualisation of the precision */
+          mpfr_set_prec (t, Nt);
+          mpfr_set_prec (ti, Nt);
+          mpfr_set_prec (te, Nt);
+
+          /* compute exp(y*ln(x)) */
+          mpfr_log (ti, x, GMP_RNDU);         /* ln(n) */
+          mpfr_mul (te, y, ti, GMP_RNDU);       /* y*ln(n) */
+          mpfr_exp (t, te, GMP_RNDN);         /* exp(x*ln(n))*/
 
 	/* estimation of the error -- see pow function in algorithms.ps*/
-	err = Nt - (MPFR_EXP(te)+3);
+          err = Nt - (MPFR_EXP(te) + 3);
 
 	/* actualisation of the precision */
-	Nt += 10;
+          Nt += 10;
 
-      } while (err<0 || !mpfr_can_round(t,err,GMP_RNDN,rnd_mode,Ny));
-      inexact = mpfr_set(z,t,rnd_mode);
-      mpfr_clear(t);
-      mpfr_clear(ti);
-      mpfr_clear(te);
+          ok = mpfr_can_round (t, err, GMP_RNDN, rnd_mode, Ny);
+
+          /* check exact power */
+          if (ok == 0 && loop == 1)
+            ok = mpfr_pow_is_exact (x, y);
+
+        }
+      while (err < 0 || ok == 0);
+      
+      inexact = mpfr_set (z, t, rnd_mode);
+
+      mpfr_clear (t);
+      mpfr_clear (ti);
+      mpfr_clear (te);
     }
     return inexact;
 }