Fix pow of negative numbers to integer exponents (bugs 369, 2678, 3866).

diff --git a/ChangeLog b/ChangeLog
index ca2fa68..f75e121 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,5 +1,20 @@
 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.
+
        [BZ #3868]
        [BZ #13879]
        [BZ #13910]

diff --git a/NEWS b/NEWS
index c2f56a8..3ebc0ff 100644 (file)
--- a/NEWS
+++ b/NEWS
@@ -9,16 +9,16 @@ Version 2.16
 
 * The following bugs are resolved with this release:
 
-  174, 350, 411, 2541, 2547, 2548, 2551, 2552, 2553, 2554, 2562, 2563, 2565,
-  2566, 2576, 3335, 3868, 3976, 3992, 4026, 4108, 4596, 4822, 5077, 5461,
-  5805, 5993, 6471, 6730, 6884, 6907, 6911, 9739, 9902, 10110, 10135, 10140,
-  10210, 10545, 10716, 11174, 11322, 11365, 11451, 11494, 12047, 13058,
-  13525, 13526, 13527, 13528, 13529, 13530, 13531, 13532, 13533, 13547,
-  13551, 13552, 13553, 13555, 13559, 13566, 13583, 13618, 13637, 13656,
-  13658, 13673, 13695, 13704, 13706, 13726, 13738, 13786, 13792, 13806,
-  13824, 13840, 13841, 13844, 13846, 13851, 13852, 13854, 13871, 13879,
-  13883, 13892, 13910, 13911, 13912, 13913, 13915, 13916, 13917, 13918,
-  13919, 13920, 13921
+  174, 350, 369, 411, 2541, 2547, 2548, 2551, 2552, 2553, 2554, 2562, 2563,
+  2565, 2566, 2576, 2678, 3335, 3866, 3868, 3976, 3992, 4026, 4108, 4596,
+  4822, 5077, 5461, 5805, 5993, 6471, 6730, 6884, 6907, 6911, 9739, 9902,
+  10110, 10135, 10140, 10210, 10545, 10716, 11174, 11322, 11365, 11451,
+  11494, 12047, 13058, 13525, 13526, 13527, 13528, 13529, 13530, 13531,
+  13532, 13533, 13547, 13551, 13552, 13553, 13555, 13559, 13566, 13583,
+  13618, 13637, 13656, 13658, 13673, 13695, 13704, 13706, 13726, 13738,
+  13786, 13792, 13806, 13824, 13840, 13841, 13844, 13846, 13851, 13852,
+  13854, 13871, 13879, 13883, 13892, 13910, 13911, 13912, 13913, 13915,
+  13916, 13917, 13918, 13919, 13920, 13921
 
 * ISO C11 support:
 

diff --git a/math/libm-test.inc b/math/libm-test.inc
index 68f6ef2..32bce45 100644 (file)
--- 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 (file)
--- 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 (file)
--- 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 c0aa194..0e7c05b 100644 (file)
--- 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 (file)
--- 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 (file)
--- 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)