This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
Fix pow of negative numbers to integer exponents (bugs 369, 2678,3866)
- From: "Joseph S. Myers" <joseph at codesourcery dot com>
- To: libc-alpha at sourceware dot org
- Date: Wed, 28 Mar 2012 13:41:21 +0000 (UTC)
- Subject: Fix pow of negative numbers to integer exponents (bugs 369, 2678,3866)
This patch fixes various problems with pow of negative numbers to
integer exponents, as reported in bugs 369, 2678 and 3866 (which this
patch completes fixing). The assembly implementations for x86 and
x86_64 needed fixing to handle these cases (as noted in bug 369, the
patch there was incorrect); powl was the most complicated case because
of the need to determine parity rather than presume that "large"
integer exponents must be even, though pow will need such
complications as well later as part of fixing bug 706.
__kernel_standard needed to handle getting the sign of zero right for
underflow. __kernel_standard_l needed to handle pow overflow and
underflow directly itself, rather than passing to __kernel_standard,
because the sign of the result depends on the parity of an integer
exponent that may not round exactly to double.
Tested x86 and x86_64.
2012-03-28 Joseph Myers <joseph@codesourcery.com>
[BZ #369]
[BZ #2678]
[BZ #3866]
* sysdeps/i386/fpu/e_pow.S (__ieee754_pow): Take absolute value of
x for large integer exponent.
* sysdeps/i386/fpu/e_powf.S (__ieee754_powf): Likewise.
* sysdeps/i386/fpu/e_powl.S (__ieee754_powl): Likewise. Adjust
sign of result as needed afterwards.
* sysdeps/x86_64/fpu/e_powl.S (__ieee754_powl): Likewise.
* sysdeps/ieee754/k_standard.c (__kernel_standard): Handle sign of
result for underflowing pow the same as for overflow.
(__kernel_standard_l): Handle powl overflow and underflow here
rather than calling __kernel_standard.
* math/libm-test.inc (pow_test): Add more tests.
diff --git a/math/libm-test.inc b/math/libm-test.inc
index 68f6ef2..32bce45 100644
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -5818,6 +5818,259 @@ pow_test (void)
TEST_ff_f (pow, -7.49321e+133, -9.80818e+16, 0);
#endif
+ TEST_ff_f (pow, -1.0, -0xffffff, -1.0);
+ TEST_ff_f (pow, -1.0, -0x1fffffe, 1.0);
+#ifndef TEST_FLOAT
+ TEST_ff_f (pow, -1.0, -0x1.fffffffffffffp+52L, -1.0);
+ TEST_ff_f (pow, -1.0, -0x1.fffffffffffffp+53L, 1.0);
+#endif
+#ifdef TEST_LDOUBLE
+# if LDBL_MANT_DIG >= 64
+ TEST_ff_f (pow, -1.0, -0x1.fffffffffffffffep+63L, -1.0);
+ TEST_ff_f (pow, -1.0, -0x1.fffffffffffffffep+64L, 1.0);
+# endif
+# if LDBL_MANT_DIG >= 106
+ TEST_ff_f (pow, -1.0, -0x1.ffffffffffffffffffffffffff8p+105L, -1.0);
+ TEST_ff_f (pow, -1.0, -0x1.ffffffffffffffffffffffffff8p+106L, 1.0);
+# endif
+# if LDBL_MANT_DIG >= 113
+ TEST_ff_f (pow, -1.0, -0x1.ffffffffffffffffffffffffffffp+112L, -1.0);
+ TEST_ff_f (pow, -1.0, -0x1.ffffffffffffffffffffffffffffp+113L, 1.0);
+# endif
+#endif
+ TEST_ff_f (pow, -1.0, -max_value, 1.0);
+
+ TEST_ff_f (pow, -1.0, 0xffffff, -1.0);
+ TEST_ff_f (pow, -1.0, 0x1fffffe, 1.0);
+#ifndef TEST_FLOAT
+ TEST_ff_f (pow, -1.0, 0x1.fffffffffffffp+52L, -1.0);
+ TEST_ff_f (pow, -1.0, 0x1.fffffffffffffp+53L, 1.0);
+#endif
+#ifdef TEST_LDOUBLE
+# if LDBL_MANT_DIG >= 64
+ TEST_ff_f (pow, -1.0, 0x1.fffffffffffffffep+63L, -1.0);
+ TEST_ff_f (pow, -1.0, 0x1.fffffffffffffffep+64L, 1.0);
+# endif
+# if LDBL_MANT_DIG >= 106
+ TEST_ff_f (pow, -1.0, 0x1.ffffffffffffffffffffffffff8p+105L, -1.0);
+ TEST_ff_f (pow, -1.0, 0x1.ffffffffffffffffffffffffff8p+106L, 1.0);
+# endif
+# if LDBL_MANT_DIG >= 113
+ TEST_ff_f (pow, -1.0, 0x1.ffffffffffffffffffffffffffffp+112L, -1.0);
+ TEST_ff_f (pow, -1.0, 0x1.ffffffffffffffffffffffffffffp+113L, 1.0);
+# endif
+#endif
+ TEST_ff_f (pow, -1.0, max_value, 1.0);
+
+ TEST_ff_f (pow, -2.0, 126, 0x1p126);
+ TEST_ff_f (pow, -2.0, 127, -0x1p127);
+ TEST_ff_f (pow, -2.0, -126, 0x1p-126);
+ TEST_ff_f (pow, -2.0, -127, -0x1p-127);
+
+ TEST_ff_f (pow, -2.0, -0xffffff, minus_zero);
+ TEST_ff_f (pow, -2.0, -0x1fffffe, plus_zero);
+#ifndef TEST_FLOAT
+ TEST_ff_f (pow, -2.0, -0x1.fffffffffffffp+52L, minus_zero);
+ TEST_ff_f (pow, -2.0, -0x1.fffffffffffffp+53L, plus_zero);
+#endif
+#ifdef TEST_LDOUBLE
+# if LDBL_MANT_DIG >= 64
+ TEST_ff_f (pow, -2.0, -0x1.fffffffffffffffep+63L, minus_zero);
+ TEST_ff_f (pow, -2.0, -0x1.fffffffffffffffep+64L, plus_zero);
+# endif
+# if LDBL_MANT_DIG >= 106
+ TEST_ff_f (pow, -2.0, -0x1.ffffffffffffffffffffffffff8p+105L, minus_zero);
+ TEST_ff_f (pow, -2.0, -0x1.ffffffffffffffffffffffffff8p+106L, plus_zero);
+# endif
+# if LDBL_MANT_DIG >= 113
+ TEST_ff_f (pow, -2.0, -0x1.ffffffffffffffffffffffffffffp+112L, minus_zero);
+ TEST_ff_f (pow, -2.0, -0x1.ffffffffffffffffffffffffffffp+113L, plus_zero);
+# endif
+#endif
+ TEST_ff_f (pow, -2.0, -max_value, plus_zero);
+
+ TEST_ff_f (pow, -2.0, 0xffffff, minus_infty, OVERFLOW_EXCEPTION);
+ TEST_ff_f (pow, -2.0, 0x1fffffe, plus_infty, OVERFLOW_EXCEPTION);
+#ifndef TEST_FLOAT
+ TEST_ff_f (pow, -2.0, 0x1.fffffffffffffp+52L, minus_infty, OVERFLOW_EXCEPTION);
+ TEST_ff_f (pow, -2.0, 0x1.fffffffffffffp+53L, plus_infty, OVERFLOW_EXCEPTION);
+#endif
+#ifdef TEST_LDOUBLE
+# if LDBL_MANT_DIG >= 64
+ TEST_ff_f (pow, -2.0, 0x1.fffffffffffffffep+63L, minus_infty, OVERFLOW_EXCEPTION);
+ TEST_ff_f (pow, -2.0, 0x1.fffffffffffffffep+64L, plus_infty, OVERFLOW_EXCEPTION);
+# endif
+# if LDBL_MANT_DIG >= 106
+ TEST_ff_f (pow, -2.0, 0x1.ffffffffffffffffffffffffff8p+105L, minus_infty, OVERFLOW_EXCEPTION);
+ TEST_ff_f (pow, -2.0, 0x1.ffffffffffffffffffffffffff8p+106L, plus_infty, OVERFLOW_EXCEPTION);
+# endif
+# if LDBL_MANT_DIG >= 113
+ TEST_ff_f (pow, -2.0, 0x1.ffffffffffffffffffffffffffffp+112L, minus_infty, OVERFLOW_EXCEPTION);
+ TEST_ff_f (pow, -2.0, 0x1.ffffffffffffffffffffffffffffp+113L, plus_infty, OVERFLOW_EXCEPTION);
+# endif
+#endif
+ /* Bug 13873: OVERFLOW exception may be missing. */
+ TEST_ff_f (pow, -2.0, max_value, plus_infty, OVERFLOW_EXCEPTION_OK);
+
+ TEST_ff_f (pow, -max_value, 0.5, nan_value, INVALID_EXCEPTION);
+ TEST_ff_f (pow, -max_value, 1.5, nan_value, INVALID_EXCEPTION);
+ TEST_ff_f (pow, -max_value, 1000.5, nan_value, INVALID_EXCEPTION);
+ TEST_ff_f (pow, -max_value, -2, plus_zero);
+ TEST_ff_f (pow, -max_value, -3, minus_zero);
+ TEST_ff_f (pow, -max_value, 2, plus_infty, OVERFLOW_EXCEPTION);
+ TEST_ff_f (pow, -max_value, 3, minus_infty, OVERFLOW_EXCEPTION);
+
+ TEST_ff_f (pow, -max_value, -0xffffff, minus_zero);
+ TEST_ff_f (pow, -max_value, -0x1fffffe, plus_zero);
+#ifndef TEST_FLOAT
+ TEST_ff_f (pow, -max_value, -0x1.fffffffffffffp+52L, minus_zero);
+ TEST_ff_f (pow, -max_value, -0x1.fffffffffffffp+53L, plus_zero);
+#endif
+#ifdef TEST_LDOUBLE
+# if LDBL_MANT_DIG >= 64
+ TEST_ff_f (pow, -max_value, -0x1.fffffffffffffffep+63L, minus_zero);
+ TEST_ff_f (pow, -max_value, -0x1.fffffffffffffffep+64L, plus_zero);
+# endif
+# if LDBL_MANT_DIG >= 106
+ TEST_ff_f (pow, -max_value, -0x1.ffffffffffffffffffffffffff8p+105L, minus_zero);
+ TEST_ff_f (pow, -max_value, -0x1.ffffffffffffffffffffffffff8p+106L, plus_zero);
+# endif
+# if LDBL_MANT_DIG >= 113
+ TEST_ff_f (pow, -max_value, -0x1.ffffffffffffffffffffffffffffp+112L, minus_zero);
+ TEST_ff_f (pow, -max_value, -0x1.ffffffffffffffffffffffffffffp+113L, plus_zero);
+# endif
+#endif
+ /* Bug 13872: spurious OVERFLOW exception may be present. */
+ TEST_ff_f (pow, -max_value, -max_value, plus_zero, OVERFLOW_EXCEPTION_OK);
+
+ TEST_ff_f (pow, -max_value, 0xffffff, minus_infty, OVERFLOW_EXCEPTION);
+ TEST_ff_f (pow, -max_value, 0x1fffffe, plus_infty, OVERFLOW_EXCEPTION);
+#ifndef TEST_FLOAT
+ TEST_ff_f (pow, -max_value, 0x1.fffffffffffffp+52L, minus_infty, OVERFLOW_EXCEPTION);
+ TEST_ff_f (pow, -max_value, 0x1.fffffffffffffp+53L, plus_infty, OVERFLOW_EXCEPTION);
+#endif
+#ifdef TEST_LDOUBLE
+# if LDBL_MANT_DIG >= 64
+ TEST_ff_f (pow, -max_value, 0x1.fffffffffffffffep+63L, minus_infty, OVERFLOW_EXCEPTION);
+ TEST_ff_f (pow, -max_value, 0x1.fffffffffffffffep+64L, plus_infty, OVERFLOW_EXCEPTION);
+# endif
+# if LDBL_MANT_DIG >= 106
+ TEST_ff_f (pow, -max_value, 0x1.ffffffffffffffffffffffffff8p+105L, minus_infty, OVERFLOW_EXCEPTION);
+ TEST_ff_f (pow, -max_value, 0x1.ffffffffffffffffffffffffff8p+106L, plus_infty, OVERFLOW_EXCEPTION);
+# endif
+# if LDBL_MANT_DIG >= 113
+ TEST_ff_f (pow, -max_value, 0x1.ffffffffffffffffffffffffffffp+112L, minus_infty, OVERFLOW_EXCEPTION);
+ TEST_ff_f (pow, -max_value, 0x1.ffffffffffffffffffffffffffffp+113L, plus_infty, OVERFLOW_EXCEPTION);
+# endif
+#endif
+ /* Bug 13873: OVERFLOW exception may be missing. */
+ TEST_ff_f (pow, -max_value, max_value, plus_infty, OVERFLOW_EXCEPTION_OK);
+
+ TEST_ff_f (pow, -0.5, 126, 0x1p-126);
+ TEST_ff_f (pow, -0.5, 127, -0x1p-127);
+ TEST_ff_f (pow, -0.5, -126, 0x1p126);
+ TEST_ff_f (pow, -0.5, -127, -0x1p127);
+
+ TEST_ff_f (pow, -0.5, -0xffffff, minus_infty, OVERFLOW_EXCEPTION);
+ TEST_ff_f (pow, -0.5, -0x1fffffe, plus_infty, OVERFLOW_EXCEPTION);
+#ifndef TEST_FLOAT
+ TEST_ff_f (pow, -0.5, -0x1.fffffffffffffp+52L, minus_infty, OVERFLOW_EXCEPTION);
+ TEST_ff_f (pow, -0.5, -0x1.fffffffffffffp+53L, plus_infty, OVERFLOW_EXCEPTION);
+#endif
+#ifdef TEST_LDOUBLE
+# if LDBL_MANT_DIG >= 64
+ TEST_ff_f (pow, -0.5, -0x1.fffffffffffffffep+63L, minus_infty, OVERFLOW_EXCEPTION);
+ TEST_ff_f (pow, -0.5, -0x1.fffffffffffffffep+64L, plus_infty, OVERFLOW_EXCEPTION);
+# endif
+# if LDBL_MANT_DIG >= 106
+ TEST_ff_f (pow, -0.5, -0x1.ffffffffffffffffffffffffff8p+105L, minus_infty, OVERFLOW_EXCEPTION);
+ TEST_ff_f (pow, -0.5, -0x1.ffffffffffffffffffffffffff8p+106L, plus_infty, OVERFLOW_EXCEPTION);
+# endif
+# if LDBL_MANT_DIG >= 113
+ TEST_ff_f (pow, -0.5, -0x1.ffffffffffffffffffffffffffffp+112L, minus_infty, OVERFLOW_EXCEPTION);
+ TEST_ff_f (pow, -0.5, -0x1.ffffffffffffffffffffffffffffp+113L, plus_infty, OVERFLOW_EXCEPTION);
+# endif
+#endif
+ /* Bug 13873: OVERFLOW exception may be missing. */
+ TEST_ff_f (pow, -0.5, -max_value, plus_infty, OVERFLOW_EXCEPTION_OK);
+
+ TEST_ff_f (pow, -0.5, 0xffffff, minus_zero);
+ TEST_ff_f (pow, -0.5, 0x1fffffe, plus_zero);
+#ifndef TEST_FLOAT
+ TEST_ff_f (pow, -0.5, 0x1.fffffffffffffp+52L, minus_zero);
+ TEST_ff_f (pow, -0.5, 0x1.fffffffffffffp+53L, plus_zero);
+#endif
+#ifdef TEST_LDOUBLE
+# if LDBL_MANT_DIG >= 64
+ TEST_ff_f (pow, -0.5, 0x1.fffffffffffffffep+63L, minus_zero);
+ TEST_ff_f (pow, -0.5, 0x1.fffffffffffffffep+64L, plus_zero);
+# endif
+# if LDBL_MANT_DIG >= 106
+ TEST_ff_f (pow, -0.5, 0x1.ffffffffffffffffffffffffff8p+105L, minus_zero);
+ TEST_ff_f (pow, -0.5, 0x1.ffffffffffffffffffffffffff8p+106L, plus_zero);
+# endif
+# if LDBL_MANT_DIG >= 113
+ TEST_ff_f (pow, -0.5, 0x1.ffffffffffffffffffffffffffffp+112L, minus_zero);
+ TEST_ff_f (pow, -0.5, 0x1.ffffffffffffffffffffffffffffp+113L, plus_zero);
+# endif
+#endif
+ TEST_ff_f (pow, -0.5, max_value, plus_zero);
+
+ TEST_ff_f (pow, -min_value, 0.5, nan_value, INVALID_EXCEPTION);
+ TEST_ff_f (pow, -min_value, 1.5, nan_value, INVALID_EXCEPTION);
+ TEST_ff_f (pow, -min_value, 1000.5, nan_value, INVALID_EXCEPTION);
+ TEST_ff_f (pow, -min_value, -2, plus_infty, OVERFLOW_EXCEPTION);
+ TEST_ff_f (pow, -min_value, -3, minus_infty, OVERFLOW_EXCEPTION);
+ TEST_ff_f (pow, -min_value, 1, -min_value);
+ TEST_ff_f (pow, -min_value, 2, plus_zero);
+ TEST_ff_f (pow, -min_value, 3, minus_zero);
+
+ TEST_ff_f (pow, -min_value, -0xffffff, minus_infty, OVERFLOW_EXCEPTION);
+ TEST_ff_f (pow, -min_value, -0x1fffffe, plus_infty, OVERFLOW_EXCEPTION);
+#ifndef TEST_FLOAT
+ TEST_ff_f (pow, -min_value, -0x1.fffffffffffffp+52L, minus_infty, OVERFLOW_EXCEPTION);
+ TEST_ff_f (pow, -min_value, -0x1.fffffffffffffp+53L, plus_infty, OVERFLOW_EXCEPTION);
+#endif
+#ifdef TEST_LDOUBLE
+# if LDBL_MANT_DIG >= 64
+ TEST_ff_f (pow, -min_value, -0x1.fffffffffffffffep+63L, minus_infty, OVERFLOW_EXCEPTION);
+ TEST_ff_f (pow, -min_value, -0x1.fffffffffffffffep+64L, plus_infty, OVERFLOW_EXCEPTION);
+# endif
+# if LDBL_MANT_DIG >= 106
+ TEST_ff_f (pow, -min_value, -0x1.ffffffffffffffffffffffffff8p+105L, minus_infty, OVERFLOW_EXCEPTION);
+ TEST_ff_f (pow, -min_value, -0x1.ffffffffffffffffffffffffff8p+106L, plus_infty, OVERFLOW_EXCEPTION);
+# endif
+# if LDBL_MANT_DIG >= 113
+ TEST_ff_f (pow, -min_value, -0x1.ffffffffffffffffffffffffffffp+112L, minus_infty, OVERFLOW_EXCEPTION);
+ TEST_ff_f (pow, -min_value, -0x1.ffffffffffffffffffffffffffffp+113L, plus_infty, OVERFLOW_EXCEPTION);
+# endif
+#endif
+ /* Bug 13873: OVERFLOW exception may be missing. */
+ TEST_ff_f (pow, -min_value, -max_value, plus_infty, OVERFLOW_EXCEPTION_OK);
+
+ TEST_ff_f (pow, -min_value, 0xffffff, minus_zero);
+ TEST_ff_f (pow, -min_value, 0x1fffffe, plus_zero);
+#ifndef TEST_FLOAT
+ TEST_ff_f (pow, -min_value, 0x1.fffffffffffffp+52L, minus_zero);
+ TEST_ff_f (pow, -min_value, 0x1.fffffffffffffp+53L, plus_zero);
+#endif
+#ifdef TEST_LDOUBLE
+# if LDBL_MANT_DIG >= 64
+ TEST_ff_f (pow, -min_value, 0x1.fffffffffffffffep+63L, minus_zero);
+ TEST_ff_f (pow, -min_value, 0x1.fffffffffffffffep+64L, plus_zero);
+# endif
+# if LDBL_MANT_DIG >= 106
+ TEST_ff_f (pow, -min_value, 0x1.ffffffffffffffffffffffffff8p+105L, minus_zero);
+ TEST_ff_f (pow, -min_value, 0x1.ffffffffffffffffffffffffff8p+106L, plus_zero);
+# endif
+# if LDBL_MANT_DIG >= 113
+ TEST_ff_f (pow, -min_value, 0x1.ffffffffffffffffffffffffffffp+112L, minus_zero);
+ TEST_ff_f (pow, -min_value, 0x1.ffffffffffffffffffffffffffffp+113L, plus_zero);
+# endif
+#endif
+ /* Bug 13872: spurious OVERFLOW exception may be present. */
+ TEST_ff_f (pow, -min_value, max_value, plus_zero, OVERFLOW_EXCEPTION_OK);
+
END (pow);
}
diff --git a/sysdeps/i386/fpu/e_pow.S b/sysdeps/i386/fpu/e_pow.S
index 1abedf6..b61a946 100644
--- a/sysdeps/i386/fpu/e_pow.S
+++ b/sysdeps/i386/fpu/e_pow.S
@@ -114,7 +114,7 @@ ENTRY(__ieee754_pow)
fucomp %st(1) // y : x
fnstsw
sahf
- jne 2f
+ jne 3f
/* OK, we have an integer value for y. */
popl %eax
@@ -157,7 +157,12 @@ ENTRY(__ieee754_pow)
cfi_adjust_cfa_offset (8)
.align ALIGNARG(4)
-2: /* y is a real number. */
+2: /* y is a large integer (so even). */
+ fxch // x : y
+ fabs // |x| : y
+ fxch // y : x
+ .align ALIGNARG(4)
+3: /* y is a real number. */
fxch // x : y
fldl MO(one) // 1.0 : x : y
fldl MO(limit) // 0.29 : 1.0 : x : y
diff --git a/sysdeps/i386/fpu/e_powf.S b/sysdeps/i386/fpu/e_powf.S
index aa58ed2..529a96f 100644
--- a/sysdeps/i386/fpu/e_powf.S
+++ b/sysdeps/i386/fpu/e_powf.S
@@ -114,7 +114,7 @@ ENTRY(__ieee754_powf)
fucomp %st(1) // y : x
fnstsw
sahf
- jne 2f
+ jne 3f
/* OK, we have an integer value for y. */
popl %edx
@@ -151,7 +151,12 @@ ENTRY(__ieee754_powf)
cfi_adjust_cfa_offset (4)
.align ALIGNARG(4)
-2: /* y is a real number. */
+2: /* y is a large integer (so even). */
+ fxch // x : y
+ fabs // |x| : y
+ fxch // y : x
+ .align ALIGNARG(4)
+3: /* y is a real number. */
fxch // x : y
fldl MO(one) // 1.0 : x : y
fldl MO(limit) // 0.29 : 1.0 : x : y
diff --git a/sysdeps/i386/fpu/e_powl.S b/sysdeps/i386/fpu/e_powl.S
index c0aa194c..0e7c05b 100644
--- a/sysdeps/i386/fpu/e_powl.S
+++ b/sysdeps/i386/fpu/e_powl.S
@@ -117,7 +117,7 @@ ENTRY(__ieee754_powl)
fucomp %st(1) // y : x
fnstsw
sahf
- jne 2f
+ jne 3f
/* OK, we have an integer value for y. */
popl %eax
@@ -160,7 +160,14 @@ ENTRY(__ieee754_powl)
cfi_adjust_cfa_offset (8)
.align ALIGNARG(4)
-2: /* y is a real number. */
+2: // y is a large integer (absolute value at least 1L<<63), but
+ // may be odd unless at least 1L<<64. So it may be necessary
+ // to adjust the sign of a negative result afterwards.
+ fxch // x : y
+ fabs // |x| : y
+ fxch // y : |x|
+ .align ALIGNARG(4)
+3: /* y is a real number. */
fxch // x : y
fldl MO(one) // 1.0 : x : y
fldl MO(limit) // 0.29 : 1.0 : x : y
@@ -190,18 +197,51 @@ ENTRY(__ieee754_powl)
f2xm1 // 2^fract(y*log2(x))-1 : int(y*log2(x))
faddl MO(one) // 2^fract(y*log2(x)) : int(y*log2(x))
fscale // 2^fract(y*log2(x))*2^int(y*log2(x)) : int(y*log2(x))
- addl $8, %esp
- cfi_adjust_cfa_offset (-8)
fstp %st(1) // 2^fract(y*log2(x))*2^int(y*log2(x))
- ret
+ jmp 29f
- cfi_adjust_cfa_offset (8)
28: fstp %st(1) // y*log2(x)
fldl MO(one) // 1 : y*log2(x)
fscale // 2^(y*log2(x)) : y*log2(x)
- addl $8, %esp
- cfi_adjust_cfa_offset (-8)
fstp %st(1) // 2^(y*log2(x))
+29: testb $2, %dh
+ jz 292f
+ // x is negative. If y is an odd integer, negate the result.
+ fldt 24(%esp) // y : abs(result)
+ fld %st // y : y : abs(result)
+ fabs // |y| : y : abs(result)
+ fcompl MO(p64) // y : abs(result)
+ fnstsw
+ sahf
+ jnc 291f
+ fldl MO(p63) // p63 : y : abs(result)
+ fxch // y : p63 : abs(result)
+ fprem // y%p63 : p63 : abs(result)
+ fstp %st(1) // y%p63 : abs(result)
+
+ // We must find out whether y is an odd integer.
+ fld %st // y : y : abs(result)
+ fistpll (%esp) // y : abs(result)
+ fildll (%esp) // int(y) : y : abs(result)
+ fucompp // abs(result)
+ fnstsw
+ sahf
+ jne 292f
+
+ // OK, the value is an integer, but is it odd?
+ popl %eax
+ cfi_adjust_cfa_offset (-4)
+ popl %edx
+ cfi_adjust_cfa_offset (-4)
+ andb $1, %al
+ jz 290f // jump if not odd
+ // It's an odd integer.
+ fchs
+290: ret
+ cfi_adjust_cfa_offset (8)
+291: fstp %st(0) // abs(result)
+292: addl $8, %esp
+ cfi_adjust_cfa_offset (-8)
ret
// pow(x,±0) = 1
diff --git a/sysdeps/ieee754/k_standard.c b/sysdeps/ieee754/k_standard.c
index c3326d9..4e65bb1 100644
--- a/sysdeps/ieee754/k_standard.c
+++ b/sysdeps/ieee754/k_standard.c
@@ -508,6 +508,9 @@ __kernel_standard(double x, double y, int type)
exc.type = UNDERFLOW;
exc.name = type < 100 ? "pow" : (type < 200 ? "powf" : "powl");
exc.retval = zero;
+ y *= 0.5;
+ if (x < zero && __rint (y) != y)
+ exc.retval = -zero;
if (_LIB_VERSION == _POSIX_)
__set_errno (ERANGE);
else if (!matherr(&exc)) {
@@ -1004,6 +1007,8 @@ long double
__kernel_standard_l (long double x, long double y, int type)
{
double dx, dy;
+ struct exception exc;
+
if (isfinite (x))
{
long double ax = fabsl (x);
@@ -1028,5 +1033,52 @@ __kernel_standard_l (long double x, long double y, int type)
}
else
dy = y;
- return __kernel_standard (dx, dy, type);
+
+ switch (type)
+ {
+ case 221:
+ /* powl (x, y) overflow. */
+ exc.arg1 = dx;
+ exc.arg2 = dy;
+ exc.type = OVERFLOW;
+ exc.name = "powl";
+ if (_LIB_VERSION == _SVID_)
+ {
+ exc.retval = HUGE;
+ y *= 0.5;
+ if (x < zero && __rintl (y) != y)
+ exc.retval = -HUGE;
+ }
+ else
+ {
+ exc.retval = HUGE_VAL;
+ y *= 0.5;
+ if (x < zero && __rintl (y) != y)
+ exc.retval = -HUGE_VAL;
+ }
+ if (_LIB_VERSION == _POSIX_)
+ __set_errno (ERANGE);
+ else if (!matherr (&exc))
+ __set_errno (ERANGE);
+ return exc.retval;
+
+ case 222:
+ /* powl (x, y) underflow. */
+ exc.arg1 = dx;
+ exc.arg2 = dy;
+ exc.type = UNDERFLOW;
+ exc.name = "powl";
+ exc.retval = zero;
+ y *= 0.5;
+ if (x < zero && __rintl (y) != y)
+ exc.retval = -zero;
+ if (_LIB_VERSION == _POSIX_)
+ __set_errno (ERANGE);
+ else if (!matherr (&exc))
+ __set_errno (ERANGE);
+ return exc.retval;
+
+ default:
+ return __kernel_standard (dx, dy, type);
+ }
}
diff --git a/sysdeps/x86_64/fpu/e_powl.S b/sysdeps/x86_64/fpu/e_powl.S
index 562791d..0626ce4 100644
--- a/sysdeps/x86_64/fpu/e_powl.S
+++ b/sysdeps/x86_64/fpu/e_powl.S
@@ -107,7 +107,7 @@ ENTRY(__ieee754_powl)
fistpll -8(%rsp) // y : x
fildll -8(%rsp) // int(y) : y : x
fucomip %st(1),%st // y : x
- jne 2f
+ jne 3f
/* OK, we have an integer value for y. */
mov -8(%rsp),%eax
@@ -145,7 +145,14 @@ ENTRY(__ieee754_powl)
ret
.align ALIGNARG(4)
-2: /* y is a real number. */
+2: // y is a large integer (absolute value at least 1L<<63), but
+ // may be odd unless at least 1L<<64. So it may be necessary
+ // to adjust the sign of a negative result afterwards.
+ fxch // x : y
+ fabs // |x| : y
+ fxch // y : |x|
+ .align ALIGNARG(4)
+3: /* y is a real number. */
fxch // x : y
fldl MO(one) // 1.0 : x : y
fldl MO(limit) // 0.29 : 1.0 : x : y
@@ -176,13 +183,45 @@ ENTRY(__ieee754_powl)
faddl MO(one) // 2^fract(y*log2(x)) : int(y*log2(x))
fscale // 2^fract(y*log2(x))*2^int(y*log2(x)) : int(y*log2(x))
fstp %st(1) // 2^fract(y*log2(x))*2^int(y*log2(x))
- ret
+ jmp 29f
28: fstp %st(1) // y*log2(x)
fldl MO(one) // 1 : y*log2(x)
fscale // 2^(y*log2(x)) : y*log2(x)
fstp %st(1) // 2^(y*log2(x))
- ret
+29: testb $2, %dh
+ jz 292f
+ // x is negative. If y is an odd integer, negate the result.
+ fldt 24(%rsp) // y : abs(result)
+ fldl MO(p64) // 1L<<64 : y : abs(result)
+ fld %st(1) // y : 1L<<64 : y : abs(result)
+ fabs // |y| : 1L<<64 : y : abs(result)
+ fcomip %st(1), %st // 1L<<64 : y : abs(result)
+ fstp %st(0) // y : abs(result)
+ jnc 291f
+ fldl MO(p63) // p63 : y : abs(result)
+ fxch // y : p63 : abs(result)
+ fprem // y%p63 : p63 : abs(result)
+ fstp %st(1) // y%p63 : abs(result)
+
+ // We must find out whether y is an odd integer.
+ fld %st // y : y : abs(result)
+ fistpll -8(%rsp) // y : abs(result)
+ fildll -8(%rsp) // int(y) : y : abs(result)
+ fucomip %st(1),%st // y : abs(result)
+ ffreep %st // abs(result)
+ jne 292f
+
+ // OK, the value is an integer, but is it odd?
+ mov -8(%rsp), %eax
+ mov -4(%rsp), %edx
+ andb $1, %al
+ jz 290f // jump if not odd
+ // It's an odd integer.
+ fchs
+290: ret
+291: fstp %st(0) // abs(result)
+292: ret
// pow(x,±0) = 1
.align ALIGNARG(4)
--
Joseph S. Myers
joseph@codesourcery.com