[git commit master 1/1] libm: fix powf testcase failures

[Denys Vlasenko]

commit: http://git.uclibc.org/uClibc/commit/?id=ae73aafe99fa6bb5e7422f2bdedea39f03ead72c
branch: http://git.uclibc.org/uClibc/commit/?id=refs/heads/master

Fixed failures:
Failure: Test: pow (1, NaN) == 1
Failure: Test: pow (1, inf) == 1
Failure: Test: pow (-1, inf) == 1
Failure: Test: pow (1, -inf) == 1
Failure: Test: pow (-1, -inf) == 1

Signed-off-by: Denys Vlasenko <vda.linux at googlemail.com>
---
 libm/e_pow.c |   63 ++++++++++++++++++++++++++++++++-------------------------
 1 files changed, 35 insertions(+), 28 deletions(-)

diff --git a/libm/e_pow.c b/libm/e_pow.c
index 137f600..3be9003 100644
--- a/libm/e_pow.c
+++ b/libm/e_pow.c
@@ -21,25 +21,26 @@
  *	3. Return x**y = 2**n*exp(y'*log2)
  *
  * Special cases:
- *	1.  (anything) ** 0  is 1
- *	2.  (anything) ** 1  is itself
- *	3.  (anything) ** 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
- *	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
- *	13. -0 ** (-anything except 0, NAN, odd integer)  is +INF
- *	14. -0 ** (odd integer) = -( +0 ** (odd integer) )
- *	15. +INF ** (+anything except 0,NAN) is +INF
- *	16. +INF ** (-anything except 0,NAN) is +0
- *	17. -INF ** (anything)  = -0 ** (-anything)
- *	18. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer)
- *	19. (-anything except 0 and inf) ** (non-integer) is NAN
+ *	1.  +-1 ** anything  is 1.0
+ *	2.  +-1 ** +-INF     is 1.0
+ *	3.  (anything) ** 0  is 1
+ *	4.  (anything) ** 1  is itself
+ *	5.  (anything) ** NAN is NAN
+ *	6.  NAN ** (anything except 0) is NAN
+ *	7.  +-(|x| > 1) **  +INF is +INF
+ *	8.  +-(|x| > 1) **  -INF is +0
+ *	9.  +-(|x| < 1) **  +INF is +0
+ *	10  +-(|x| < 1) **  -INF is +INF
+ *	11. +0 ** (+anything except 0, NAN)               is +0
+ *	12. -0 ** (+anything except 0, NAN, odd integer)  is +0
+ *	13. +0 ** (-anything except 0, NAN)               is +INF
+ *	14. -0 ** (-anything except 0, NAN, odd integer)  is +INF
+ *	15. -0 ** (odd integer) = -( +0 ** (odd integer) )
+ *	16. +INF ** (+anything except 0,NAN) is +INF
+ *	17. +INF ** (-anything except 0,NAN) is +0
+ *	18. -INF ** (anything)  = -0 ** (-anything)
+ *	19. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer)
+ *	20. (-anything except 0 and inf) ** (non-integer) is NAN
  *
  * Accuracy:
  *	pow(x,y) returns x**y nearly rounded. In particular
@@ -99,8 +100,14 @@ double attribute_hidden __ieee754_pow(double x, double y)
 	u_int32_t lx,ly;
 
 	EXTRACT_WORDS(hx,lx,x);
+    /* x==1: 1**y = 1 (even if y is NaN) */
+	if (hx==0x3ff00000 && lx==0) {
+		return x;
+	}
+	ix = hx&0x7fffffff;
+
 	EXTRACT_WORDS(hy,ly,y);
-	ix = hx&0x7fffffff;  iy = hy&0x7fffffff;
+	iy = hy&0x7fffffff;
 
     /* y==zero: x**0 = 1 */
 	if((iy|ly)==0) return one;
@@ -132,13 +139,13 @@ double attribute_hidden __ieee754_pow(double x, double y)
 
     /* special value of y */
 	if(ly==0) {
-	    if (iy==0x7ff00000) {	/* y is +-inf */
-	        if(((ix-0x3ff00000)|lx)==0)
-		    return  y - y;	/* inf**+-1 is NaN */
-	        else if (ix >= 0x3ff00000)/* (|x|>1)**+-inf = inf,0 */
-		    return (hy>=0)? y: zero;
-	        else			/* (|x|<1)**-,+inf = inf,0 */
-		    return (hy<0)?-y: zero;
+	    if (iy==0x7ff00000) {       /* y is +-inf */
+	        if (((ix-0x3ff00000)|lx)==0)
+		    return one;	        /* +-1**+-inf is 1 (yes, weird rule) */
+	        if (ix >= 0x3ff00000)   /* (|x|>1)**+-inf = inf,0 */
+		    return (hy>=0) ? y : zero;
+	        /* (|x|<1)**-,+inf = inf,0 */
+		return (hy<0) ? -y : zero;
 	    }
 	    if(iy==0x3ff00000) {	/* y is  +-1 */
 		if(hy<0) return one/x; else return x;
@@ -146,7 +153,7 @@ double attribute_hidden __ieee754_pow(double x, double y)
 	    if(hy==0x40000000) return x*x; /* y is  2 */
 	    if(hy==0x3fe00000) {	/* y is  0.5 */
 		if(hx>=0)	/* x >= +0 */
-		return __ieee754_sqrt(x);
+		    return __ieee754_sqrt(x);
 	    }
 	}
 
-- 
1.7.1