Fix overflows problems.

Clean up overflow handling.
Maybe some bugs remain...


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

3 files changed, 103 insertions, 68 deletions

diff --git a/pow_si.c b/pow_si.c
index db201ca86..6b904ac5b 100644
--- a/pow_si.c
+++ b/pow_si.c
@@ -70,7 +70,16 @@ mpfr_pow_si (mpfr_ptr y, mpfr_srcptr x, long int n, mp_rnd_t rnd_mode)
           mp_exp_t expx = MPFR_EXP (x); /* warning: x and y may be the same
                                             variable */
           mpfr_set_si (y, (n % 2) ? MPFR_INT_SIGN(x) : 1, rnd_mode);
-          MPFR_EXP(y) += n * (expx - 1);
+	  expx --;
+	  MPFR_ASSERTD (n < 0);
+	  /* Warning n*expx may overflow! 
+	     Many systems abort with LONG_MIN / 1 or LONG_MIN/-1*/
+	  if (n != -1 && expx > 0 && -expx < MPFR_EXP_MIN / (-n))
+	    MPFR_EXP (y) = MPFR_EMIN_MIN - 1; /* Underflow */
+	  else if (n != -1 && expx < 0 && -expx > MPFR_EXP_MAX / (-n))
+	    MPFR_EXP (y) = MPFR_EMAX_MAX + 1; /* Overflow */
+          else
+	    MPFR_EXP(y) += n * expx;
           return mpfr_check_range (y, 0, rnd_mode);
         }
 
@@ -79,7 +88,7 @@ mpfr_pow_si (mpfr_ptr y, mpfr_srcptr x, long int n, mp_rnd_t rnd_mode)
       /* General case */
       {
         /* Declaration of the intermediary variable */
-        mpfr_t t, ti;
+        mpfr_t t;
 
         /* Declaration of the size variable */
         mp_prec_t Nx = MPFR_PREC(x);   /* Precision of input variable */
@@ -99,32 +108,27 @@ mpfr_pow_si (mpfr_ptr y, mpfr_srcptr x, long int n, mp_rnd_t rnd_mode)
 
         /* initialise of intermediary	variable */
         mpfr_init2 (t, Nt);
-        mpfr_init2 (ti, Nt);
 
-        do
-          {
-            /* reactualisation of the precision */
-            mpfr_set_prec (t, Nt);
-            mpfr_set_prec (ti, Nt);
-
-            /* compute 1/(x^n) n>0*/
-            mpfr_pow_ui (ti, x, (unsigned long int) n, GMP_RNDN);
-            mpfr_ui_div (t, 1, ti, GMP_RNDN);
-
-            /* error estimate -- see pow function in algorithms.ps */
-            err = Nt - 3;
-
-            /* actualisation of the precision */
-            Nt += 10;
-          }
-        while (err < 0 || (!mpfr_can_round (t, err, GMP_RNDN, GMP_RNDZ,
-					    Ny + (rnd_mode == GMP_RNDN))
-			   /* An overflow can occurs, producing an underflow */
-			   && !MPFR_IS_ZERO(t) ));
+        do {
+	  /* reactualisation of the precision */
+	  mpfr_set_prec (t, Nt);
+	  
+	  /* compute 1/(x^n) n>0*/
+	  mpfr_pow_ui (t, x, (unsigned long int) n, GMP_RNDN);
+	  inexact = MPFR_IS_ZERO (t) || MPFR_IS_INF (t);
+	  mpfr_ui_div (t, 1, t, GMP_RNDN);
+	  inexact = inexact || MPFR_IS_ZERO (t) || MPFR_IS_INF (t);
+	  /* error estimate -- see pow function in algorithms.ps */
+	  err = Nt - 3;
+	  
+	  /* actualisation of the precision */
+	  Nt += BITS_PER_MP_LIMB;
+	} while (inexact == 0 &&
+		 !mpfr_can_round (t, err, GMP_RNDN, GMP_RNDZ,
+				  Ny + (rnd_mode == GMP_RNDN)));
 
         inexact = mpfr_set (y, t, rnd_mode);
         mpfr_clear (t);
-        mpfr_clear (ti);
 	MPFR_SAVE_EXPO_FREE (expo);
         return mpfr_check_range (y, inexact, rnd_mode);
       }
diff --git a/pow_ui.c b/pow_ui.c
index 91f2be7ec..875ee59c8 100644
--- a/pow_ui.c
+++ b/pow_ui.c
@@ -61,6 +61,10 @@ mpfr_pow_ui (mpfr_ptr x, mpfr_srcptr y, unsigned long int n, mp_rnd_t rnd)
           MPFR_ASSERTD (MPFR_IS_ZERO (y));
           /* 0^n = 0 for any n */
 	  MPFR_SET_ZERO (x);
+	  if (MPFR_IS_POS (y) || ((n & 1) == 0))
+            MPFR_SET_POS (x);
+          else
+            MPFR_SET_NEG (x);
 	  MPFR_RET (0);
 	}
     }
@@ -79,58 +83,56 @@ mpfr_pow_ui (mpfr_ptr x, mpfr_srcptr y, unsigned long int n, mp_rnd_t rnd)
   /* MPFR_CLEAR_FLAGS useless due to mpfr_set */
 
   MPFR_SAVE_EXPO_MARK (expo);
+  __gmpfr_emin -= 3;  /* So that we can check for underflow properly */
 
   prec = MPFR_PREC (x);
   mpfr_init2 (res, prec + 9);
 
   rnd1 = MPFR_IS_POS (y) ? GMP_RNDU : GMP_RNDD; /* away */
 
-  do
-    {
-      int i;
-
-      prec += 3;
-      for (m = n, i = 0; m; i++, m >>= 1)
-        prec++;
-      /* now 2^(i-1) <= n < 2^i */
-      mpfr_set_prec (res, prec);
-      err = prec <= (mpfr_prec_t) i ? 0 : prec - (mpfr_prec_t) i;
-      MPFR_ASSERTD (i >= 1);
-      /* First step: compute square from y */
-      inexact = mpfr_sqr (res, y, GMP_RNDU);
-      if (n & (1UL << (i-2)))
-	inexact |= mpfr_mul (res, res, y, rnd1);
-      for (i -= 3; i >= 0; i--)
-        {
-          inexact |= mpfr_sqr (res, res, GMP_RNDU);
-          if (n & (1UL << i))
-            inexact |= mpfr_mul (res, res, y, rnd1);
-        }
-
-      /* Check Overflow (can't use MPFR_GET_EXP since there can be an INF )*/
-      if (MPFR_UNLIKELY (MPFR_IS_INF (res) 
-			 || MPFR_EXP (res) >= __gmpfr_emax))
-	{
-          mpfr_clear (res);
-          MPFR_SAVE_EXPO_FREE (expo);
-          return mpfr_set_overflow (x, rnd, 
-				    (n % 2) ? MPFR_SIGN (y) : MPFR_SIGN_POS);
-	}
-      /* Check Underflow (can't use MPFR_GET_EXP since there can be a ZERO )*/
-      if (MPFR_UNLIKELY (MPFR_IS_ZERO (res)
-			 || MPFR_EXP (res) <= __gmpfr_emin))
-        {
-          mpfr_clear (res);
-	  MPFR_SAVE_EXPO_FREE (expo);
-          return mpfr_set_underflow (x, rnd, 
-				     (n % 2) ? MPFR_SIGN(y) : MPFR_SIGN_POS);
-        }
-    }
-  while (inexact && !mpfr_can_round (res, err, GMP_RNDN, GMP_RNDZ,
-                                     MPFR_PREC(x) + (rnd == GMP_RNDN)));
-
+  do {
+    int i;
+    prec += 3;
+    for (m = n, i = 0; m; i++, m >>= 1)
+      prec++;
+    /* now 2^(i-1) <= n < 2^i */
+    mpfr_set_prec (res, prec);
+    err = prec <= (mpfr_prec_t) i ? 0 : prec - (mpfr_prec_t) i;
+    MPFR_ASSERTD (i >= 1);
+    mpfr_clear_overflow ();
+    mpfr_clear_underflow ();
+    /* First step: compute square from y */
+    inexact = mpfr_sqr (res, y, GMP_RNDU);
+    if (n & (1UL << (i-2)))
+      inexact |= mpfr_mul (res, res, y, rnd1);
+    for (i -= 3; i >= 0 && !mpfr_overflow_p () && !mpfr_underflow_p (); i--) 
+      {
+	inexact |= mpfr_sqr (res, res, GMP_RNDU);
+	if (n & (1UL << i))
+	  inexact |= mpfr_mul (res, res, y, rnd1);
+      }
+  }
+  while (inexact && !mpfr_overflow_p () && !mpfr_underflow_p ()
+	 && !mpfr_can_round (res, err, GMP_RNDN, GMP_RNDZ,
+			     MPFR_PREC(x) + (rnd == GMP_RNDN)));
   inexact = mpfr_set (x, res, rnd);
   mpfr_clear (res);
+
+  /* Check Overflow */
+  if (MPFR_UNLIKELY (mpfr_overflow_p ())) {
+    MPFR_SAVE_EXPO_FREE (expo);
+    return mpfr_set_overflow (x, rnd,
+			      (n % 2) ? MPFR_SIGN (y) : MPFR_SIGN_POS);
+  }
+  /* Check Underflow  */
+  else if (MPFR_UNLIKELY (mpfr_underflow_p ()))
+    {
+      if (rnd == GMP_RNDN)
+	rnd = GMP_RNDZ;
+      MPFR_SAVE_EXPO_FREE (expo);
+      return mpfr_set_underflow (x, rnd,
+				 (n % 2) ? MPFR_SIGN(y) : MPFR_SIGN_POS);
+    }
   MPFR_SAVE_EXPO_FREE (expo);
   return mpfr_check_range (x, inexact, rnd);
 }
diff --git a/tests/tpow.c b/tests/tpow.c
index d31f90351..510575632 100644
--- a/tests/tpow.c
+++ b/tests/tpow.c
@@ -93,6 +93,35 @@ check_pow_ui (void)
       exit (1);
     }
 
+  /* Check 0 */
+  MPFR_SET_ZERO (a);
+  MPFR_SET_POS (a);
+  mpfr_set_si (b, -1, GMP_RNDN);
+  res = mpfr_pow_ui (b, a, 1, GMP_RNDN);
+  if (res != 0 || MPFR_IS_NEG (b))
+    {
+      printf ("Error for (0+)^1\n");
+      exit (1);
+    }
+  MPFR_SET_ZERO (a);
+  MPFR_SET_NEG (a);
+  mpfr_set_ui (b, 1, GMP_RNDN);
+  res = mpfr_pow_ui (b, a, 5, GMP_RNDN);
+  if (res != 0 || MPFR_IS_POS (b))
+    {
+      printf ("Error for (0-)^5\n");
+      exit (1);
+    }
+  MPFR_SET_ZERO (a);
+  MPFR_SET_NEG (a);
+  mpfr_set_si (b, -1, GMP_RNDN);
+  res = mpfr_pow_ui (b, a, 6, GMP_RNDN);
+  if (res != 0 || MPFR_IS_NEG (b))
+    {
+      printf ("Error for (0-)^6\n");
+      exit (1);
+    }
+
   mpfr_clear (a);
   mpfr_clear (b);
 }