Fix a couple of corner cases in the implementation of pow(3) to make it

compatible with C99.  Most notably:

  - 1**y == 1, even if y is NaN
  - (-1)**+-Inf == 1

and adjust the cephes testsuite to test for the right thing here.

ok martynas@

3 files changed, 15 insertions, 9 deletions

diff --git a/lib/libm/src/e_pow.c b/lib/libm/src/e_pow.c
index 5fce17c8c7b..2af0e921b3e 100644
--- a/lib/libm/src/e_pow.c
+++ b/lib/libm/src/e_pow.c
@@ -24,13 +24,13 @@
  * Special cases:
  *	1.  (anything) ** 0  is 1
  *	2.  (anything) ** 1  is itself
- *	3.  (anything) ** NAN is NAN
+ *	3.  (anything except 1) ** NAN is NAN
  *	4.  NAN ** (anything except 0) is NAN
  *	5.  +-(|x| > 1) **  +INF is +INF
  *	6.  +-(|x| > 1) **  -INF is +0
  *	7.  +-(|x| < 1) **  +INF is +0
  *	8.  +-(|x| < 1) **  -INF is +INF
- *	9.  +-1         ** +-INF is NAN
+ *	9.  +-1         ** +-INF is 1
  *	10. +0 ** (+anything except 0, NAN)               is +0
  *	11. -0 ** (+anything except 0, NAN, odd integer)  is +0
  *	12. +0 ** (-anything except 0, NAN)               is +INF
@@ -109,6 +109,9 @@ pow(double x, double y)
     /* y==zero: x**0 = 1 */
 	if((iy|ly)==0) return one;
 
+    /* x==1: 1**y = 1, even if y is NaN */
+	if (hx==0x3ff00000 && lx == 0) return one;
+
     /* +-NaN return x+y */
 	if(ix > 0x7ff00000 || ((ix==0x7ff00000)&&(lx!=0)) ||
 	   iy > 0x7ff00000 || ((iy==0x7ff00000)&&(ly!=0)))
@@ -138,7 +141,7 @@ pow(double x, double y)
 	if(ly==0) {
 	    if (iy==0x7ff00000) {	/* y is +-inf */
 	        if(((ix-0x3ff00000)|lx)==0)
-		    return  y - y;	/* inf**+-1 is NaN */
+		    return  one;	/* (-1)**+-inf is 1 */
 	        else if (ix >= 0x3ff00000)/* (|x|>1)**+-inf = inf,0 */
 		    return (hy>=0)? y: zero;
 	        else			/* (|x|<1)**-,+inf = inf,0 */
diff --git a/lib/libm/src/e_powf.c b/lib/libm/src/e_powf.c
index d6e52d8b99b..3c4269d7146 100644
--- a/lib/libm/src/e_powf.c
+++ b/lib/libm/src/e_powf.c
@@ -64,6 +64,9 @@ powf(float x, float y)
     /* y==zero: x**0 = 1 */
 	if(iy==0) return one; 	
 
+    /* x==1: 1**y = 1, even if y is NaN */
+	if (hx==0x3f800000) return one;
+
     /* +-NaN return x+y */
 	if(ix > 0x7f800000 ||
 	   iy > 0x7f800000)
@@ -87,7 +90,7 @@ powf(float x, float y)
     /* special value of y */
 	if (iy==0x7f800000) {	/* y is +-inf */
 	    if (ix==0x3f800000)
-	        return  y - y;	/* inf**+-1 is NaN */
+	        return  one;	/* (-1)**+-inf is NaN */
 	    else if (ix > 0x3f800000)/* (|x|>1)**+-inf = inf,0 */
 	        return (hy>=0)? y: zero;
 	    else			/* (|x|<1)**-,+inf = inf,0 */
diff --git a/lib/libm/src/ld80/e_powl.c b/lib/libm/src/ld80/e_powl.c
index d05de530531..f72dd8d68b8 100644
--- a/lib/libm/src/ld80/e_powl.c
+++ b/lib/libm/src/ld80/e_powl.c
@@ -1,4 +1,4 @@
-/*	$OpenBSD: e_powl.c,v 1.3 2011/07/20 21:02:51 martynas Exp $	*/
+/*	$OpenBSD: e_powl.c,v 1.4 2013/08/02 20:23:28 kettenis Exp $	*/
 
 /*
  * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
@@ -211,6 +211,9 @@ long e;
 if( y == 0.0L )
 	return( 1.0L );
 
+if( x == 1.0L )
+	return( 1.0L );
+
 if( isnan(x) )
 	return( x );
 if( isnan(y) )
@@ -219,10 +222,7 @@ if( isnan(y) )
 if( y == 1.0L )
 	return( x );
 
-if( !isfinite(y) && (x == -1.0L || x == 1.0L) )
-	return y - y;	/* +-1**inf is NaN */
-
-if( x == 1.0L )
+if( !isfinite(y) && x == -1.0L )
 	return( 1.0L );
 
 if( y >= LDBL_MAX )