Fix pow of zero and infinity to large powers.

diff --git a/ChangeLog b/ChangeLog
index edac2dc..d548878 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,24 @@
+2012-03-20  Joseph Myers  <joseph@codesourcery.com>
+
+       [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  <hongjiu.lu@intel.com>
 
        * debug/backtracesymsfd.c: Include <_itoa.h> instead of

diff --git a/math/libm-test.inc b/math/libm-test.inc
index af3d645..64d12a5 100644 (file)
--- 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 (file)
--- 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 (file)
--- 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..c0aa194 100644 (file)
--- 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 (file)
--- 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 (file)
--- 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