Fix spurious overflow exceptions from x86/x86_64 powl (bug 13872).

diff --git a/ChangeLog b/ChangeLog
index 3686491a2d6dfbfb48ce5027674cbaa1afc73fc0..a85a1c06ef6715ba75c3be29b468a9ed67797159 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,5 +1,13 @@
 2012-04-09  Joseph Myers  <joseph@codesourcery.com>
 
+       [BZ #13872]
+       * sysdeps/i386/fpu/e_powl.S (p78): New object.
+       (__ieee754_powl): Saturate large exponents rather than testing for
+       overflow of y*log2(x).
+       * sysdeps/x86_64/fpu/e_powl.S: Likewise.
+       * math/libm-test.inc (pow_test): Do not permit spurious overflow
+       exceptions.
+
        [BZ #11521]
        * math/s_ctan.c: Include <float.h>.
        (__ctan): Avoid internal overflow or cancellation in calculating

diff --git a/NEWS b/NEWS
index bb303f8b9320b2887db51340e3260b9d7ebfc122..357f3b4e697ca1793171f1449525720fa54464cd 100644 (file)
--- a/NEWS
+++ b/NEWS
@@ -18,9 +18,10 @@ Version 2.16
   13530, 13531, 13532, 13533, 13547, 13551, 13552, 13553, 13555, 13559,
   13566, 13583, 13592, 13618, 13637, 13656, 13658, 13673, 13691, 13695,
   13704, 13705, 13706, 13726, 13738, 13760, 13761, 13786, 13792, 13806,
-  13824, 13840, 13841, 13844, 13846, 13851, 13852, 13854, 13871, 13873,
-  13879, 13883, 13892, 13895, 13908, 13910, 13911, 13912, 13913, 13915,
-  13916, 13917, 13918, 13919, 13920, 13921, 13926, 13928, 13938, 13963
+  13824, 13840, 13841, 13844, 13846, 13851, 13852, 13854, 13871, 13872,
+  13873, 13879, 13883, 13892, 13895, 13908, 13910, 13911, 13912, 13913,
+  13915, 13916, 13917, 13918, 13919, 13920, 13921, 13926, 13928, 13938,
+  13963
 
 * ISO C11 support:
 

diff --git a/math/libm-test.inc b/math/libm-test.inc
index a551b9f3f4626ddb54ecabfe61c20da3df623456..2809d57d094d83a03a4cdcc98481f2750a62dcc0 100644 (file)
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -5682,8 +5682,7 @@ pow_test (void)
   TEST_ff_f (pow, 0x1p72L, 0x1p72L, plus_infty, OVERFLOW_EXCEPTION);
   TEST_ff_f (pow, 10, -0x1p72L, 0);
   TEST_ff_f (pow, max_value, max_value, plus_infty, OVERFLOW_EXCEPTION);
-  /* Bug 13872: spurious OVERFLOW exception may be present.  */
-  TEST_ff_f (pow, 10, -max_value, 0, OVERFLOW_EXCEPTION_OK);
+  TEST_ff_f (pow, 10, -max_value, 0);
 
   TEST_ff_f (pow, 0, 1, 0);
   TEST_ff_f (pow, 0, 11, 0);
@@ -5997,8 +5996,7 @@ pow_test (void)
   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, -max_value, plus_zero);
 
   TEST_ff_f (pow, -max_value, 0xffffff, minus_infty, OVERFLOW_EXCEPTION);
   TEST_ff_f (pow, -max_value, 0x1fffffe, plus_infty, OVERFLOW_EXCEPTION);
@@ -6122,8 +6120,7 @@ pow_test (void)
   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);
+  TEST_ff_f (pow, -min_value, max_value, plus_zero);
 
 #ifndef TEST_LDOUBLE /* Bug 13881.  */
   TEST_ff_f (pow, 0x0.ffffffp0, 10, 0.999999403953712118183885036774764444747L);

diff --git a/sysdeps/i386/fpu/e_powl.S b/sysdeps/i386/fpu/e_powl.S
index 0e7c05bb82f4e8b88e6213faa4c23c3fd5787a31..5b166eab4b9e786969df0580ceddd28f2b755aeb 100644 (file)
--- a/sysdeps/i386/fpu/e_powl.S
+++ b/sysdeps/i386/fpu/e_powl.S
@@ -35,6 +35,9 @@ p63:  .byte 0, 0, 0, 0, 0, 0, 0xe0, 0x43
        ASM_TYPE_DIRECTIVE(p64,@object)
 p64:   .byte 0, 0, 0, 0, 0, 0, 0xf0, 0x43
        ASM_SIZE_DIRECTIVE(p64)
+       ASM_TYPE_DIRECTIVE(p78,@object)
+p78:   .byte 0, 0, 0, 0, 0, 0, 0xd0, 0x44
+       ASM_SIZE_DIRECTIVE(p78)
 
        .section .rodata.cst16,"aM",@progbits,16
 
@@ -166,6 +169,21 @@ ENTRY(__ieee754_powl)
        fxch                    // x : y
        fabs                    // |x| : y
        fxch                    // y : |x|
+       // If y has absolute value at least 1L<<78, then any finite
+       // nonzero x will result in 0 (underflow), 1 or infinity (overflow).
+       // Saturate y to those bounds to avoid overflow in the calculation
+       // of y*log2(x).
+       fld     %st             // y : y : |x|
+       fabs                    // |y| : y : |x|
+       fcompl  MO(p78)         // y : |x|
+       fnstsw
+       sahf
+       jc      3f
+       fstp    %st(0)          // pop y
+       fldl    MO(p78)         // 1L<<78 : |x|
+       testb   $2, %dl
+       jz      3f              // y > 0
+       fchs                    // -(1L<<78) : |x|
        .align ALIGNARG(4)
 3:     /* y is a real number.  */
        fxch                    // x : y
@@ -185,11 +203,6 @@ ENTRY(__ieee754_powl)
 
 7:     fyl2x                   // log2(x) : y
 8:     fmul    %st(1)          // y*log2(x) : y
-       fxam
-       fnstsw
-       andb    $0x45, %ah
-       cmpb    $0x05, %ah      // is y*log2(x) == ±inf ?
-       je      28f
        fst     %st(1)          // y*log2(x) : y*log2(x)
        frndint                 // int(y*log2(x)) : y*log2(x)
        fsubr   %st, %st(1)     // int(y*log2(x)) : fract(y*log2(x))
@@ -198,13 +211,7 @@ 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))
-       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))
-29:    testb   $2, %dh
+       testb   $2, %dh
        jz      292f
        // x is negative.  If y is an odd integer, negate the result.
        fldt    24(%esp)        // y : abs(result)

diff --git a/sysdeps/x86_64/fpu/e_powl.S b/sysdeps/x86_64/fpu/e_powl.S
index 0626ce41722f064a2abfce8d5690478b5eb8ffa9..10ede22648d70086568d9e5c5e59738893a0f5dd 100644 (file)
--- a/sysdeps/x86_64/fpu/e_powl.S
+++ b/sysdeps/x86_64/fpu/e_powl.S
@@ -35,6 +35,9 @@ p63:  .byte 0, 0, 0, 0, 0, 0, 0xe0, 0x43
        ASM_TYPE_DIRECTIVE(p64,@object)
 p64:   .byte 0, 0, 0, 0, 0, 0, 0xf0, 0x43
        ASM_SIZE_DIRECTIVE(p64)
+       ASM_TYPE_DIRECTIVE(p78,@object)
+p78:   .byte 0, 0, 0, 0, 0, 0, 0xd0, 0x44
+       ASM_SIZE_DIRECTIVE(p78)
 
        .section .rodata.cst16,"aM",@progbits,16
 
@@ -151,6 +154,21 @@ ENTRY(__ieee754_powl)
        fxch                    // x : y
        fabs                    // |x| : y
        fxch                    // y : |x|
+       // If y has absolute value at least 1L<<78, then any finite
+       // nonzero x will result in 0 (underflow), 1 or infinity (overflow).
+       // Saturate y to those bounds to avoid overflow in the calculation
+       // of y*log2(x).
+       fldl    MO(p78)         // 1L<<78 : y : |x|
+       fld     %st(1)          // y : 1L<<78 : y : |x|
+       fabs                    // |y| : 1L<<78 : y : |x|
+       fcomip  %st(1), %st     // 1L<<78 : y : |x|
+       fstp    %st(0)          // y : |x|
+       jc      3f
+       fstp    %st(0)          // pop y
+       fldl    MO(p78)         // 1L<<78 : |x|
+       testb   $2, %dl
+       jz      3f              // y > 0
+       fchs                    // -(1L<<78) : |x|
        .align ALIGNARG(4)
 3:     /* y is a real number.  */
        fxch                    // x : y
@@ -170,11 +188,6 @@ ENTRY(__ieee754_powl)
 
 7:     fyl2x                   // log2(x) : y
 8:     fmul    %st(1)          // y*log2(x) : y
-       fxam
-       fnstsw
-       andb    $0x45, %ah
-       cmpb    $0x05, %ah      // is y*log2(x) == ±inf ?
-       je      28f
        fst     %st(1)          // y*log2(x) : y*log2(x)
        frndint                 // int(y*log2(x)) : y*log2(x)
        fsubr   %st, %st(1)     // int(y*log2(x)) : fract(y*log2(x))
@@ -183,13 +196,7 @@ 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))
-       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))
-29:    testb   $2, %dh
+       testb   $2, %dh
        jz      292f
        // x is negative.  If y is an odd integer, negate the result.
        fldt    24(%rsp)        // y : abs(result)