From f88acd39da2a509081e541b84ecbf204ef20f9e8 Mon Sep 17 00:00:00 2001 From: Joseph Myers Date: Thu, 19 Dec 2013 13:36:10 +0000 Subject: Fix x86/x86_64 expm1 inaccuracy near 0 in directed rounding modes (bug 16293). Bug 16293 is inaccuracy of x86/x86_64 versions of expm1, near 0 in directed rounding modes, that arises from frndint rounding the exponent to 1 or -1 instead of 0, resulting in large cancellation error. This inaccuracy in turn affects other functions such as sinh that use expm1. This patch fixes the problem by setting round-to-nearest mode temporarily around the affected calls to frndint. I don't think this is needed for other uses of frndint, such as in exp itself, as only for expm1 is the cancellation error significant. Tested x86_64 and x86 and ulps updated accordingly. * sysdeps/i386/fpu/e_expl.S (IEEE754_EXPL) [USE_AS_EXPM1L]: Set round-to-nearest mode when using frndint. * sysdeps/i386/fpu/s_expm1.S (__expm1): Likewise. * sysdeps/i386/fpu/s_expm1f.S (__expm1f): Likewise. * sysdeps/x86_64/fpu/e_expl.S (IEEE754_EXPL) [USE_AS_EXPM1L]: Likewise. * math/auto-libm-test-in: Add more tests of expm1. Do not expect sinh test to fail. * math/auto-libm-test-out: Regenerated. * math/libm-test.inc (TEST_COND_x86_64): Remove macro. (TEST_COND_x86): Likewise. (expm1_tonearest_test_data): New array. (expm1_test_tonearest): New function. (expm1_towardzero_test_data): New array. (expm1_test_towardzero): New function. (expm1_downward_test_data): New array. (expm1_test_downward): New function. (expm1_upward_test_data): New array. (expm1_test_upward): New function. (main): Run the new test functions. * sysdeps/i386/fpu/libm-test-ulps: Update. * sysdeps/x86_64/fpu/libm-test-ulps: Likewise. --- sysdeps/x86_64/fpu/e_expl.S | 11 ++ sysdeps/x86_64/fpu/libm-test-ulps | 319 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 330 insertions(+) (limited to 'sysdeps/x86_64') diff --git a/sysdeps/x86_64/fpu/e_expl.S b/sysdeps/x86_64/fpu/e_expl.S index a919780390..1c21f03ddc 100644 --- a/sysdeps/x86_64/fpu/e_expl.S +++ b/sysdeps/x86_64/fpu/e_expl.S @@ -127,9 +127,20 @@ ENTRY(IEEE754_EXPL) #endif 3: FLDLOG /* 1 log2(base) */ fmul %st(1), %st /* 1 x log2(base) */ +#ifdef USE_AS_EXPM1L + /* Set round-to-nearest temporarily. */ + fstcw -4(%rsp) + movl $0xf3ff, %edx + andl -4(%rsp), %edx + movl %edx, -8(%rsp) + fldcw -8(%rsp) +#endif frndint /* 1 i */ fld %st(1) /* 2 x */ frndint /* 2 xi */ +#ifdef USE_AS_EXPM1L + fldcw -4(%rsp) +#endif fld %st(1) /* 3 i */ fldt MO(c0) /* 4 c0 */ fld %st(2) /* 5 xi */ diff --git a/sysdeps/x86_64/fpu/libm-test-ulps b/sysdeps/x86_64/fpu/libm-test-ulps index b45ce1dfcd..89b4bc8289 100644 --- a/sysdeps/x86_64/fpu/libm-test-ulps +++ b/sysdeps/x86_64/fpu/libm-test-ulps @@ -7011,9 +7011,15 @@ float: 1 ifloat: 1 # expm1 +Test "expm1 (-0x1p-64)": +ildouble: 1 +ldouble: 1 Test "expm1 (-0x2.dp+4)": ildouble: 1 ldouble: 1 +Test "expm1 (-0x4p-12)": +ildouble: 1 +ldouble: 1 Test "expm1 (-45.0)": ildouble: 1 ldouble: 1 @@ -7048,6 +7054,281 @@ Test "expm1 (500.0)": double: 1 idouble: 1 +# expm1_downward +Test "expm1_downward (-0x1p-100)": +ildouble: 1 +ldouble: 1 +Test "expm1_downward (-0x2.ep+4)": +ildouble: 1 +ldouble: 1 +Test "expm1_downward (-0x4.9p+4)": +ildouble: 1 +ldouble: 1 +Test "expm1_downward (-0x4.bp+4)": +ildouble: 1 +ldouble: 1 +Test "expm1_downward (-0x4p-4)": +ildouble: 1 +ldouble: 1 +Test "expm1_downward (-0x5p+4)": +ildouble: 1 +ldouble: 1 +Test "expm1_downward (-0x6.4p+4)": +ildouble: 1 +ldouble: 1 +Test "expm1_downward (0x1.f4p+8)": +double: 1 +idouble: 1 +Test "expm1_downward (0x1p+0)": +ildouble: 1 +ldouble: 1 +Test "expm1_downward (0x1p-100)": +ildouble: 1 +ldouble: 1 +Test "expm1_downward (0x1p-32)": +ildouble: 1 +ldouble: 1 +Test "expm1_downward (0x3.2p+4)": +double: 1 +float: 1 +idouble: 1 +ifloat: 1 +ildouble: 1 +ldouble: 1 +Test "expm1_downward (0x4p-12)": +ildouble: 1 +ldouble: 1 +Test "expm1_downward (0x4p-52)": +ildouble: 1 +ldouble: 1 +Test "expm1_downward (0x7.fp+4)": +double: 1 +idouble: 1 +ildouble: 1 +ldouble: 1 +Test "expm1_downward (0x8p-32)": +ildouble: 1 +ldouble: 1 + +# expm1_tonearest +Test "expm1_tonearest (-0x1p-64)": +ildouble: 1 +ldouble: 1 +Test "expm1_tonearest (-0x2.dp+4)": +ildouble: 1 +ldouble: 1 +Test "expm1_tonearest (-0x4p-12)": +ildouble: 1 +ldouble: 1 +Test "expm1_tonearest (0x1.f4p+8)": +double: 1 +idouble: 1 +Test "expm1_tonearest (0x1p+0)": +double: 1 +float: 1 +idouble: 1 +ifloat: 1 +ildouble: 1 +ldouble: 1 +Test "expm1_tonearest (0x2.c5c4p+12)": +ildouble: 1 +ldouble: 1 +Test "expm1_tonearest (0xcp-4)": +double: 1 +idouble: 1 + +# expm1_towardzero +Test "expm1_towardzero (-0x1.86ap+16)": +ildouble: 1 +ldouble: 1 +Test "expm1_towardzero (-0x1p-100)": +double: 1 +float: 1 +idouble: 1 +ifloat: 1 +ildouble: 1 +ldouble: 1 +Test "expm1_towardzero (-0x1p-20)": +ildouble: 2 +ldouble: 2 +Test "expm1_towardzero (-0x1p-32)": +float: 1 +ifloat: 1 +ildouble: 1 +ldouble: 1 +Test "expm1_towardzero (-0x1p-64)": +double: 1 +float: 1 +idouble: 1 +ifloat: 1 +ildouble: 1 +ldouble: 1 +Test "expm1_towardzero (-0x2.71p+12)": +ildouble: 1 +ldouble: 1 +Test "expm1_towardzero (-0x2.dp+4)": +ildouble: 1 +ldouble: 1 +Test "expm1_towardzero (-0x3.e8p+8)": +ildouble: 1 +ldouble: 1 +Test "expm1_towardzero (-0x4.ap+4)": +ildouble: 1 +ldouble: 1 +Test "expm1_towardzero (-0x4.ep+4)": +ildouble: 1 +ldouble: 1 +Test "expm1_towardzero (-0x4.fp+4)": +ildouble: 1 +ldouble: 1 +Test "expm1_towardzero (-0x4p-12)": +ildouble: 1 +ldouble: 1 +Test "expm1_towardzero (-0x4p-52)": +float: 1 +ifloat: 1 +ildouble: 2 +ldouble: 2 +Test "expm1_towardzero (-0x8p-32)": +float: 1 +ifloat: 1 +ildouble: 1 +ldouble: 1 +Test "expm1_towardzero (-0xf.ffffffffffff8p+1020)": +ildouble: 1 +ldouble: 1 +Test "expm1_towardzero (-0xf.fffffffffffffffp+16380)": +ildouble: 1 +ldouble: 1 +Test "expm1_towardzero (-0xf.fffffp+124)": +ildouble: 1 +ldouble: 1 +Test "expm1_towardzero (0x1.f4p+8)": +double: 1 +idouble: 1 +Test "expm1_towardzero (0x1p+0)": +ildouble: 1 +ldouble: 1 +Test "expm1_towardzero (0x1p-100)": +ildouble: 1 +ldouble: 1 +Test "expm1_towardzero (0x1p-32)": +ildouble: 1 +ldouble: 1 +Test "expm1_towardzero (0x3.2p+4)": +double: 1 +float: 1 +idouble: 1 +ifloat: 1 +ildouble: 1 +ldouble: 1 +Test "expm1_towardzero (0x4p-12)": +ildouble: 1 +ldouble: 1 +Test "expm1_towardzero (0x4p-52)": +ildouble: 1 +ldouble: 1 +Test "expm1_towardzero (0x7.fp+4)": +double: 1 +idouble: 1 +ildouble: 1 +ldouble: 1 +Test "expm1_towardzero (0x8p-32)": +ildouble: 1 +ldouble: 1 + +# expm1_upward +Test "expm1_upward (-0x1.86ap+16)": +ildouble: 1 +ldouble: 1 +Test "expm1_upward (-0x1p-100)": +double: 1 +float: 1 +idouble: 1 +ifloat: 1 +ildouble: 1 +ldouble: 1 +Test "expm1_upward (-0x1p-20)": +ildouble: 2 +ldouble: 2 +Test "expm1_upward (-0x1p-32)": +float: 1 +ifloat: 1 +ildouble: 1 +ldouble: 1 +Test "expm1_upward (-0x1p-64)": +double: 1 +float: 1 +idouble: 1 +ifloat: 1 +ildouble: 1 +ldouble: 1 +Test "expm1_upward (-0x2.71p+12)": +ildouble: 1 +ldouble: 1 +Test "expm1_upward (-0x2.dp+4)": +ildouble: 1 +ldouble: 1 +Test "expm1_upward (-0x3.e8p+8)": +ildouble: 1 +ldouble: 1 +Test "expm1_upward (-0x4.ap+4)": +ildouble: 1 +ldouble: 1 +Test "expm1_upward (-0x4.ep+4)": +ildouble: 1 +ldouble: 1 +Test "expm1_upward (-0x4.fp+4)": +ildouble: 1 +ldouble: 1 +Test "expm1_upward (-0x4p-12)": +ildouble: 1 +ldouble: 1 +Test "expm1_upward (-0x4p-52)": +float: 1 +ifloat: 1 +ildouble: 2 +ldouble: 2 +Test "expm1_upward (-0x8p-32)": +float: 1 +ifloat: 1 +ildouble: 1 +ldouble: 1 +Test "expm1_upward (-0xf.ffffffffffff8p+1020)": +ildouble: 1 +ldouble: 1 +Test "expm1_upward (-0xf.fffffffffffffffp+16380)": +ildouble: 1 +ldouble: 1 +Test "expm1_upward (-0xf.fffffp+124)": +ildouble: 1 +ldouble: 1 +Test "expm1_upward (0x1.f4p+8)": +ildouble: 1 +ldouble: 1 +Test "expm1_upward (0x1p-100)": +double: 1 +float: 1 +idouble: 1 +ifloat: 1 +Test "expm1_upward (0x1p-32)": +float: 1 +ifloat: 1 +Test "expm1_upward (0x1p-64)": +double: 1 +float: 1 +idouble: 1 +ifloat: 1 +Test "expm1_upward (0x4p-4)": +ildouble: 1 +ldouble: 1 +Test "expm1_upward (0x4p-52)": +float: 1 +ifloat: 1 +Test "expm1_upward (0x8p-32)": +float: 1 +ifloat: 1 + # gamma Test "gamma (-0.5)": ildouble: 1 @@ -8864,6 +9145,9 @@ ldouble: 1 Test "sinh_downward (0x1.8p+4)": ildouble: 1 ldouble: 1 +Test "sinh_downward (0x8p-32)": +ildouble: 1 +ldouble: 1 Test "sinh_downward (22)": float: 1 ifloat: 1 @@ -8894,6 +9178,9 @@ ldouble: 1 Test "sinh_towardzero (0x1.8p+4)": ildouble: 1 ldouble: 1 +Test "sinh_towardzero (0x8p-32)": +ildouble: 1 +ldouble: 1 Test "sinh_towardzero (22)": float: 1 ifloat: 1 @@ -12014,6 +12301,38 @@ ifloat: 1 ildouble: 1 ldouble: 1 +Function: "expm1_downward": +double: 1 +float: 1 +idouble: 1 +ifloat: 1 +ildouble: 1 +ldouble: 1 + +Function: "expm1_tonearest": +double: 1 +float: 1 +idouble: 1 +ifloat: 1 +ildouble: 1 +ldouble: 1 + +Function: "expm1_towardzero": +double: 1 +float: 1 +idouble: 1 +ifloat: 1 +ildouble: 2 +ldouble: 2 + +Function: "expm1_upward": +double: 1 +float: 1 +idouble: 1 +ifloat: 1 +ildouble: 2 +ldouble: 2 + Function: "gamma": double: 1 float: 2 -- cgit 1.4.1