Fix powf (close to -1, large) (bug 18647) [committed]

[Joseph Myers <joseph at codesourcery dot com>]

The flt-32 implementation of powf wrongly uses x-1 instead of |x|-1
when computing log (x) for the case where |x| is close to 1 and y is
large.  This patch fixes the logic accordingly.  Relevant tests
existed for x close to 1, and corresponding tests are added for x
close to -1, as well as for some new variant cases.

Tested for x86_64 and x86.  Committed.

(auto-libm-test-out diffs omitted below.)

2015-08-05  Joseph Myers  <joseph@codesourcery.com>

	[BZ #18647]
	* sysdeps/ieee754/flt-32/e_powf.c (__ieee754_powf): For large y
	and |x| close to 1, use absolute value of x when computing log.
	* math/auto-libm-test-in: Add more tests of pow.
	* math/auto-libm-test-out: Regenerated.

diff --git a/math/auto-libm-test-in b/math/auto-libm-test-in
index f783c5b..c17f6f5 100644
--- a/math/auto-libm-test-in
+++ b/math/auto-libm-test-in
@@ -2361,6 +2361,46 @@ pow 0x0.ffffffffffffffffffffffffffff8p0 -0x1.23456789abcdef0123456789abcdp126
 pow 0x1.0000000000000000000000000001p0 0x1.23456789abcdef0123456789abcdp125
 pow 0x1.0000000000000000000000000001p0 -0x1.23456789abcdef0123456789abcdp125
 
+pow -0x0.ffffffp0 10
+pow -0x0.ffffffp0 100
+pow -0x0.ffffffp0 1000
+pow -0x0.ffffffp0 0x1p24
+pow -0x0.ffffffp0 0x1p30
+pow -0x0.ffffffp0 0x1.234566p30
+pow -0x0.ffffffp0 -10
+pow -0x0.ffffffp0 -100
+pow -0x0.ffffffp0 -1000
+pow -0x0.ffffffp0 -0x1p24
+pow -0x0.ffffffp0 -0x1p30
+pow -0x0.ffffffp0 -0x1.234566p30
+pow -0x1.000002p0 0x1p24
+pow -0x1.000002p0 0x1.234566p29
+pow -0x1.000002p0 -0x1.234566p29
+
+pow -0x0.fffffffffffff8p0 0x1.23456789abcdfp62
+pow -0x0.fffffffffffff8p0 -0x1.23456789abcdfp62
+pow -0x1.0000000000001p0 0x1.23456789abcdfp61
+pow -0x1.0000000000001p0 -0x1.23456789abcdfp61
+
+pow -0x0.ffffffffffffffffp0 0x1.23456789abcdef0ep77
+pow -0x0.ffffffffffffffffp0 -0x1.23456789abcdef0ep77
+pow -0x1.0000000000000002p0 0x1.23456789abcdef0ep76
+pow -0x1.0000000000000002p0 -0x1.23456789abcdef0ep76
+
+pow -0x0.ffffffffffffffffffffffffffff8p0 0x1.23456789abcdef0123456789abcdp126
+pow -0x0.ffffffffffffffffffffffffffff8p0 -0x1.23456789abcdef0123456789abcdp126
+pow -0x1.0000000000000000000000000001p0 0x1.23456789abcdef0123456789abcdp125
+pow -0x1.0000000000000000000000000001p0 -0x1.23456789abcdef0123456789abcdp125
+
+pow 0x1.000002p0 0x1p30
+pow -0x1.000002p0 0x1p30
+pow 0x1.000002p0 max
+pow -0x1.000002p0 max
+pow 0x1.00000ep0 0x1p30
+pow -0x1.00000ep0 0x1p30
+pow 0x1.00000ep0 max
+pow -0x1.00000ep0 max
+
 pow 1e4932 0.75
 pow 1e4928 0.75
 pow 1e4924 0.75
diff --git a/sysdeps/ieee754/flt-32/e_powf.c b/sysdeps/ieee754/flt-32/e_powf.c
index 12c408f..8e8d918 100644
--- a/sysdeps/ieee754/flt-32/e_powf.c
+++ b/sysdeps/ieee754/flt-32/e_powf.c
@@ -131,7 +131,7 @@ __ieee754_powf(float x, float y)
 	    if(ix>0x3f800007) return (hy>0)? huge*huge:tiny*tiny;
 	/* now |1-x| is tiny <= 2**-20, suffice to compute
 	   log(x) by x-x^2/2+x^3/3-x^4/4 */
-	    t = x-1;		/* t has 20 trailing zeros */
+	    t = ax-1;		/* t has 20 trailing zeros */
 	    w = (t*t)*((float)0.5-t*((float)0.333333333333-t*(float)0.25));
 	    u = ivln2_h*t;	/* ivln2_h has 16 sig. bits */
 	    v = t*ivln2_l-w*ivln2;

-- 
Joseph S. Myers
joseph@codesourcery.com