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.

4 files changed, 412 insertions, 0 deletions

diff --git a/sysdeps/i386/fpu/e_expl.S b/sysdeps/i386/fpu/e_expl.S
index af339f0..a8a5e70 100644
--- a/sysdeps/i386/fpu/e_expl.S
+++ b/sysdeps/i386/fpu/e_expl.S
@@ -130,9 +130,24 @@ 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.  */
+	subl	$8, %esp
+	cfi_adjust_cfa_offset (8)
+	fstcw	4(%esp)
+	movl	$0xf3ff, %edx
+	andl	4(%esp), %edx
+	movl	%edx, (%esp)
+	fldcw	(%esp)
+#endif
 	frndint			/* 1  i               */
 	fld	%st(1)		/* 2  x               */
 	frndint			/* 2  xi              */
+#ifdef USE_AS_EXPM1L
+	fldcw	4(%esp)
+	addl	$8, %esp
+	cfi_adjust_cfa_offset (-8)
+#endif
 	fld	%st(1)		/* 3  i               */
 	fldt	MO(c0)		/* 4  c0              */
 	fld	%st(2)		/* 5  xi              */
diff --git a/sysdeps/i386/fpu/libm-test-ulps b/sysdeps/i386/fpu/libm-test-ulps
index a471e5f..bb3adcb 100644
--- a/sysdeps/i386/fpu/libm-test-ulps
+++ b/sysdeps/i386/fpu/libm-test-ulps
@@ -6092,9 +6092,15 @@ idouble: 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
@@ -6108,6 +6114,338 @@ Test "expm1 (11356.25)":
 ildouble: 1
 ldouble: 1
 
+# expm1_downward
+Test "expm1_downward (-0x1p-100)":
+double: 1
+float: 1
+idouble: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+Test "expm1_downward (-0x2.ep+4)":
+idouble: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+Test "expm1_downward (-0x4.9p+4)":
+idouble: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+Test "expm1_downward (-0x4.bp+4)":
+idouble: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+Test "expm1_downward (-0x4p-4)":
+ildouble: 1
+ldouble: 1
+Test "expm1_downward (-0x5p+4)":
+idouble: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+Test "expm1_downward (-0x6.4p+4)":
+idouble: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+Test "expm1_downward (0x1p+0)":
+ildouble: 1
+ldouble: 1
+Test "expm1_downward (0x1p-100)":
+ildouble: 1
+ldouble: 1
+Test "expm1_downward (0x1p-32)":
+idouble: 1
+ildouble: 1
+ldouble: 1
+Test "expm1_downward (0x3.2p+4)":
+ildouble: 1
+ldouble: 1
+Test "expm1_downward (0x4p-12)":
+ildouble: 1
+ldouble: 1
+Test "expm1_downward (0x4p-52)":
+idouble: 1
+ildouble: 1
+ldouble: 1
+Test "expm1_downward (0x7.fp+4)":
+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 (0x1p+0)":
+ildouble: 1
+ldouble: 1
+Test "expm1_tonearest (0x2.c5c4p+12)":
+ildouble: 1
+ldouble: 1
+
+# expm1_towardzero
+Test "expm1_towardzero (-0x1.2p+4)":
+float: 1
+Test "expm1_towardzero (-0x1.86ap+16)":
+double: 1
+float: 1
+idouble: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+Test "expm1_towardzero (-0x1p-100)":
+double: 1
+float: 1
+ildouble: 1
+ldouble: 1
+Test "expm1_towardzero (-0x1p-20)":
+ildouble: 2
+ldouble: 2
+Test "expm1_towardzero (-0x1p-32)":
+idouble: 1
+ildouble: 1
+ldouble: 1
+Test "expm1_towardzero (-0x1p-64)":
+ildouble: 1
+ldouble: 1
+Test "expm1_towardzero (-0x2.4p+4)":
+float: 1
+Test "expm1_towardzero (-0x2.5p+4)":
+float: 1
+Test "expm1_towardzero (-0x2.6p+4)":
+double: 1
+float: 1
+Test "expm1_towardzero (-0x2.71p+12)":
+double: 1
+float: 1
+idouble: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+Test "expm1_towardzero (-0x2.cp+4)":
+double: 1
+float: 1
+Test "expm1_towardzero (-0x2.dp+4)":
+double: 1
+float: 1
+ildouble: 1
+ldouble: 1
+Test "expm1_towardzero (-0x2.ep+4)":
+double: 1
+float: 1
+Test "expm1_towardzero (-0x3.e8p+8)":
+double: 1
+float: 1
+idouble: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+Test "expm1_towardzero (-0x4.9p+4)":
+double: 1
+float: 1
+Test "expm1_towardzero (-0x4.ap+4)":
+double: 1
+float: 1
+ildouble: 1
+ldouble: 1
+Test "expm1_towardzero (-0x4.bp+4)":
+double: 1
+float: 1
+Test "expm1_towardzero (-0x4.ep+4)":
+double: 1
+float: 1
+ildouble: 1
+ldouble: 1
+Test "expm1_towardzero (-0x4.fp+4)":
+double: 1
+float: 1
+ildouble: 1
+ldouble: 1
+Test "expm1_towardzero (-0x4p-12)":
+ildouble: 1
+ldouble: 1
+Test "expm1_towardzero (-0x4p-52)":
+idouble: 1
+ildouble: 2
+ldouble: 2
+Test "expm1_towardzero (-0x5p+4)":
+double: 1
+float: 1
+Test "expm1_towardzero (-0x6.4p+4)":
+double: 1
+float: 1
+Test "expm1_towardzero (-0x8p-32)":
+ildouble: 1
+ldouble: 1
+Test "expm1_towardzero (-0xf.ffffffffffff8p+1020)":
+double: 1
+idouble: 1
+ildouble: 1
+ldouble: 1
+Test "expm1_towardzero (-0xf.fffffffffffffffp+16380)":
+ildouble: 1
+ldouble: 1
+Test "expm1_towardzero (-0xf.fffffp+124)":
+double: 1
+float: 1
+idouble: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+Test "expm1_towardzero (0x1p+0)":
+ildouble: 1
+ldouble: 1
+Test "expm1_towardzero (0x1p-100)":
+ildouble: 1
+ldouble: 1
+Test "expm1_towardzero (0x1p-32)":
+idouble: 1
+ildouble: 1
+ldouble: 1
+Test "expm1_towardzero (0x3.2p+4)":
+ildouble: 1
+ldouble: 1
+Test "expm1_towardzero (0x4p-12)":
+ildouble: 1
+ldouble: 1
+Test "expm1_towardzero (0x4p-52)":
+idouble: 1
+ildouble: 1
+ldouble: 1
+Test "expm1_towardzero (0x7.fp+4)":
+ildouble: 1
+ldouble: 1
+Test "expm1_towardzero (0x8p-32)":
+ildouble: 1
+ldouble: 1
+
+# expm1_upward
+Test "expm1_upward (-0x1.2p+4)":
+float: 1
+Test "expm1_upward (-0x1.86ap+16)":
+double: 1
+float: 1
+idouble: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+Test "expm1_upward (-0x1p-100)":
+double: 1
+float: 1
+ildouble: 1
+ldouble: 1
+Test "expm1_upward (-0x1p-20)":
+ildouble: 2
+ldouble: 2
+Test "expm1_upward (-0x1p-32)":
+idouble: 1
+ildouble: 1
+ldouble: 1
+Test "expm1_upward (-0x1p-64)":
+ildouble: 1
+ldouble: 1
+Test "expm1_upward (-0x2.4p+4)":
+float: 1
+Test "expm1_upward (-0x2.5p+4)":
+float: 1
+Test "expm1_upward (-0x2.6p+4)":
+double: 1
+float: 1
+Test "expm1_upward (-0x2.71p+12)":
+double: 1
+float: 1
+idouble: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+Test "expm1_upward (-0x2.cp+4)":
+double: 1
+float: 1
+Test "expm1_upward (-0x2.dp+4)":
+double: 1
+float: 1
+ildouble: 1
+ldouble: 1
+Test "expm1_upward (-0x2.ep+4)":
+double: 1
+float: 1
+Test "expm1_upward (-0x3.e8p+8)":
+double: 1
+float: 1
+idouble: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+Test "expm1_upward (-0x4.9p+4)":
+double: 1
+float: 1
+Test "expm1_upward (-0x4.ap+4)":
+double: 1
+float: 1
+ildouble: 1
+ldouble: 1
+Test "expm1_upward (-0x4.bp+4)":
+double: 1
+float: 1
+Test "expm1_upward (-0x4.ep+4)":
+double: 1
+float: 1
+ildouble: 1
+ldouble: 1
+Test "expm1_upward (-0x4.fp+4)":
+double: 1
+float: 1
+ildouble: 1
+ldouble: 1
+Test "expm1_upward (-0x4p-12)":
+ildouble: 1
+ldouble: 1
+Test "expm1_upward (-0x4p-52)":
+idouble: 1
+ildouble: 2
+ldouble: 2
+Test "expm1_upward (-0x5p+4)":
+double: 1
+float: 1
+Test "expm1_upward (-0x6.4p+4)":
+double: 1
+float: 1
+Test "expm1_upward (-0x8p-32)":
+ildouble: 1
+ldouble: 1
+Test "expm1_upward (-0xf.ffffffffffff8p+1020)":
+double: 1
+idouble: 1
+ildouble: 1
+ldouble: 1
+Test "expm1_upward (-0xf.fffffffffffffffp+16380)":
+ildouble: 1
+ldouble: 1
+Test "expm1_upward (-0xf.fffffp+124)":
+double: 1
+float: 1
+idouble: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+Test "expm1_upward (0x1.f4p+8)":
+ildouble: 1
+ldouble: 1
+Test "expm1_upward (0x4p-4)":
+ildouble: 1
+ldouble: 1
+
 # gamma
 Test "gamma (-0.5)":
 double: 1
@@ -7773,6 +8111,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 (0xcp-4)":
 float: 1
 ildouble: 1
@@ -7810,6 +8151,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 (0xcp-4)":
 float: 1
 ildouble: 1
@@ -7840,6 +8184,9 @@ ldouble: 1
 Test "sinh_upward (0x1.7p+4)":
 ildouble: 1
 ldouble: 1
+Test "sinh_upward (0x8p-32)":
+double: 1
+float: 1
 Test "sinh_upward (0xcp-4)":
 float: 1
 ildouble: 2
@@ -10663,6 +11010,34 @@ Function: "expm1":
 ildouble: 1
 ldouble: 1
 
+Function: "expm1_downward":
+double: 1
+float: 1
+idouble: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+
+Function: "expm1_tonearest":
+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
diff --git a/sysdeps/i386/fpu/s_expm1.S b/sysdeps/i386/fpu/s_expm1.S
index 91617eb..3980552 100644
--- a/sysdeps/i386/fpu/s_expm1.S
+++ b/sysdeps/i386/fpu/s_expm1.S
@@ -78,7 +78,18 @@ ENTRY(__expm1)
 5:	fldt	MO(l2e)		// log2(e) : x
 	fmulp			// log2(e)*x
 	fld	%st		// log2(e)*x : log2(e)*x
+	// Set round-to-nearest temporarily.
+	subl	$8, %esp
+	cfi_adjust_cfa_offset (8)
+	fstcw	4(%esp)
+	movl	$0xf3ff, %ecx
+	andl	4(%esp), %ecx
+	movl	%ecx, (%esp)
+	fldcw	(%esp)
 	frndint			// int(log2(e)*x) : log2(e)*x
+	fldcw	4(%esp)
+	addl	$8, %esp
+	cfi_adjust_cfa_offset (-8)
 	fsubr	%st, %st(1)	// int(log2(e)*x) : fract(log2(e)*x)
 	fxch			// fract(log2(e)*x) : int(log2(e)*x)
 	f2xm1			// 2^fract(log2(e)*x)-1 : int(log2(e)*x)
diff --git a/sysdeps/i386/fpu/s_expm1f.S b/sysdeps/i386/fpu/s_expm1f.S
index 2f1b6f6..e5e52a0 100644
--- a/sysdeps/i386/fpu/s_expm1f.S
+++ b/sysdeps/i386/fpu/s_expm1f.S
@@ -78,7 +78,18 @@ ENTRY(__expm1f)
 5:	fldt	MO(l2e)		// log2(e) : x
 	fmulp			// log2(e)*x
 	fld	%st		// log2(e)*x : log2(e)*x
+	// Set round-to-nearest temporarily.
+	subl	$8, %esp
+	cfi_adjust_cfa_offset (8)
+	fstcw	4(%esp)
+	movl	$0xf3ff, %ecx
+	andl	4(%esp), %ecx
+	movl	%ecx, (%esp)
+	fldcw	(%esp)
 	frndint			// int(log2(e)*x) : log2(e)*x
+	fldcw	4(%esp)
+	addl	$8, %esp
+	cfi_adjust_cfa_offset (-8)
 	fsubr	%st, %st(1)	// int(log2(e)*x) : fract(log2(e)*x)
 	fxch			// fract(log2(e)*x) : int(log2(e)*x)
 	f2xm1			// 2^fract(log2(e)*x)-1 : int(log2(e)*x)