Merged -r4319:4320 -r4331:4332 -r4337:4343 -r4344:4345 -r4355:4358

from the trunk: Fixed several special cases in pow_si.c and added
corresponding testcases in tests/tpow.c. NaN could be returned for
pow_si(2, LONG_MIN+1), and there were integer overflows for some
input values (but it seems that results were correct as long as the
implementation (C compiler) guaranteed two's complement wrapping).
Moreover, as a consequence of the above fixes, this patch avoids a
bug in Sun's compiler for Solaris/x86 (with optimizations).


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

2 files changed, 99 insertions, 21 deletions

diff --git a/pow_si.c b/pow_si.c
index 35f3e7019..5b9b894ef 100644
--- a/pow_si.c
+++ b/pow_si.c
@@ -68,24 +68,20 @@ mpfr_pow_si (mpfr_ptr y, mpfr_srcptr x, long int n, mp_rnd_t rnd)
       /* detect exact powers: x^(-n) is exact iff x is a power of 2 */
       if (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_set_si (y, (n % 2) ? MPFR_INT_SIGN(x) : 1, rnd);
-          expx --;
+          mp_exp_t expx = MPFR_EXP (x) - 1, expy;
           MPFR_ASSERTD (n < 0);
-          /* Warning n*expx may overflow!
-             Some 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);
+          /* Warning: n * expx may overflow!
+             Some systems (apparently alpha-freebsd) abort with
+             LONG_MIN / 1, and LONG_MIN / -1 is undefined. */
+          expy =
+            n != -1 && expx > 0 && expx > (__gmpfr_emin - 1) / n ?
+            MPFR_EMIN_MIN - 2 /* Underflow */ :
+            n != -1 && expx < 0 && expx < (__gmpfr_emax - 1) / n ?
+            MPFR_EMAX_MAX /* Overflow */ : n * expx;
+          return mpfr_set_si_2exp (y, n % 2 ? MPFR_INT_SIGN (x) : 1,
+                                   expy, rnd);
         }
 
-      n = -n;
-
       /* General case */
       {
         /* Declaration of the intermediary variable */
@@ -95,9 +91,12 @@ mpfr_pow_si (mpfr_ptr y, mpfr_srcptr x, long int n, mp_rnd_t rnd)
         mp_prec_t Nt;                              /* working precision */
         mp_exp_t  err;                             /* error */
         int inexact;
+        unsigned long abs_n;
         MPFR_SAVE_EXPO_DECL (expo);
         MPFR_ZIV_DECL (loop);
 
+        abs_n = - (unsigned long) n;
+
         /* compute the precision of intermediary variable */
         /* the optimal number of bits : see algorithms.tex */
         Nt = Ny + 3 + MPFR_INT_CEIL_LOG2 (Ny);
@@ -110,17 +109,17 @@ mpfr_pow_si (mpfr_ptr y, mpfr_srcptr x, long int n, mp_rnd_t rnd)
         MPFR_ZIV_INIT (loop, Nt);
         for (;;)
           {
-            /* compute 1/(x^n) n>0*/
-            mpfr_pow_ui (t, x, (unsigned long int) n, GMP_RNDN);
+            /* 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. */
+            /* FIXME: old code improved, but I think this is still incorrect. */
             if (MPFR_UNLIKELY (MPFR_IS_ZERO (t)))
               {
                 MPFR_ZIV_FREE (loop);
                 mpfr_clear (t);
                 MPFR_SAVE_EXPO_FREE (expo);
                 return mpfr_underflow (y, rnd == GMP_RNDN ? GMP_RNDZ : rnd,
-                                       (unsigned) n & 1 ? MPFR_SIGN (x) :
+                                       abs_n & 1 ? MPFR_SIGN (x) :
                                        MPFR_SIGN_POS);
               }
             if (MPFR_UNLIKELY (MPFR_IS_INF (t)))
@@ -128,8 +127,7 @@ mpfr_pow_si (mpfr_ptr y, mpfr_srcptr x, long int n, mp_rnd_t rnd)
                 MPFR_ZIV_FREE (loop);
                 mpfr_clear (t);
                 MPFR_SAVE_EXPO_FREE (expo);
-                return mpfr_overflow (y, rnd,
-                                      (unsigned) n & 1 ? MPFR_SIGN (x) :
+                return mpfr_overflow (y, rnd, abs_n & 1 ? MPFR_SIGN (x) :
                                       MPFR_SIGN_POS);
               }
             /* error estimate -- see pow function in algorithms.ps */
diff --git a/tests/tpow.c b/tests/tpow.c
index d1abd846a..6bd9eee01 100644
--- a/tests/tpow.c
+++ b/tests/tpow.c
@@ -24,6 +24,7 @@ MA 02110-1301, USA. */
 #include <stdlib.h>
 #include <float.h>
 #include <math.h>
+#include <limits.h>
 
 #include "mpfr-test.h"
 
@@ -232,6 +233,55 @@ check_pow_si (void)
   mpfr_pow_si (x, x, -2, GMP_RNDN);
   MPFR_ASSERTN(mpfr_inf_p (x) && mpfr_sgn (x) > 0);
 
+  mpfr_set_si (x, 2, GMP_RNDN);
+  mpfr_pow_si (x, x, LONG_MAX, GMP_RNDN);  /* 2^LONG_MAX */
+  if (LONG_MAX > mpfr_get_emax () - 1)  /* LONG_MAX + 1 > emax */
+    {
+      MPFR_ASSERTN (mpfr_inf_p (x));
+    }
+  else
+    {
+      MPFR_ASSERTN (mpfr_cmp_si_2exp (x, 1, LONG_MAX));
+    }
+
+  mpfr_set_si (x, 2, GMP_RNDN);
+  mpfr_pow_si (x, x, LONG_MIN, GMP_RNDN);  /* 2^LONG_MIN */
+  if (LONG_MIN + 1 < mpfr_get_emin ())
+    {
+      MPFR_ASSERTN (mpfr_zero_p (x));
+    }
+  else
+    {
+      MPFR_ASSERTN (mpfr_cmp_si_2exp (x, 1, LONG_MIN));
+    }
+
+  mpfr_set_si (x, 2, GMP_RNDN);
+  mpfr_pow_si (x, x, LONG_MIN + 1, GMP_RNDN);  /* 2^(LONG_MIN+1) */
+  if (mpfr_nan_p (x))
+    {
+      printf ("Error in pow_si(2, LONG_MIN+1): got NaN\n");
+      exit (1);
+    }
+  if (LONG_MIN + 2 < mpfr_get_emin ())
+    {
+      MPFR_ASSERTN (mpfr_zero_p (x));
+    }
+  else
+    {
+      MPFR_ASSERTN (mpfr_cmp_si_2exp (x, 1, LONG_MIN + 1));
+    }
+
+  mpfr_set_si_2exp (x, 1, -1, GMP_RNDN);  /* 0.5 */
+  mpfr_pow_si (x, x, LONG_MIN, GMP_RNDN);  /* 2^(-LONG_MIN) */
+  if (LONG_MIN < 1 - mpfr_get_emax ())  /* 1 - LONG_MIN > emax */
+    {
+      MPFR_ASSERTN (mpfr_inf_p (x));
+    }
+  else
+    {
+      MPFR_ASSERTN (mpfr_cmp_si_2exp (x, 2, - (LONG_MIN + 1)));
+    }
+
   mpfr_clear (x);
 }
 
@@ -271,6 +321,35 @@ check_special_pow_si ()
 }
 
 static void
+pow_si_long_min (void)
+{
+  mpfr_t x, y, z;
+  int inex;
+
+  mpfr_inits2 (sizeof(long) * CHAR_BIT + 32, x, y, z, (void *) 0);
+  mpfr_set_si_2exp (x, 3, -1, GMP_RNDN);  /* 1.5 */
+
+  inex = mpfr_set_si (y, LONG_MIN, GMP_RNDN);
+  MPFR_ASSERTN (inex == 0);
+  mpfr_nextbelow (y);
+  mpfr_pow (y, x, y, GMP_RNDD);
+
+  inex = mpfr_set_si (z, LONG_MIN, GMP_RNDN);
+  MPFR_ASSERTN (inex == 0);
+  mpfr_nextabove (z);
+  mpfr_pow (z, x, z, GMP_RNDU);
+
+  mpfr_pow_si (x, x, LONG_MIN, GMP_RNDN);  /* 1.5^LONG_MIN */
+  if (mpfr_cmp (x, y) < 0 || mpfr_cmp (x, z) > 0)
+    {
+      printf ("Error in pow_si_long_min\n");
+      exit (1);
+    }
+
+  mpfr_clears (x, y, z, (void *) 0);
+}
+
+static void
 check_inexact (mp_prec_t p)
 {
   mpfr_t x, y, z, t;
@@ -752,6 +831,7 @@ main (void)
   check_pow_ui ();
   check_pow_si ();
   check_special_pow_si ();
+  pow_si_long_min ();
   for (p = 2; p < 100; p++)
     check_inexact (p);
   underflows ();