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.

2 files changed, 330 insertions, 0 deletions

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