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.

1 files changed, 11 insertions, 0 deletions

diff --git a/sysdeps/i386/fpu/s_expm1.S b/sysdeps/i386/fpu/s_expm1.S
index 91617eb74c..398055276e 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)