From 2460d3aa21f04cdf28497683bd3e29183189f779 Mon Sep 17 00:00:00 2001 From: Joseph Myers Date: Wed, 21 Mar 2012 12:14:57 +0000 Subject: [PATCH] Fix pow of zero and infinity to large powers. --- ChangeLog | 21 +++++ math/libm-test.inc | 207 +++++++++++++++++++++++++++++++++++++++++ sysdeps/i386/fpu/e_pow.S | 47 ++++++---- sysdeps/i386/fpu/e_powf.S | 47 ++++++---- sysdeps/i386/fpu/e_powl.S | 41 ++++++++ sysdeps/ieee754/dbl-64/e_pow.c | 2 +- sysdeps/x86_64/fpu/e_powl.S | 41 ++++++++ 7 files changed, 373 insertions(+), 33 deletions(-) diff --git a/ChangeLog b/ChangeLog index edac2dc..d548878 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,24 @@ +2012-03-20 Joseph Myers + + [BZ #3866] + * sysdeps/i386/fpu/e_pow.S (__ieee754_pow): Test for y outside the + range of signed 64-bit integers before using fistpll. Remove + checks for whether integers fit in mantissa bits. + * sysdeps/i386/fpu/e_powf.S (__ieee754_powf): Test for y outside + the range of signed 32-bit integers before using fistpl. Remove + checks for whether integers fit in mantissa bits. + * sysdeps/i386/fpu/e_powl.S (p64): New object. + (__ieee754_powl): Test for y outside the range of signed 64-bit + integers before using fistpll. Reduce 64-bit values to 63-bit + ones as needed. + * sysdeps/ieee754/dbl-64/e_pow.c (__ieee754_pow): Ensure + divide-by-zero is raised for zero to large negative powers. + * sysdeps/x86_64/fpu/e_powl.S (p64): New object. + (__ieee754_powl): Test for y outside the range of signed 64-bit + integers before using fistpll. Reduce 64-bit values to 63-bit + ones as needed. + * math/libm-test.inc (pow_test): Add more tests. + 2012-03-20 H.J. Lu * debug/backtracesymsfd.c: Include <_itoa.h> instead of diff --git a/math/libm-test.inc b/math/libm-test.inc index af3d645..64d12a5 100644 --- a/math/libm-test.inc +++ b/math/libm-test.inc @@ -5431,11 +5431,75 @@ pow_test (void) TEST_ff_f (pow, 0, -11, plus_infty, DIVIDE_BY_ZERO_EXCEPTION); check_int ("errno for pow(0,-odd) == ERANGE", errno, ERANGE, 0, 0, 0); errno = 0; + TEST_ff_f (pow, 0, -0xffffff, plus_infty, DIVIDE_BY_ZERO_EXCEPTION); + check_int ("errno for pow(0,-odd) == ERANGE", errno, ERANGE, 0, 0, 0); +#ifndef TEST_FLOAT + errno = 0; + TEST_ff_f (pow, 0, -0x1.fffffffffffffp+52L, plus_infty, DIVIDE_BY_ZERO_EXCEPTION); + check_int ("errno for pow(0,-odd) == ERANGE", errno, ERANGE, 0, 0, 0); +#endif +#ifdef TEST_LDOUBLE +# if LDBL_MANT_DIG >= 64 + errno = 0; + TEST_ff_f (pow, 0, -0x1.fffffffffffffffep+63L, plus_infty, DIVIDE_BY_ZERO_EXCEPTION); + check_int ("errno for pow(0,-odd) == ERANGE", errno, ERANGE, 0, 0, 0); +# endif +# if LDBL_MANT_DIG >= 106 + errno = 0; + TEST_ff_f (pow, 0, -0x1.ffffffffffffffffffffffffff8p+105L, plus_infty, DIVIDE_BY_ZERO_EXCEPTION); + check_int ("errno for pow(0,-odd) == ERANGE", errno, ERANGE, 0, 0, 0); +# endif +# if LDBL_MANT_DIG >= 113 + errno = 0; + TEST_ff_f (pow, 0, -0x1.ffffffffffffffffffffffffffffp+112L, plus_infty, DIVIDE_BY_ZERO_EXCEPTION); + check_int ("errno for pow(0,-odd) == ERANGE", errno, ERANGE, 0, 0, 0); +# endif +#endif TEST_ff_f (pow, minus_zero, -1, minus_infty, DIVIDE_BY_ZERO_EXCEPTION); check_int ("errno for pow(-0,-odd) == ERANGE", errno, ERANGE, 0, 0, 0); errno = 0; TEST_ff_f (pow, minus_zero, -11L, minus_infty, DIVIDE_BY_ZERO_EXCEPTION); check_int ("errno for pow(-0,-odd) == ERANGE", errno, ERANGE, 0, 0, 0); + errno = 0; + TEST_ff_f (pow, minus_zero, -0xffffff, minus_infty, DIVIDE_BY_ZERO_EXCEPTION); + check_int ("errno for pow(-0,-odd) == ERANGE", errno, ERANGE, 0, 0, 0); + errno = 0; + TEST_ff_f (pow, minus_zero, -0x1fffffe, plus_infty, DIVIDE_BY_ZERO_EXCEPTION); + check_int ("errno for pow(-0,-even) == ERANGE", errno, ERANGE, 0, 0, 0); +#ifndef TEST_FLOAT + errno = 0; + TEST_ff_f (pow, minus_zero, -0x1.fffffffffffffp+52L, minus_infty, DIVIDE_BY_ZERO_EXCEPTION); + check_int ("errno for pow(-0,-odd) == ERANGE", errno, ERANGE, 0, 0, 0); + errno = 0; + TEST_ff_f (pow, minus_zero, -0x1.fffffffffffffp+53L, plus_infty, DIVIDE_BY_ZERO_EXCEPTION); + check_int ("errno for pow(-0,-even) == ERANGE", errno, ERANGE, 0, 0, 0); +#endif +#ifdef TEST_LDOUBLE +# if LDBL_MANT_DIG >= 64 + errno = 0; + TEST_ff_f (pow, minus_zero, -0x1.fffffffffffffffep+63L, minus_infty, DIVIDE_BY_ZERO_EXCEPTION); + check_int ("errno for pow(-0,-odd) == ERANGE", errno, ERANGE, 0, 0, 0); + errno = 0; + TEST_ff_f (pow, minus_zero, -0x1.fffffffffffffffep+64L, plus_infty, DIVIDE_BY_ZERO_EXCEPTION); + check_int ("errno for pow(-0,-even) == ERANGE", errno, ERANGE, 0, 0, 0); +# endif +# if LDBL_MANT_DIG >= 106 + errno = 0; + TEST_ff_f (pow, minus_zero, -0x1.ffffffffffffffffffffffffff8p+105L, minus_infty, DIVIDE_BY_ZERO_EXCEPTION); + check_int ("errno for pow(-0,-odd) == ERANGE", errno, ERANGE, 0, 0, 0); + errno = 0; + TEST_ff_f (pow, minus_zero, -0x1.ffffffffffffffffffffffffff8p+106L, plus_infty, DIVIDE_BY_ZERO_EXCEPTION); + check_int ("errno for pow(-0,-even) == ERANGE", errno, ERANGE, 0, 0, 0); +# endif +# if LDBL_MANT_DIG >= 113 + errno = 0; + TEST_ff_f (pow, minus_zero, -0x1.ffffffffffffffffffffffffffffp+112L, minus_infty, DIVIDE_BY_ZERO_EXCEPTION); + check_int ("errno for pow(-0,-odd) == ERANGE", errno, ERANGE, 0, 0, 0); + errno = 0; + TEST_ff_f (pow, minus_zero, -0x1.ffffffffffffffffffffffffffffp+113L, plus_infty, DIVIDE_BY_ZERO_EXCEPTION); + check_int ("errno for pow(-0,-even) == ERANGE", errno, ERANGE, 0, 0, 0); +# endif +#endif errno = 0; TEST_ff_f (pow, 0, -2, plus_infty, DIVIDE_BY_ZERO_EXCEPTION); @@ -5444,11 +5508,31 @@ pow_test (void) TEST_ff_f (pow, 0, -11.1L, plus_infty, DIVIDE_BY_ZERO_EXCEPTION); check_int ("errno for pow(0,-num) == ERANGE", errno, ERANGE, 0, 0, 0); errno = 0; + TEST_ff_f (pow, 0, -0x1p24, plus_infty, DIVIDE_BY_ZERO_EXCEPTION); + check_int ("errno for pow(0,-num) == ERANGE", errno, ERANGE, 0, 0, 0); + errno = 0; + TEST_ff_f (pow, 0, -0x1p127, plus_infty, DIVIDE_BY_ZERO_EXCEPTION); + check_int ("errno for pow(0,-num) == ERANGE", errno, ERANGE, 0, 0, 0); + errno = 0; + /* Bug 13879: spurious OVERFLOW exception may be present. */ + TEST_ff_f (pow, 0, -max_value, plus_infty, DIVIDE_BY_ZERO_EXCEPTION|OVERFLOW_EXCEPTION_OK); + check_int ("errno for pow(0,-num) == ERANGE", errno, ERANGE, 0, 0, 0); + errno = 0; TEST_ff_f (pow, minus_zero, -2, plus_infty, DIVIDE_BY_ZERO_EXCEPTION); check_int ("errno for pow(-0,-even) == ERANGE", errno, ERANGE, 0, 0, 0); errno = 0; TEST_ff_f (pow, minus_zero, -11.1L, plus_infty, DIVIDE_BY_ZERO_EXCEPTION); check_int ("errno for pow(-0,-num) == ERANGE", errno, ERANGE, 0, 0, 0); + errno = 0; + TEST_ff_f (pow, minus_zero, -0x1p24, plus_infty, DIVIDE_BY_ZERO_EXCEPTION); + check_int ("errno for pow(-0,-num) == ERANGE", errno, ERANGE, 0, 0, 0); + errno = 0; + TEST_ff_f (pow, minus_zero, -0x1p127, plus_infty, DIVIDE_BY_ZERO_EXCEPTION); + check_int ("errno for pow(-0,-num) == ERANGE", errno, ERANGE, 0, 0, 0); + errno = 0; + /* Bug 13879: spurious OVERFLOW exception may be present. */ + TEST_ff_f (pow, minus_zero, -max_value, plus_infty, DIVIDE_BY_ZERO_EXCEPTION|OVERFLOW_EXCEPTION_OK); + check_int ("errno for pow(-0,-num) == ERANGE", errno, ERANGE, 0, 0, 0); TEST_ff_f (pow, 0x1p72L, 0x1p72L, plus_infty, OVERFLOW_EXCEPTION); TEST_ff_f (pow, 10, -0x1p72L, 0); @@ -5487,32 +5571,155 @@ pow_test (void) /* pow (+inf, y) == +inf for y > 0. */ TEST_ff_f (pow, plus_infty, 2, plus_infty); + TEST_ff_f (pow, plus_infty, 0xffffff, plus_infty); +#ifndef TEST_FLOAT + TEST_ff_f (pow, plus_infty, 0x1.fffffffffffffp+52L, plus_infty); +#endif +#ifdef TEST_LDOUBLE +# if LDBL_MANT_DIG >= 64 + TEST_ff_f (pow, plus_infty, 0x1.fffffffffffffffep+63L, plus_infty); +# endif +# if LDBL_MANT_DIG >= 106 + TEST_ff_f (pow, plus_infty, 0x1.ffffffffffffffffffffffffff8p+105L, plus_infty); +# endif +# if LDBL_MANT_DIG >= 113 + TEST_ff_f (pow, plus_infty, 0x1.ffffffffffffffffffffffffffffp+112L, plus_infty); +# endif +#endif + TEST_ff_f (pow, plus_infty, 0x1p24, plus_infty); + TEST_ff_f (pow, plus_infty, 0x1p127, plus_infty); + TEST_ff_f (pow, plus_infty, max_value, plus_infty); /* pow (+inf, y) == +0 for y < 0. */ TEST_ff_f (pow, plus_infty, -1, 0.0); + TEST_ff_f (pow, plus_infty, -0xffffff, 0.0); +#ifndef TEST_FLOAT + TEST_ff_f (pow, plus_infty, -0x1.fffffffffffffp+52L, 0.0); +#endif +#ifdef TEST_LDOUBLE +# if LDBL_MANT_DIG >= 64 + TEST_ff_f (pow, plus_infty, -0x1.fffffffffffffffep+63L, 0.0); +# endif +# if LDBL_MANT_DIG >= 106 + TEST_ff_f (pow, plus_infty, -0x1.ffffffffffffffffffffffffff8p+105L, 0.0); +# endif +# if LDBL_MANT_DIG >= 113 + TEST_ff_f (pow, plus_infty, -0x1.ffffffffffffffffffffffffffffp+112L, 0.0); +# endif +#endif + TEST_ff_f (pow, plus_infty, -0x1p24, 0.0); + TEST_ff_f (pow, plus_infty, -0x1p127, 0.0); + TEST_ff_f (pow, plus_infty, -max_value, 0.0); /* pow (-inf, y) == -inf for y an odd integer > 0. */ TEST_ff_f (pow, minus_infty, 27, minus_infty); + TEST_ff_f (pow, minus_infty, 0xffffff, minus_infty); + TEST_ff_f (pow, minus_infty, 0x1fffffe, plus_infty); +#ifndef TEST_FLOAT + TEST_ff_f (pow, minus_infty, 0x1.fffffffffffffp+52L, minus_infty); + TEST_ff_f (pow, minus_infty, 0x1.fffffffffffffp+53L, plus_infty); +#endif +#ifdef TEST_LDOUBLE +# if LDBL_MANT_DIG >= 64 + TEST_ff_f (pow, minus_infty, 0x1.fffffffffffffffep+63L, minus_infty); + TEST_ff_f (pow, minus_infty, 0x1.fffffffffffffffep+64L, plus_infty); +# endif +# if LDBL_MANT_DIG >= 106 + TEST_ff_f (pow, minus_infty, 0x1.ffffffffffffffffffffffffff8p+105L, minus_infty); + TEST_ff_f (pow, minus_infty, 0x1.ffffffffffffffffffffffffff8p+106L, plus_infty); +# endif +# if LDBL_MANT_DIG >= 113 + TEST_ff_f (pow, minus_infty, 0x1.ffffffffffffffffffffffffffffp+112L, minus_infty); + TEST_ff_f (pow, minus_infty, 0x1.ffffffffffffffffffffffffffffp+113L, plus_infty); +# endif +#endif /* pow (-inf, y) == +inf for y > 0 and not an odd integer. */ TEST_ff_f (pow, minus_infty, 28, plus_infty); + TEST_ff_f (pow, minus_infty, 0x1p24, plus_infty); + TEST_ff_f (pow, minus_infty, 0x1p127, plus_infty); + TEST_ff_f (pow, minus_infty, max_value, plus_infty); /* pow (-inf, y) == -0 for y an odd integer < 0. */ TEST_ff_f (pow, minus_infty, -3, minus_zero); + TEST_ff_f (pow, minus_infty, -0xffffff, minus_zero); + TEST_ff_f (pow, minus_infty, -0x1fffffe, plus_zero); +#ifndef TEST_FLOAT + TEST_ff_f (pow, minus_infty, -0x1.fffffffffffffp+52L, minus_zero); + TEST_ff_f (pow, minus_infty, -0x1.fffffffffffffp+53L, plus_zero); +#endif +#ifdef TEST_LDOUBLE +# if LDBL_MANT_DIG >= 64 + TEST_ff_f (pow, minus_infty, -0x1.fffffffffffffffep+63L, minus_zero); + TEST_ff_f (pow, minus_infty, -0x1.fffffffffffffffep+64L, plus_zero); +# endif +# if LDBL_MANT_DIG >= 106 + TEST_ff_f (pow, minus_infty, -0x1.ffffffffffffffffffffffffff8p+105L, minus_zero); + TEST_ff_f (pow, minus_infty, -0x1.ffffffffffffffffffffffffff8p+106L, plus_zero); +# endif +# if LDBL_MANT_DIG >= 113 + TEST_ff_f (pow, minus_infty, -0x1.ffffffffffffffffffffffffffffp+112L, minus_zero); + TEST_ff_f (pow, minus_infty, -0x1.ffffffffffffffffffffffffffffp+113L, plus_zero); +# endif +#endif /* pow (-inf, y) == +0 for y < 0 and not an odd integer. */ TEST_ff_f (pow, minus_infty, -2.0, 0.0); + TEST_ff_f (pow, minus_infty, -0x1p24, 0.0); + TEST_ff_f (pow, minus_infty, -0x1p127, 0.0); + TEST_ff_f (pow, minus_infty, -max_value, 0.0); /* pow (+0, y) == +0 for y an odd integer > 0. */ TEST_ff_f (pow, 0.0, 27, 0.0); + TEST_ff_f (pow, 0.0, 0xffffff, 0.0); +#ifndef TEST_FLOAT + TEST_ff_f (pow, 0.0, 0x1.fffffffffffffp+52L, 0.0); +#endif +#ifdef TEST_LDOUBLE +# if LDBL_MANT_DIG >= 64 + TEST_ff_f (pow, 0.0, 0x1.fffffffffffffffep+63L, 0.0); +# endif +# if LDBL_MANT_DIG >= 106 + TEST_ff_f (pow, 0.0, 0x1.ffffffffffffffffffffffffff8p+105L, 0.0); +# endif +# if LDBL_MANT_DIG >= 113 + TEST_ff_f (pow, 0.0, 0x1.ffffffffffffffffffffffffffffp+112L, 0.0); +# endif +#endif /* pow (-0, y) == -0 for y an odd integer > 0. */ TEST_ff_f (pow, minus_zero, 27, minus_zero); + TEST_ff_f (pow, minus_zero, 0xffffff, minus_zero); + TEST_ff_f (pow, minus_zero, 0x1fffffe, plus_zero); +#ifndef TEST_FLOAT + TEST_ff_f (pow, minus_zero, 0x1.fffffffffffffp+52L, minus_zero); + TEST_ff_f (pow, minus_zero, 0x1.fffffffffffffp+53L, plus_zero); +#endif +#ifdef TEST_LDOUBLE +# if LDBL_MANT_DIG >= 64 + TEST_ff_f (pow, minus_zero, 0x1.fffffffffffffffep+63L, minus_zero); + TEST_ff_f (pow, minus_zero, 0x1.fffffffffffffffep+64L, plus_zero); +# endif +# if LDBL_MANT_DIG >= 106 + TEST_ff_f (pow, minus_zero, 0x1.ffffffffffffffffffffffffff8p+105L, minus_zero); + TEST_ff_f (pow, minus_zero, 0x1.ffffffffffffffffffffffffff8p+106L, plus_zero); +# endif +# if LDBL_MANT_DIG >= 113 + TEST_ff_f (pow, minus_zero, 0x1.ffffffffffffffffffffffffffffp+112L, minus_zero); + TEST_ff_f (pow, minus_zero, 0x1.ffffffffffffffffffffffffffffp+113L, plus_zero); +# endif +#endif /* pow (+0, y) == +0 for y > 0 and not an odd integer. */ TEST_ff_f (pow, 0.0, 4, 0.0); + TEST_ff_f (pow, 0.0, 0x1p24, 0.0); + TEST_ff_f (pow, 0.0, 0x1p127, 0.0); + TEST_ff_f (pow, 0.0, max_value, 0.0); /* pow (-0, y) == +0 for y > 0 and not an odd integer. */ TEST_ff_f (pow, minus_zero, 4, 0.0); + TEST_ff_f (pow, minus_zero, 0x1p24, 0.0); + TEST_ff_f (pow, minus_zero, 0x1p127, 0.0); + TEST_ff_f (pow, minus_zero, max_value, 0.0); TEST_ff_f (pow, 16, 0.25L, 2); TEST_ff_f (pow, 0x1p64L, 0.125L, 256); diff --git a/sysdeps/i386/fpu/e_pow.S b/sysdeps/i386/fpu/e_pow.S index 63c44f1..1abedf6 100644 --- a/sysdeps/i386/fpu/e_pow.S +++ b/sysdeps/i386/fpu/e_pow.S @@ -230,6 +230,16 @@ ENTRY(__ieee754_pow) testb $2, %dh jz 16f // jump if x == +inf + // fistpll raises invalid exception for |y| >= 1L<<63, so test + // that (in which case y is certainly even) before testing + // whether y is odd. + fld %st // y : y + fabs // |y| : y + fcompl MO(p63) // y + fnstsw + sahf + jnc 16f + // We must find out whether y is an odd integer. fld %st // y : y fistpll (%esp) // y @@ -239,20 +249,13 @@ ENTRY(__ieee754_pow) sahf jne 17f - // OK, the value is an integer, but is the number of bits small - // enough so that all are coming from the mantissa? + // OK, the value is an integer. popl %eax cfi_adjust_cfa_offset (-4) popl %edx cfi_adjust_cfa_offset (-4) andb $1, %al jz 18f // jump if not odd - movl %edx, %eax - orl %edx, %edx - jns 155f - negl %eax -155: cmpl $0x00200000, %eax - ja 18f // does not fit in mantissa bits // It's an odd integer. shrl $31, %edx fldl MOX(minf_mzero, %edx, 8) @@ -289,6 +292,16 @@ ENTRY(__ieee754_pow) testb $2, %dh jz 25f + // fistpll raises invalid exception for |y| >= 1L<<63, so test + // that (in which case y is certainly even) before testing + // whether y is odd. + fld %st // y : y + fabs // |y| : y + fcompl MO(p63) // y + fnstsw + sahf + jnc 25f + fld %st // y : y fistpll (%esp) // y fildll (%esp) // int(y) : y @@ -297,16 +310,13 @@ ENTRY(__ieee754_pow) sahf jne 26f - // OK, the value is an integer, but is the number of bits small - // enough so that all are coming from the mantissa? + // OK, the value is an integer. popl %eax cfi_adjust_cfa_offset (-4) popl %edx cfi_adjust_cfa_offset (-4) andb $1, %al jz 27f // jump if not odd - cmpl $0xffe00000, %edx - jbe 27f // does not fit in mantissa bits // It's an odd integer. // Raise divide-by-zero exception and get minus infinity value. fldl MO(one) @@ -329,6 +339,14 @@ ENTRY(__ieee754_pow) 21: testb $2, %dh jz 22f + // fistpll raises invalid exception for |y| >= 1L<<63, so test + // that (in which case y is certainly even) before testing + // whether y is odd. + fcoml MO(p63) // y + fnstsw + sahf + jnc 22f + fld %st // y : y fistpll (%esp) // y fildll (%esp) // int(y) : y @@ -337,16 +355,13 @@ ENTRY(__ieee754_pow) sahf jne 23f - // OK, the value is an integer, but is the number of bits small - // enough so that all are coming from the mantissa? + // OK, the value is an integer. popl %eax cfi_adjust_cfa_offset (-4) popl %edx cfi_adjust_cfa_offset (-4) andb $1, %al jz 24f // jump if not odd - cmpl $0xffe00000, %edx - jae 24f // does not fit in mantissa bits // It's an odd integer. fldl MO(mzero) ret diff --git a/sysdeps/i386/fpu/e_powf.S b/sysdeps/i386/fpu/e_powf.S index cc8456d..aa58ed2 100644 --- a/sysdeps/i386/fpu/e_powf.S +++ b/sysdeps/i386/fpu/e_powf.S @@ -224,6 +224,16 @@ ENTRY(__ieee754_powf) testb $2, %dh jz 16f // jump if x == +inf + // fistpl raises invalid exception for |y| >= 1L<<31, so test + // that (in which case y is certainly even) before testing + // whether y is odd. + fld %st // y : y + fabs // |y| : y + fcompl MO(p31) // y + fnstsw + sahf + jnc 16f + // We must find out whether y is an odd integer. fld %st // y : y fistpl (%esp) // y @@ -233,18 +243,11 @@ ENTRY(__ieee754_powf) sahf jne 17f - // OK, the value is an integer, but is the number of bits small - // enough so that all are coming from the mantissa? + // OK, the value is an integer. popl %edx cfi_adjust_cfa_offset (-4) testb $1, %dl jz 18f // jump if not odd - movl %edx, %eax - orl %edx, %edx - jns 155f - negl %eax -155: cmpl $0x01000000, %eax - ja 18f // does not fit in mantissa bits // It's an odd integer. shrl $31, %edx fldl MOX(minf_mzero, %edx, 8) @@ -281,6 +284,16 @@ ENTRY(__ieee754_powf) testb $2, %dh jz 25f + // fistpl raises invalid exception for |y| >= 1L<<31, so test + // that (in which case y is certainly even) before testing + // whether y is odd. + fld %st // y : y + fabs // |y| : y + fcompl MO(p31) // y + fnstsw + sahf + jnc 25f + fld %st // y : y fistpl (%esp) // y fildl (%esp) // int(y) : y @@ -289,14 +302,11 @@ ENTRY(__ieee754_powf) sahf jne 26f - // OK, the value is an integer, but is the number of bits small - // enough so that all are coming from the mantissa? + // OK, the value is an integer. popl %edx cfi_adjust_cfa_offset (-4) testb $1, %dl jz 27f // jump if not odd - cmpl $0xff000000, %edx - jbe 27f // does not fit in mantissa bits // It's an odd integer. // Raise divide-by-zero exception and get minus infinity value. fldl MO(one) @@ -319,6 +329,14 @@ ENTRY(__ieee754_powf) 21: testb $2, %dh jz 22f + // fistpl raises invalid exception for |y| >= 1L<<31, so test + // that (in which case y is certainly even) before testing + // whether y is odd. + fcoml MO(p31) // y + fnstsw + sahf + jnc 22f + fld %st // y : y fistpl (%esp) // y fildl (%esp) // int(y) : y @@ -327,14 +345,11 @@ ENTRY(__ieee754_powf) sahf jne 23f - // OK, the value is an integer, but is the number of bits small - // enough so that all are coming from the mantissa? + // OK, the value is an integer. popl %edx cfi_adjust_cfa_offset (-4) testb $1, %dl jz 24f // jump if not odd - cmpl $0xff000000, %edx - jae 24f // does not fit in mantissa bits // It's an odd integer. fldl MO(mzero) ret diff --git a/sysdeps/i386/fpu/e_powl.S b/sysdeps/i386/fpu/e_powl.S index 5d85089..c0aa194c 100644 --- a/sysdeps/i386/fpu/e_powl.S +++ b/sysdeps/i386/fpu/e_powl.S @@ -32,6 +32,9 @@ limit: .double 0.29 ASM_TYPE_DIRECTIVE(p63,@object) p63: .byte 0, 0, 0, 0, 0, 0, 0xe0, 0x43 ASM_SIZE_DIRECTIVE(p63) + ASM_TYPE_DIRECTIVE(p64,@object) +p64: .byte 0, 0, 0, 0, 0, 0, 0xf0, 0x43 + ASM_SIZE_DIRECTIVE(p64) .section .rodata.cst16,"aM",@progbits,16 @@ -243,6 +246,19 @@ ENTRY(__ieee754_powl) testb $2, %dh jz 16f // jump if x == +inf + // fistpll raises invalid exception for |y| >= 1L<<63, but y + // may be odd unless we know |y| >= 1L<<64. + fld %st // y : y + fabs // |y| : y + fcompl MO(p64) // y + fnstsw + sahf + jnc 16f + fldl MO(p63) // p63 : y + fxch // y : p63 + fprem // y%p63 : p63 + fstp %st(1) // y%p63 + // We must find out whether y is an odd integer. fld %st // y : y fistpll (%esp) // y @@ -295,6 +311,19 @@ ENTRY(__ieee754_powl) testb $2, %dh jz 25f + // fistpll raises invalid exception for |y| >= 1L<<63, but y + // may be odd unless we know |y| >= 1L<<64. + fld %st // y : y + fabs // |y| : y + fcompl MO(p64) // y + fnstsw + sahf + jnc 25f + fldl MO(p63) // p63 : y + fxch // y : p63 + fprem // y%p63 : p63 + fstp %st(1) // y%p63 + fld %st // y : y fistpll (%esp) // y fildll (%esp) // int(y) : y @@ -332,6 +361,18 @@ ENTRY(__ieee754_powl) 21: testb $2, %dh jz 22f + // fistpll raises invalid exception for |y| >= 1L<<63, but y + // may be odd unless we know |y| >= 1L<<64. + fld %st // y : y + fcompl MO(p64) // y + fnstsw + sahf + jnc 22f + fldl MO(p63) // p63 : y + fxch // y : p63 + fprem // y%p63 : p63 + fstp %st(1) // y%p63 + fld %st // y : y fistpll (%esp) // y fildll (%esp) // int(y) : y diff --git a/sysdeps/ieee754/dbl-64/e_pow.c b/sysdeps/ieee754/dbl-64/e_pow.c index f936a72..26ffaad 100644 --- a/sysdeps/ieee754/dbl-64/e_pow.c +++ b/sysdeps/ieee754/dbl-64/e_pow.c @@ -111,7 +111,7 @@ __ieee754_pow(double x, double y) { if (((v.i[HIGH_HALF] & 0x7fffffff) == 0x7ff00000 && v.i[LOW_HALF] != 0) || (v.i[HIGH_HALF] & 0x7fffffff) > 0x7ff00000) return y; - if (ABS(y) > 1.0e20) return (y>0)?0:INF.x; + if (ABS(y) > 1.0e20) return (y>0)?0:1.0/ABS(x); k = checkint(y); if (k == -1) return y < 0 ? 1.0/x : x; diff --git a/sysdeps/x86_64/fpu/e_powl.S b/sysdeps/x86_64/fpu/e_powl.S index bd6d828..562791d 100644 --- a/sysdeps/x86_64/fpu/e_powl.S +++ b/sysdeps/x86_64/fpu/e_powl.S @@ -32,6 +32,9 @@ limit: .double 0.29 ASM_TYPE_DIRECTIVE(p63,@object) p63: .byte 0, 0, 0, 0, 0, 0, 0xe0, 0x43 ASM_SIZE_DIRECTIVE(p63) + ASM_TYPE_DIRECTIVE(p64,@object) +p64: .byte 0, 0, 0, 0, 0, 0, 0xf0, 0x43 + ASM_SIZE_DIRECTIVE(p64) .section .rodata.cst16,"aM",@progbits,16 @@ -227,6 +230,19 @@ ENTRY(__ieee754_powl) testb $2, %dh jz 16f // jump if x == +inf + // fistpll raises invalid exception for |y| >= 1L<<63, but y + // may be odd unless we know |y| >= 1L<<64. + fldl MO(p64) // 1L<<64 : y + fld %st(1) // y : 1L<<64 : y + fabs // |y| : 1L<<64 : y + fcomip %st(1), %st // 1L<<64 : y + fstp %st(0) // y + jnc 16f + fldl MO(p63) // p63 : y + fxch // y : p63 + fprem // y%p63 : p63 + fstp %st(1) // y%p63 + // We must find out whether y is an odd integer. fld %st // y : y fistpll -8(%rsp) // y @@ -284,6 +300,19 @@ ENTRY(__ieee754_powl) testb $2, %dh jz 25f + // fistpll raises invalid exception for |y| >= 1L<<63, but y + // may be odd unless we know |y| >= 1L<<64. + fldl MO(p64) // 1L<<64 : y + fld %st(1) // y : 1L<<64 : y + fabs // |y| : 1L<<64 : y + fcomip %st(1), %st // 1L<<64 : y + fstp %st(0) // y + jnc 25f + fldl MO(p63) // p63 : y + fxch // y : p63 + fprem // y%p63 : p63 + fstp %st(1) // y%p63 + fld %st // y : y fistpll -8(%rsp) // y fildll -8(%rsp) // int(y) : y @@ -315,6 +344,18 @@ ENTRY(__ieee754_powl) 21: testb $2, %dh jz 22f + // fistpll raises invalid exception for |y| >= 1L<<63, but y + // may be odd unless we know |y| >= 1L<<64. + fldl MO(p64) // 1L<<64 : y + fxch // y : 1L<<64 + fcomi %st(1), %st // y : 1L<<64 + fstp %st(1) // y + jnc 22f + fldl MO(p63) // p63 : y + fxch // y : p63 + fprem // y%p63 : p63 + fstp %st(1) // y%p63 + fld %st // y : y fistpll -8(%rsp) // y fildll -8(%rsp) // int(y) : y -- 2.7.4