Fix numerous errors of >= 1 ulp for cosf(x) and sinf(x) (1 line)

and add a comment about related magic (many lines)).

__kernel_cos[f]() needs a trick to reduce the error to below 1 ulp
when |x| >= 0.3 for the range-reduced x.  Modulo other bugs, naive
code that doesn't use the trick would have an error of >= 1 ulp
in about 0.00006% of cases when |x| >= 0.3 for the unreduced x,
with a maximum relative error of about 1.03 ulps.  Mistransation
of the trick from the double precision case resulted in errors in
about 0.2% of cases, with a maximum relative error of about 1.3 ulps.

The mistranslation involved not doing implicit masking of the 32-bit
float word corresponding to to implicit masking of the lower 32-bit
double word by clearing it.

sinf() uses __kernel_cosf() for half of all cases so its errors from
this bug are similar.  tanf() is not affected.

The error bounds in the above and in my other recent commit messages
are for amd64.  Extra precision for floats on i386's accidentally masks
this bug, but only if k_cosf.c is compiled with -O.  Although the extra
precision helps here, this is accidental and depends on longstanding
gcc precision bugs (not clipping extra precision on assignment...),
and the gcc bugs are mostly avoided by compiling without -O.  I now
develop libm mainly on amd64 systems to simplify error detection and
debugging.

1 files changed, 12 insertions, 1 deletions

diff --git a/lib/msun/src/k_cosf.c b/lib/msun/src/k_cosf.c
index fb021ed..4239ae1 100644
--- a/lib/msun/src/k_cosf.c
+++ b/lib/msun/src/k_cosf.c
@@ -44,10 +44,21 @@ __kernel_cosf(float x, float y)
 	if(ix < 0x3e99999a) 			/* if |x| < 0.3 */
 	    return one - ((float)0.5*z - (z*r - x*y));
 	else {
+	    /*
+	     * qx is an approximation to x*x/2 such that 1-qx and x*x/2-qx
+	     * are both exact.  Its implementation is optimized for
+	     * efficiency in preference to accuracy.  We use x*x/2 ~ x/4 for
+	     * x near 0.5 and mask off just enough low bits (3) for both of
+	     * the above differences to be exact.  We use a constant for
+	     * x > 0.78125 to keep using the same algorithm as k_cos.c,
+	     * although this gives only a small improvement in accuracy, at
+	     * least here.  Using x*x/2 to approximate itself (masking off
+	     * 3 low bits again) would give better accuracy.
+	     */
 	    if(ix > 0x3f480000) {		/* x > 0.78125 */
 		qx = (float)0.28125;
 	    } else {
-	        SET_FLOAT_WORD(qx,ix-0x01000000);	/* x/4 */
+	        SET_FLOAT_WORD(qx,(ix-0x01000000)&0xfffffff8);
 	    }
 	    hz = (float)0.5*z-qx;
 	    a  = one-qx;