Use only double precision for "kernel" cosf and sinf (except for

returning float).  The functions are renamed from __kernel_{cos,sin}f()
to __kernel_{cos,sin}df() so that misuses of them will cause link errors
and not crashes.

This version is an almost-routine translation with no special optimizations
for accuracy or efficiency.  The not-quite-routine part is that in
__kernel_cosf(), regenerating the minimax polynomial with double
precision coefficients gives a coefficient for the x**2 term that is
not quite -0.5, so the literal 0.5 in the code and the related `hz'
variable need to be modified; also, the special code for reducing the
error in 1.0-x**2*0.5 is no longer needed, so it is convenient to
adjust all the logic for the x**2 term a little.  Note that without
extra precision, it would be very bad to use a coefficient of other
than -0.5 for the x**2 term -- the old version depends on multiplication
by -0.5 being infinitely precise so as not to need even more special
code for reducing the error in 1-x**2*0.5.

This gives an unimportant increase in accuracy, from ~0.8 to ~0.501
ulps.  Almost all of the error is from the final rounding step, since
the choice of the minimax polynomials so that their contribution to the
error is a bit less than 0.5 ulps just happens to give contributions that
are significantly less (~.001 ulps).

An Athlons, for uniformly distributed args in [-2pi, 2pi], this gives
overall speed increases in the 10-20% range, despite giving a speed
decrease of typically 19% (from 31 cycles up to 37) for sinf() on args
in [-pi/4, pi/4].

1 files changed, 7 insertions, 19 deletions

diff --git a/lib/msun/src/s_sinf.c b/lib/msun/src/s_sinf.c
index 1e4270e..4986d9a 100644
--- a/lib/msun/src/s_sinf.c
+++ b/lib/msun/src/s_sinf.c
@@ -18,8 +18,8 @@ static char rcsid[] = "$FreeBSD$";
 #endif
 
 #include "math.h"
-#define	INLINE_KERNEL_COSF
-#define	INLINE_KERNEL_SINF
+#define	INLINE_KERNEL_COSDF
+#define	INLINE_KERNEL_SINDF
 #include "math_private.h"
 #include "k_cosf.c"
 #include "k_sinf.c"
@@ -31,18 +31,6 @@ s2pio2 = 2*M_PI_2,			/* 0x400921FB, 0x54442D18 */
 s3pio2 = 3*M_PI_2,			/* 0x4012D97C, 0x7F3321D2 */
 s4pio2 = 4*M_PI_2;			/* 0x401921FB, 0x54442D18 */
 
-static inline float
-__kernel_cosdf(double x)
-{
-	return __kernel_cosf((float)x, x - (float)x);
-}
-
-static inline float
-__kernel_sindf(double x)
-{
-	return __kernel_sinf((float)x, x - (float)x, 1);
-}
-
 float
 sinf(float x)
 {
@@ -55,7 +43,7 @@ sinf(float x)
 	if(ix <= 0x3f490fda) {		/* |x| ~<= pi/4 */
 	    if(ix<0x39800000)		/* |x| < 2**-12 */
 		if(((int)x)==0) return x;	/* x with inexact if x != 0 */
-	    return __kernel_sinf(x,0.0,0);
+	    return __kernel_sindf(x);
 	}
 	if(ix<=0x407b53d1) {		/* |x| <= ~5*pi/4 */
 	    if(ix<=0x4016cbe3) {	/* |x| <= ~3pi/4 */
@@ -83,11 +71,11 @@ sinf(float x)
 	else {
 	    n = __ieee754_rem_pio2f(x,y);
 	    switch(n&3) {
-		case 0: return  __kernel_sinf(y[0],y[1],1);
-		case 1: return  __kernel_cosf(y[0],y[1]);
-		case 2: return -__kernel_sinf(y[0],y[1],1);
+		case 0: return  __kernel_sindf((double)y[0]+y[1]);
+		case 1: return  __kernel_cosdf((double)y[0]+y[1]);
+		case 2: return -__kernel_sindf((double)y[0]+y[1]);
 		default:
-			return -__kernel_cosf(y[0],y[1]);
+			return -__kernel_cosdf((double)y[0]+y[1]);
 	    }
 	}
 }