[PATCH] math/libm-test.inc: Fix cos and sincos tests against zero.

["Carlos O'Donell" <carlos at redhat dot com>]

Andreas, Joseph,

This is a precursor to my cleanup of the definition of ulp in the math
test suite. This patch brings two tests to 1/2 ulp for x86-64 and x86, 
otherwise when I turn on precise ulp calculations these tests show 
up as errors. The problem here is independent of the definition of ulp.

The value of PI is never exactly PI in any floating point representation,
and the value of PI/2 is never PI/2. It is wrong to expect cos(M_PI_2l)
to return 0, instead it will return an answer that is  non-zero because
M_PI_2l doesn't round to 0 in the type used.

That is to say that the correct answer is to do the following:
* Take PI or PI/2.
* Round to the floating point representation.
* Take the rounded value and compute an infinite precision cos or sin.
* Use the rounded result of the infinite precision cos or sin as the answer.

I used printf to do the type rounding e.g.
printf ("%e, %.100e, %a\n", (float)(2.0L*M_PIl), (float)(2.0L*M_PIl), (float)(2.0L*M_PIl));
printf ("%e, %.100e, %a\n", (double)(2.0L*M_PIl), (double)(2.0L*M_PIl), (double)(2.0L*M_PIl));
printf ("%Le, %.100Le, %La\n", (long double)(2.0L*M_PIl), (long double)(2.0L*M_PIl), (long double)(2.0L*M_PIl));
printf ("%e, %.100e, %a\n", (float)(M_El), (float)(M_El), (float)(M_El));
printf ("%e, %.100e, %a\n", (double)(M_El), (double)(M_El), (double)(M_El));
printf ("%Le, %.100Le, %La\n", (long double)(M_El), (long double)(M_El), (long double)(M_El));

Then I used Wolfram's Alpha here to get an infinite precision sin/cos.
http://www.wolframalpha.com/

The following changes bring x86-64 and x86 to 1/2 ulp for two tests.
It shows that the x86 cos implementation is quite good, and that
our test are flawed.

Unfortunately given that the rounding errors are type dependent we
need to fix this for each type. No regressions on x86-64 or x86.

Is this the right way forward, even if it is logically the correct
answer? I would say yes. I would say the test suite as it is right
now is testing against the "naive" expected result. It works because
ulp(0) for the testsuite is a huge value that allows all sorts of
inaccurate zero values.

Comments?

2013-04-09  Carlos O'Donell <carlos@redhat.com>

	* math/libm-test.inc: 
	* sysdeps/x86_64/fpu/libm-test-ulps: Regenerate.
	* sysdeps/i386/fpu/libm-test-ulps: Regenerate.

diff --git a/math/libm-test.inc b/math/libm-test.inc
index 78d2107..de42102 100644
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -5330,7 +5330,20 @@ cos_test (void)
 
   TEST_f_f (cos, M_PI_6l * 2.0, 0.5);
   TEST_f_f (cos, M_PI_6l * 4.0, -0.5);
-  TEST_f_f (cos, M_PI_2l, 0);
+
+  /* The value of M_PI_2l is never exactly PI/2, and therefore the
+     answer is never exactly zero. The answer is equal to the error
+     in rounding PI/2 for the type used.  Thus the answer is unique
+     to each type.  */
+#ifdef TEST_FLOAT
+  TEST_f_f (cos, M_PI_2l, -4.371139000186241438857289400265215e-8L);
+#endif
+#if defined TEST_DOUBLE
+  TEST_f_f (cos, M_PI_2l, 6.123233995736765886130329661375001e-17L);
+#endif
+#if defined TEST_LDOUBLE
+  TEST_f_f (cos, M_PI_2l, -2.50827880633416601177866354016537e-20L);
+#endif
 
   TEST_f_f (cos, 0.75L, 0.731688868873820886311838753000084544L);
 
@@ -12134,7 +12146,20 @@ sincos_test (void)
   TEST_extra (sincos, minus_infty, qnan_value, qnan_value, INVALID_EXCEPTION);
   TEST_extra (sincos, qnan_value, qnan_value, qnan_value);
 
-  TEST_extra (sincos, M_PI_2l, 1, 0);
+  /* The value of M_PI_2l is never exactly PI/2, and therefore the
+     answer is never exactly zero. The answer is equal to the error
+     in rounding PI/2 for the type used.  Thus the answer is unique
+     to each type.  */
+#ifdef TEST_FLOAT
+  TEST_extra (sincos, M_PI_2l, 1, -4.371139000186241438857289400265215e-8L);
+#endif
+#if defined TEST_DOUBLE
+  TEST_extra (sincos, M_PI_2l, 1, 6.123233995736765886130329661375001e-17L);
+#endif
+#if defined TEST_LDOUBLE
+  TEST_extra (sincos, M_PI_2l, 1, -2.50827880633416601177866354016537e-20L);
+#endif
+
   TEST_extra (sincos, M_PI_6l, 0.5, 0.86602540378443864676372317075293616L);
   TEST_extra (sincos, M_PI_6l*2.0, 0.86602540378443864676372317075293616L, 0.5);
   TEST_extra (sincos, 0.75L, 0.681638760023334166733241952779893935L, 0.731688868873820886311838753000084544L);
diff --git a/sysdeps/x86_64/fpu/libm-test-ulps b/sysdeps/x86_64/fpu/libm-test-ulps
index 28f4bfc..9c9c473 100644
--- a/sysdeps/x86_64/fpu/libm-test-ulps
+++ b/sysdeps/x86_64/fpu/libm-test-ulps
@@ -4717,13 +4717,6 @@ idouble: 2
 ifloat: 1
 ildouble: 1
 ldouble: 1
-Test "cos (pi/2) == 0":
-double: 1
-float: 1
-idouble: 1
-ifloat: 1
-ildouble: 1
-ldouble: 1
 
 # cos_downward
 Test "cos_downward (1) == 0.5403023058681397174009366074429766037323":
@@ -5984,13 +5977,6 @@ idouble: 1
 ifloat: 1
 ildouble: 1
 ldouble: 1
-Test "sincos (pi/2, &sin_res, &cos_res) puts 0 in cos_res":
-double: 1
-float: 1
-idouble: 1
-ifloat: 1
-ildouble: 1
-ldouble: 1
 Test "sincos (pi/6, &sin_res, &cos_res) puts 0.86602540378443864676372317075293616 in cos_res":
 float: 1
 ifloat: 1
diff --git a/sysdeps/i386/fpu/libm-test-ulps b/sysdeps/i386/fpu/libm-test-ulps
index 4de455c..fecaa94 100644
--- a/sysdeps/i386/fpu/libm-test-ulps
+++ b/sysdeps/i386/fpu/libm-test-ulps
@@ -4149,13 +4149,6 @@ idouble: 2
 ifloat: 1
 ildouble: 1
 ldouble: 1
-Test "cos (pi/2) == 0":
-double: 1
-float: 1
-idouble: 1
-ifloat: 1
-ildouble: 1
-ldouble: 1
 
 # cos_downward
 Test "cos_downward (1) == 0.5403023058681397174009366074429766037323":
@@ -5513,13 +5506,6 @@ idouble: 1
 ifloat: 1
 ildouble: 1
 ldouble: 1
-Test "sincos (pi/2, &sin_res, &cos_res) puts 0 in cos_res":
-double: 1
-float: 1
-idouble: 1
-ifloat: 1
-ildouble: 1
-ldouble: 1
 
 # sinh
 Test "sinh (0.75) == 0.822316731935829980703661634446913849":
---

Cheers,
Carlos.