changed back to C99 standard for special values of x^y

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

3 files changed, 26 insertions, 27 deletions

diff --git a/mpfr.texi b/mpfr.texi
index d89b2b7a7..c153d1212 100644
--- a/mpfr.texi
+++ b/mpfr.texi
@@ -930,7 +930,7 @@ rounded in the direction @var{rnd}.
 Set @var{rop} to @m{@var{op1}^{op2}, @var{op1} raised to @var{op2}},
 rounded in the direction @var{rnd}.
 Special values are currently handled as described in the ISO C99 standard
-for the @code{pow} function:
+for the @code{pow} function (note this may change in future versions):
 @itemize @bullet
 @item @code{pow(@pom{}0, @var{y})} returns plus or minus infinity for @var{y} a negative odd integer.
 @item @code{pow(@pom{}0, @var{y})} returns plus infinity for @var{y} negative and not an odd integer.
diff --git a/pow.c b/pow.c
index 9142497da..60311b320 100644
--- a/pow.c
+++ b/pow.c
@@ -140,7 +140,7 @@ is_odd (mpfr_srcptr y)
 
 /* 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:
+   For the special cases, see Section F.9.4.4 of the C99 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.
@@ -226,37 +226,33 @@ 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 (even non-zero) */
+      /* 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 NaN for y non-integer
-                          s*Inf for y positive integer of sign s
-                          s*0   for y negative integer of sign s */
+             -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);
 	  MPFR_CLEAR_FLAGS(z);
           if (negative) /* x = -Inf */
             {
-              int odd = is_odd (y);
-              if (odd == 0) /* non-integer */
-                {
-                  MPFR_SET_NAN(z);
-                  MPFR_RET_NAN;
-                }
-              else 
-                {
-                  if (MPFR_IS_POS(y))
-                    MPFR_SET_INF(z);
-                  else
-                    MPFR_SET_ZERO(z);
-                  if (odd == 1) /* 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 (is_odd (y) == 1) /* y is an odd integer */
+                MPFR_SET_NEG(z);
+              else
+                MPFR_SET_POS(z);
             }
           else /* x = +Inf */
             {
diff --git a/tests/tpow.c b/tests/tpow.c
index 0d2f1a1e1..3112ab747 100644
--- a/tests/tpow.c
+++ b/tests/tpow.c
@@ -270,22 +270,25 @@ special ()
   MPFR_ASSERTN(mpfr_inf_p (z) && MPFR_IS_NEG(z));
   mpfr_set_str_binary (y, "1.000000001E8");
   mpfr_pow (z, x, y, GMP_RNDN);
-  MPFR_ASSERTN(mpfr_nan_p (z));
+  MPFR_ASSERTN(mpfr_inf_p (z) && MPFR_IS_POS(z));
 
   mpfr_set_inf (x, -1);
   mpfr_set_prec (y, 2 * mp_bits_per_limb);
   mpfr_set_ui (y, 1, GMP_RNDN);
   mpfr_mul_2exp (y, y, mp_bits_per_limb - 1, GMP_RNDN);
+  /* y = 2^(mp_bits_per_limb - 1) */
   mpfr_pow (z, x, y, GMP_RNDN);
   MPFR_ASSERTN(mpfr_inf_p (z) && MPFR_IS_POS(z));
   mpfr_nextabove (y);
   mpfr_pow (z, x, y, GMP_RNDN);
-  MPFR_ASSERTN(mpfr_nan_p (z));
+  /* y = 2^(mp_bits_per_limb - 1) + epsilon */
+  MPFR_ASSERTN(mpfr_inf_p (z) && MPFR_IS_POS(z));
   mpfr_nextbelow (y);
   mpfr_div_2exp (y, y, 1, GMP_RNDN);
   mpfr_nextabove (y);
   mpfr_pow (z, x, y, GMP_RNDN);
-  MPFR_ASSERTN(mpfr_nan_p (z));
+  /* y = 2^(mp_bits_per_limb - 2) + epsilon */
+  MPFR_ASSERTN(mpfr_inf_p (z) && MPFR_IS_POS(z));
 
   mpfr_set_si (x, -1, GMP_RNDN);
   mpfr_set_prec (y, 2);
@@ -331,7 +334,7 @@ particular_cases (void)
           /*          NaN +inf -inf  +0   -0   +1   -1   +2   -2  +0.5 -0.5 */
           /*  NaN */ { 0,   0,   0,  128, 128,  0,   0,   0,   0,   0,   0  },
           /* +inf */ { 0,   1,   2,  128, 128,  1,   2,   1,   2,   1,   2  },
-          /* -inf */ { 0,   1,   2,  128, 128, -1,  -2,   1,   2,   0,   0  },
+          /* -inf */ { 0,   1,   2,  128, 128, -1,  -2,   1,   2,   1,   2  },
           /*  +0  */ { 0,   2,   1,  128, 128,  2,   1,   2,   1,   2,   1  },
           /*  -0  */ { 0,   2,   1,  128, 128, -2,  -1,   2,   1,   2,   1  },
           /*  +1  */ {128, 128, 128, 128, 128, 128, 128, 128, 128, 128, 128 },