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, 60 insertions, 12 deletions

diff --git a/math/libm-test.inc b/math/libm-test.inc
index dea6c8b7ce..aab3ed2987 100644
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -267,18 +267,6 @@ struct ulp_data
 #define TEST_COND_before_rounding	(!TININESS_AFTER_ROUNDING)
 #define TEST_COND_after_rounding	TININESS_AFTER_ROUNDING
 
-#ifdef __x86_64__
-# define TEST_COND_x86_64	1
-#else
-# define TEST_COND_x86_64	0
-#endif
-
-#ifdef __i386__
-# define TEST_COND_x86		1
-#else
-# define TEST_COND_x86		0
-#endif
-
 /* Various constants (we must supply them precalculated for accuracy).  */
 #define M_PI_6l			.52359877559829887307710723054658383L
 #define M_PI_34l		2.356194490192344928846982537459627163L	/* 3*pi/4 */
@@ -7845,6 +7833,62 @@ expm1_test (void)
 }
 
 
+static const struct test_f_f_data expm1_tonearest_test_data[] =
+  {
+    AUTO_TESTS_f_f (expm1, tonearest),
+  };
+
+static void
+expm1_test_tonearest (void)
+{
+  START (expm1_tonearest);
+  RUN_TEST_LOOP_f_f (expm1, expm1_tonearest_test_data, FE_TONEAREST);
+  END;
+}
+
+
+static const struct test_f_f_data expm1_towardzero_test_data[] =
+  {
+    AUTO_TESTS_f_f (expm1, towardzero),
+  };
+
+static void
+expm1_test_towardzero (void)
+{
+  START (expm1_towardzero);
+  RUN_TEST_LOOP_f_f (expm1, expm1_towardzero_test_data, FE_TOWARDZERO);
+  END;
+}
+
+
+static const struct test_f_f_data expm1_downward_test_data[] =
+  {
+    AUTO_TESTS_f_f (expm1, downward),
+  };
+
+static void
+expm1_test_downward (void)
+{
+  START (expm1_downward);
+  RUN_TEST_LOOP_f_f (expm1, expm1_downward_test_data, FE_DOWNWARD);
+  END;
+}
+
+
+static const struct test_f_f_data expm1_upward_test_data[] =
+  {
+    AUTO_TESTS_f_f (expm1, upward),
+  };
+
+static void
+expm1_test_upward (void)
+{
+  START (expm1_upward);
+  RUN_TEST_LOOP_f_f (expm1, expm1_upward_test_data, FE_UPWARD);
+  END;
+}
+
+
 static const struct test_f_f_data fabs_test_data[] =
   {
     TEST_f_f (fabs, 0, 0, NO_INEXACT_EXCEPTION),
@@ -13337,6 +13381,10 @@ main (int argc, char **argv)
   exp10_test ();
   exp2_test ();
   expm1_test ();
+  expm1_test_tonearest ();
+  expm1_test_towardzero ();
+  expm1_test_downward ();
+  expm1_test_upward ();
   frexp_test ();
   ldexp_test ();
   log_test ();