Really reverted to rev. 1.54 + kept optimization + commented out

is_odd_even.


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

1 files changed, 63 insertions, 39 deletions

diff --git a/pow.c b/pow.c
index 5917591d3..e40a502bd 100644
--- a/pow.c
+++ b/pow.c
@@ -1,4 +1,4 @@
-/* mpfr_pow -- power function x^y 
+/* mpfr_pow -- power function x^y
 
 Copyright 2001, 2002, 2003, 2004 Free Software Foundation, Inc.
 
@@ -35,11 +35,6 @@ mpfr_pow_is_exact (mpfr_srcptr x, mpfr_srcptr y)
   mp_limb_t *yp;
   mp_size_t ysize;
 
-#if 0 /* useless since already checked in mpfr_pow */ 
-  if ((mpfr_sgn (x) < 0) && (mpfr_integer_p (y) == 0))
-    return 0;
-#endif
-
   if (mpfr_sgn (y) < 0)
     return mpfr_cmp_si_2exp (x, MPFR_SIGN(x), MPFR_GET_EXP (x) - 1) == 0;
 
@@ -53,7 +48,7 @@ mpfr_pow_is_exact (mpfr_srcptr x, mpfr_srcptr y)
   /* now yp[i] is not zero */
   count_trailing_zeros (c, yp[i]);
   d += c;
-  
+
   if (d < 0)
     {
       mpz_t a;
@@ -84,12 +79,13 @@ mpfr_pow_is_exact (mpfr_srcptr x, mpfr_srcptr y)
     return 1;
 }
 
+#if 0
 /* Return 1 if y is an odd integer, -1 if an even integer,
    0 otherwise (not an integer).
    Assumes y is not zero.
 */
 static int
-is_odd (mpfr_srcptr y)
+is_odd_even (mpfr_srcptr y)
 {
   mp_exp_t expo;
   mp_prec_t prec;
@@ -137,10 +133,50 @@ is_odd (mpfr_srcptr y)
       return 0; /* not an integer */
   return res;
 }
+#endif
+
+/* 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_GET_EXP (y);
+  if (expo <= 0)
+    return 0;  /* |y| < 1 and not 0 */
+
+  prec = MPFR_PREC(y);
+  if ((mpfr_prec_t) 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) != MPFR_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
-   For the special cases, see Section F.9.4.4 of the C99 standard:
+   For the special cases, see Section F.9.4.4 of the C standard:
      _ pow(±0, y) = ±inf for y an odd integer < 0.
      _ pow(±0, y) = +inf for y < 0 and not an odd integer.
      _ pow(±0, y) = ±0 for y an odd integer > 0.
@@ -169,13 +205,11 @@ mpfr_pow (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd_mode)
       /* pow(x, 0) returns 1 for any x, even a NaN. */
       if (MPFR_UNLIKELY( MPFR_IS_ZERO(y) ))
 	return mpfr_set_ui (z, 1, GMP_RNDN);
-      /* now x is NaN, Inf, or zero, or y is NaN or Inf */
       else if (MPFR_IS_NAN(x))
 	{
 	  MPFR_SET_NAN(z);
 	  MPFR_RET_NAN;
 	}
-      /* now x is Inf or zero, or y is NaN or Inf */
       else if (MPFR_IS_NAN(y))
 	{
 	  /* pow(+1, NaN) returns 1. */
@@ -184,7 +218,6 @@ mpfr_pow (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd_mode)
 	  MPFR_SET_NAN(z);
 	  MPFR_RET_NAN;
 	}
-      /* now x is Inf or zero, or y is Inf */
       else if (MPFR_IS_INF(y))
 	{
 	  if (MPFR_IS_INF(x))
@@ -226,40 +259,29 @@ mpfr_pow (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd_mode)
 		}
 	    }
 	}
-      /* x is Inf or zero, and y is now an ordinary number (non-zero) */
       MPFR_ASSERTD(MPFR_IS_FP(y));
       if (MPFR_IS_INF(x))
 	{
-          /* +Inf^y gives +Inf if y > 0, +0 if y < 0
-             -Inf^y gives -Inf for y an odd integer > 0           (a)
-                          +Inf for y > 0 and not an odd integer   (b)
-                          -0 for y an odd integer < 0             (c)
-                          +0 for y < 0 and not an odd integer     (d).
-             Warning: these are the rules of the C99 standard, which
-             may be counter-intuitive, especially (b) and (d), where
-             we might expect NaN when taking the limit of x^y when x -> -Inf.
-          */
 	  int negative;
 	  /* Determine the sign now, in case y and z are the same object */
-	  negative = MPFR_IS_NEG(x) && is_odd (y) == 1;
+	  negative = MPFR_IS_NEG(x) && is_odd(y);
 	  MPFR_CLEAR_FLAGS(z);
-          if (MPFR_IS_POS(y))
-            MPFR_SET_INF(z);
-          else
-            MPFR_SET_ZERO(z);
-          if (negative) /* x = -Inf and y is an odd integer */
-            MPFR_SET_NEG(z);
-          else
-            MPFR_SET_POS(z);
+	  if (MPFR_IS_POS(y))
+	    MPFR_SET_INF(z);
+	  else
+	    MPFR_SET_ZERO(z);
+	  if (negative)
+	    MPFR_SET_NEG(z);
+	  else
+	    MPFR_SET_POS(z);
 	  MPFR_RET(0);
 	}
-      /* x is zero, y is an ordinary (non-zero) number */
       MPFR_ASSERTD(MPFR_IS_FP(x));
-      MPFR_ASSERTD(MPFR_IS_ZERO(x));
+      if (MPFR_IS_ZERO(x))
 	{
 	  int negative;
 	  /* Determine the sign now, in case y and z are the same object */
-	  negative = MPFR_IS_NEG(x) && is_odd (y) == 1;
+	  negative = MPFR_IS_NEG(x) && is_odd (y);
 	  MPFR_CLEAR_FLAGS(z);
 	  if (MPFR_IS_NEG(y))
 	    MPFR_SET_INF(z);
@@ -271,6 +293,9 @@ mpfr_pow (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd_mode)
 	    MPFR_SET_POS(z);
 	  MPFR_RET(0);
 	}
+      /* Should never reach this code */
+      MPFR_ASSERTN(0);
+      return 0;
     }
 
   if (mpfr_cmp_ui (x, 1) == 0) /* 1^y is always 1 */
@@ -289,20 +314,20 @@ mpfr_pow (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd_mode)
     exy *= mpfr_get_d (y, GMP_RNDZ);
     if (exy >= (double) __gmpfr_emax)
       {
-        negative = MPFR_SIGN(x) < 0 && is_odd (y) == 1;
+        negative = MPFR_SIGN(x) < 0 && is_odd (y);
         return mpfr_set_overflow (z, rnd_mode, negative ? -1 : 1);
       }
-  }  
+  }
 
   if (mpfr_integer_p (y))
     {
       mpz_t zi;
       long int zii;
       int exptol;
-    
+
       mpz_init (zi);
       exptol = mpfr_get_z_exp (zi, y);
-        
+
       if (exptol > 0)
         mpz_mul_2exp (zi, zi, exptol);
       else
@@ -314,7 +339,6 @@ mpfr_pow (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd_mode)
       return mpfr_pow_si (z, x, zii, rnd_mode);
     }
 
-  /* now y is not an integer */
   if (MPFR_IS_NEG(x))
     {
       MPFR_SET_NAN(z);