Improve the accuracy from a max ULP of ~2000 to max ULP < 0.79

on i386-class hardware for sinl and cosl.  The hand-rolled argument
reduction have been replaced by e_rem_pio2l() implementations.  To
preserve history the following commands have been executed:

svn cp src/e_rem_pio2.c ld80/e_rem_pio2l.h
mv ${HOME}/bde/ld80/e_rem_pio2l.c ld80/e_rem_pio2l.h

svn cp src/e_rem_pio2.c ld128/e_rem_pio2l.h
mv ${HOME}/bde/ld128/e_rem_pio2l.c ld128/e_rem_pio2l.h

The ld80 version has been tested by bde, das, and kargl over the
last few years (bde, das) and few months (kargl).  An older ld128
version was tested by das.  The committed version has only been
compiled tested via 'make universe'.

Approved by: das (mentor)
Obtained from: bde

5 files changed, 321 insertions, 98 deletions

diff --git a/lib/msun/ld128/e_rem_pio2l.h b/lib/msun/ld128/e_rem_pio2l.h
new file mode 100644
index 0000000..d223627
--- /dev/null
+++ b/lib/msun/ld128/e_rem_pio2l.h
@@ -0,0 +1,143 @@
+/* From: @(#)e_rem_pio2.c 1.4 95/01/18 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ * Copyright (c) 2008 Steven G. Kargl, David Schultz, Bruce D. Evans.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ *
+ * Optimized by Bruce D. Evans.
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/* ld128 version of __ieee754_rem_pio2l(x,y)
+ * 
+ * return the remainder of x rem pi/2 in y[0]+y[1] 
+ * use __kernel_rem_pio2()
+ */
+
+#include <float.h>
+
+#include "math.h"
+#include "math_private.h"
+#include "fpmath.h"
+
+#define	BIAS	(LDBL_MAX_EXP - 1)
+
+/*
+ * XXX need to verify that nonzero integer multiples of pi/2 within the
+ * range get no closer to a long double than 2**-140, or that
+ * ilogb(x) + ilogb(min_delta) < 45 - -140.
+ */
+/*
+ * invpio2:  113 bits of 2/pi
+ * pio2_1:   first  68 bits of pi/2
+ * pio2_1t:  pi/2 - pio2_1
+ * pio2_2:   second 68 bits of pi/2
+ * pio2_2t:  pi/2 - (pio2_1+pio2_2)
+ * pio2_3:   third  68 bits of pi/2
+ * pio2_3t:  pi/2 - (pio2_1+pio2_2+pio2_3)
+ */
+
+static const double
+zero =  0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
+two24 =  1.67772160000000000000e+07; /* 0x41700000, 0x00000000 */
+
+static const long double
+invpio2 =  6.3661977236758134307553505349005747e-01L,	/*  0x145f306dc9c882a53f84eafa3ea6a.0p-113 */
+pio2_1  =  1.5707963267948966192292994253909555e+00L,	/*  0x1921fb54442d18469800000000000.0p-112 */
+pio2_1t =  2.0222662487959507323996846200947577e-21L,	/*  0x13198a2e03707344a4093822299f3.0p-181 */
+pio2_2  =  2.0222662487959507323994779168837751e-21L,	/*  0x13198a2e03707344a400000000000.0p-181 */
+pio2_2t =  2.0670321098263988236496903051604844e-43L,	/*  0x127044533e63a0105df531d89cd91.0p-254 */
+pio2_3  =  2.0670321098263988236499468110329591e-43L,	/*  0x127044533e63a0105e00000000000.0p-254 */
+pio2_3t = -2.5650587247459238361625433492959285e-65L;	/* -0x159c4ec64ddaeb5f78671cbfb2210.0p-327 */
+
+#ifdef INLINE_REM_PIO2L
+static inline __always_inline
+#endif
+int
+__ieee754_rem_pio2l(long double x, long double *y)
+{
+	union IEEEl2bits u,u1;
+	long double z,w,t,r,fn;
+	double tx[5],ty[3];
+	int64_t n;
+	int e0,ex,i,j,nx;
+	int16_t expsign;
+
+	u.e = x;
+	expsign = u.xbits.expsign;
+	ex = expsign & 0x7fff;
+	if (ex < BIAS + 45 || ex == BIAS + 45 &&
+	    u.bits.manh < 0x921fb54442d1LL) {
+	    /* |x| ~< 2^45*(pi/2), medium size */
+	    /* Use a specialized rint() to get fn.  Assume round-to-nearest. */
+	    fn = x*invpio2+0x1.8p112;
+	    fn = fn-0x1.8p112;
+#ifdef HAVE_EFFICIENT_I64RINT
+	    n  = i64rint(fn);
+#else
+	    n  = fn;
+#endif
+	    r  = x-fn*pio2_1;
+	    w  = fn*pio2_1t;	/* 1st round good to 180 bit */
+	    {
+		union IEEEl2bits u2;
+	        int ex1;
+	        j  = ex;
+	        y[0] = r-w; 
+		u2.e = y[0];
+		ex1 = u2.xbits.expsign & 0x7fff;
+	        i = j-ex1;
+	        if(i>51) {  /* 2nd iteration needed, good to 248 */
+		    t  = r;
+		    w  = fn*pio2_2;	
+		    r  = t-w;
+		    w  = fn*pio2_2t-((t-r)-w);	
+		    y[0] = r-w;
+		    u2.e = y[0];
+		    ex1 = u2.xbits.expsign & 0x7fff;
+		    i = j-ex1;
+		    if(i>119) {	/* 3rd iteration need, 316 bits acc */
+		    	t  = r;	/* will cover all possible cases */
+		    	w  = fn*pio2_3;	
+		    	r  = t-w;
+		    	w  = fn*pio2_3t-((t-r)-w);	
+		    	y[0] = r-w;
+		    }
+		}
+	    }
+	    y[1] = (r-y[0])-w;
+	    return n;
+	}
+    /* 
+     * all other (large) arguments
+     */
+	if(ex==0x7fff) {		/* x is inf or NaN */
+	    y[0]=y[1]=x-x; return 0;
+	}
+    /* set z = scalbn(|x|,ilogb(x)-23) */
+	u1.e = x;
+	e0 = ex - BIAS - 23;		/* e0 = ilogb(|x|)-23; */
+	u1.xbits.expsign = ex - e0;
+	z = u1.e;
+	for(i=0;i<4;i++) {
+		tx[i] = (double)((int32_t)(z));
+		z     = (z-tx[i])*two24;
+	}
+	tx[4] = z;
+	nx = 5;
+	while(tx[nx-1]==zero) nx--;	/* skip zero term */
+	n  =  __kernel_rem_pio2(tx,ty,e0,nx,3);
+	t = (long double)ty[2] + ty[1];
+	r = t + ty[0];
+	w = ty[0] - (r - t);
+	if(expsign<0) {y[0] = -r; y[1] = -w; return -n;}
+	y[0] = r; y[1] = w; return n;
+}
diff --git a/lib/msun/ld80/e_rem_pio2l.h b/lib/msun/ld80/e_rem_pio2l.h
new file mode 100644
index 0000000..cfe3fb9
--- /dev/null
+++ b/lib/msun/ld80/e_rem_pio2l.h
@@ -0,0 +1,152 @@
+/* From: @(#)e_rem_pio2.c 1.4 95/01/18 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ * Copyright (c) 2008 Steven G. Kargl, David Schultz, Bruce D. Evans.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ *
+ * Optimized by Bruce D. Evans.
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/* ld80 version of __ieee754_rem_pio2l(x,y)
+ * 
+ * return the remainder of x rem pi/2 in y[0]+y[1] 
+ * use __kernel_rem_pio2()
+ */
+
+#include <float.h>
+
+#include "math.h"
+#include "math_private.h"
+#include "fpmath.h"
+
+#define	BIAS	(LDBL_MAX_EXP - 1)
+
+/*
+ * invpio2:  64 bits of 2/pi
+ * pio2_1:   first  39 bits of pi/2
+ * pio2_1t:  pi/2 - pio2_1
+ * pio2_2:   second 39 bits of pi/2
+ * pio2_2t:  pi/2 - (pio2_1+pio2_2)
+ * pio2_3:   third  39 bits of pi/2
+ * pio2_3t:  pi/2 - (pio2_1+pio2_2+pio2_3)
+ */
+
+static const double
+zero =  0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
+two24 =  1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
+pio2_1  =  1.57079632679597125389e+00,	/* 0x3FF921FB, 0x54444000 */
+pio2_2  = -1.07463465549783099519e-12,	/* -0x12e7b967674000.0p-92 */
+pio2_3  =  6.36831716351370313614e-25;	/*  0x18a2e037074000.0p-133 */
+
+#if defined(__amd64__) || defined(__i386__)
+/* Long double constants are slow on these arches, and broken on i386. */
+static const volatile double
+invpio2hi =  6.3661977236758138e-01,	/*  0x145f306dc9c883.0p-53 */
+invpio2lo = -3.9356538861223811e-17,	/* -0x16b00000000000.0p-107 */
+pio2_1thi = -1.0746346554971943e-12,	/* -0x12e7b9676733af.0p-92 */
+pio2_1tlo =  8.8451028997905949e-29,	/*  0x1c080000000000.0p-146 */
+pio2_2thi =  6.3683171635109499e-25,	/*  0x18a2e03707344a.0p-133 */
+pio2_2tlo =  2.3183081793789774e-41,	/*  0x10280000000000.0p-187 */
+pio2_3thi = -2.7529965190440717e-37,	/* -0x176b7ed8fbbacc.0p-174 */
+pio2_3tlo = -4.2006647512740502e-54;	/* -0x19c00000000000.0p-230 */
+#define	invpio2	((long double)invpio2hi + invpio2lo)
+#define	pio2_1t	((long double)pio2_1thi + pio2_1tlo)
+#define	pio2_2t	((long double)pio2_2thi + pio2_2tlo)
+#define	pio2_3t	((long double)pio2_3thi + pio2_3tlo)
+#else
+static const long double
+invpio2 =  6.36619772367581343076e-01L,	/*  0xa2f9836e4e44152a.0p-64 */
+pio2_1t = -1.07463465549719416346e-12L,	/* -0x973dcb3b399d747f.0p-103 */
+pio2_2t =  6.36831716351095013979e-25L,	/*  0xc51701b839a25205.0p-144 */
+pio2_3t = -2.75299651904407171810e-37L;	/* -0xbb5bf6c7ddd660ce.0p-185 */
+#endif
+
+#ifdef INLINE_REM_PIO2L
+static inline __always_inline
+#endif
+int
+__ieee754_rem_pio2l(long double x, long double *y)
+{
+	union IEEEl2bits u,u1;
+	long double z,w,t,r,fn;
+	double tx[3],ty[2];
+	int e0,ex,i,j,nx,n;
+	int16_t expsign;
+
+	u.e = x;
+	expsign = u.xbits.expsign;
+	ex = expsign & 0x7fff;
+	if (ex < BIAS + 25 || ex == BIAS + 25 && u.bits.manh < 0xc90fdaa2) {
+	    /* |x| ~< 2^25*(pi/2), medium size */
+	    /* Use a specialized rint() to get fn.  Assume round-to-nearest. */
+	    fn = x*invpio2+0x1.8p63;
+	    fn = fn-0x1.8p63;
+#ifdef HAVE_EFFICIENT_IRINT
+	    n  = irint(fn);
+#else
+	    n  = fn;
+#endif
+	    r  = x-fn*pio2_1;
+	    w  = fn*pio2_1t;	/* 1st round good to 102 bit */
+	    {
+		union IEEEl2bits u2;
+	        int ex1;
+	        j  = ex;
+	        y[0] = r-w; 
+		u2.e = y[0];
+		ex1 = u2.xbits.expsign & 0x7fff;
+	        i = j-ex1;
+	        if(i>22) {  /* 2nd iteration needed, good to 141 */
+		    t  = r;
+		    w  = fn*pio2_2;	
+		    r  = t-w;
+		    w  = fn*pio2_2t-((t-r)-w);	
+		    y[0] = r-w;
+		    u2.e = y[0];
+		    ex1 = u2.xbits.expsign & 0x7fff;
+		    i = j-ex1;
+		    if(i>61) {	/* 3rd iteration need, 180 bits acc */
+		    	t  = r;	/* will cover all possible cases */
+		    	w  = fn*pio2_3;	
+		    	r  = t-w;
+		    	w  = fn*pio2_3t-((t-r)-w);	
+		    	y[0] = r-w;
+		    }
+		}
+	    }
+	    y[1] = (r-y[0])-w;
+	    return n;
+	}
+    /* 
+     * all other (large) arguments
+     */
+	if(ex==0x7fff) {		/* x is inf or NaN */
+	    y[0]=y[1]=x-x; return 0;
+	}
+    /* set z = scalbn(|x|,ilogb(x)-23) */
+	u1.e = x;
+	e0 = ex - BIAS - 23;		/* e0 = ilogb(|x|)-23; */
+	u1.xbits.expsign = ex - e0;
+	z = u1.e;
+	for(i=0;i<2;i++) {
+		tx[i] = (double)((int32_t)(z));
+		z     = (z-tx[i])*two24;
+	}
+	tx[2] = z;
+	nx = 3;
+	while(tx[nx-1]==zero) nx--;	/* skip zero term */
+	n  =  __kernel_rem_pio2(tx,ty,e0,nx,2);
+	r = (long double)ty[0] + ty[1];
+	w = ty[1] - (r - ty[0]);
+	if(expsign<0) {y[0] = -r; y[1] = -w; return -n;}
+	y[0] = r; y[1] = w; return n;
+}
diff --git a/lib/msun/src/s_cosl.c b/lib/msun/src/s_cosl.c
index d184aa9..77611c1 100644
--- a/lib/msun/src/s_cosl.c
+++ b/lib/msun/src/s_cosl.c
@@ -28,7 +28,6 @@
 __FBSDID("$FreeBSD$");
 
 /*
- * Compute cos(x) for x where x is reduced to y = x - k * pi / 2.
  * Limited testing on pseudorandom numbers drawn within [-2e8:4e8] shows
  * an accuracy of <= 0.7412 ULP.
  */
@@ -36,27 +35,22 @@ __FBSDID("$FreeBSD$");
 #include <float.h>
 
 #include "math.h"
+#define INLINE_REM_PIO2L
 #include "math_private.h"
-#include "fpmath.h"
-
 #if LDBL_MANT_DIG == 64
-#define	NX	3
-#define	PREC	2
+#include "../ld80/e_rem_pio2l.h"
 #elif LDBL_MANT_DIG == 113
-#define	NX	5
-#define	PREC	3
+#include "../ld128/e_rem_pio2l.h"
 #else
 #error "Unsupported long double format"
 #endif
 
-static const long double two24 = 1.67772160000000000000e+07L;
-
 long double
 cosl(long double x)
 {
 	union IEEEl2bits z;
-	int i, e0;
-	double xd[NX], yd[PREC];
+	int e0;
+	long double y[2];
 	long double hi, lo;
 
 	z.e = x;
@@ -74,26 +68,9 @@ cosl(long double x)
 	if (z.e < M_PI_4)
 		return (__kernel_cosl(z.e, 0));
 
-	/* Split z.e into a 24-bit representation. */
-	e0 = ilogbl(z.e) - 23;
-	z.e = scalbnl(z.e, -e0);
-	for (i = 0; i < NX; i++) {
-		xd[i] = (double)((int32_t)z.e);
-		z.e = (z.e - xd[i]) * two24;
-	}
-	
-	/* yd contains the pieces of xd rem pi/2 such that |yd| < pi/4. */
-	e0 = __kernel_rem_pio2(xd, yd, e0, NX, PREC);
-
-#if PREC == 2
-	hi = (long double)yd[0] + yd[1];
-	lo = yd[1] - (hi - yd[0]);
-#else /* PREC == 3 */
-	long double t;
-	t = (long double)yd[2] + yd[1];
-	hi = t + yd[0];
-	lo = yd[0] - (hi - t);
-#endif
+	e0 = __ieee754_rem_pio2l(x, y);
+	hi = y[0];
+	lo = y[1];
 
 	switch (e0 & 3) {
 	case 0:
diff --git a/lib/msun/src/s_sinl.c b/lib/msun/src/s_sinl.c
index 65d7026..19d9572 100644
--- a/lib/msun/src/s_sinl.c
+++ b/lib/msun/src/s_sinl.c
@@ -27,34 +27,25 @@
 #include <sys/cdefs.h>
 __FBSDID("$FreeBSD$");
 
-/*
- * Compute sin(x) for x where x is reduced to y = x - k * pi / 2.
- */
-
 #include <float.h>
 
 #include "math.h"
+#define INLINE_REM_PIO2L
 #include "math_private.h"
-#include "fpmath.h"
-
 #if LDBL_MANT_DIG == 64
-#define	NX	3
-#define	PREC	2
+#include "../ld80/e_rem_pio2l.h"
 #elif LDBL_MANT_DIG == 113
-#define	NX	5
-#define	PREC	3
+#include "../ld128/e_rem_pio2l.h"
 #else
 #error "Unsupported long double format"
 #endif
 
-static const long double two24 = 1.67772160000000000000e+07L;
-
 long double
 sinl(long double x)
 {
 	union IEEEl2bits z;
-	int i, e0, s;
-	double xd[NX], yd[PREC];
+	int e0, s;
+	long double y[2];
 	long double hi, lo;
 
 	z.e = x;
@@ -75,26 +66,9 @@ sinl(long double x)
 		return  (s ? -hi : hi);
 	}
 
-	/* Split z.e into a 24-bit representation. */
-	e0 = ilogbl(z.e) - 23;
-	z.e = scalbnl(z.e, -e0);
-	for (i = 0; i < NX; i++) {
-		xd[i] = (double)((int32_t)z.e);
-		z.e = (z.e - xd[i]) * two24;
-	}
-
-	/* yd contains the pieces of xd rem pi/2 such that |yd| < pi/4. */
-	e0 = __kernel_rem_pio2(xd, yd, e0, NX, PREC);
-
-#if PREC == 2
-	hi = (long double)yd[0] + yd[1];
-	lo = yd[1] - (hi - yd[0]);
-#else /* PREC == 3 */
-	long double t;
-	t = (long double)yd[2] + yd[1];
-	hi = t + yd[0];
-	lo = yd[0] - (hi - t);
-#endif
+	e0 = __ieee754_rem_pio2l(x, y);
+	hi = y[0];
+	lo = y[1];
 
 	switch (e0 & 3) {
 	case 0:
@@ -111,5 +85,5 @@ sinl(long double x)
 	    break;
 	}
 	
-	return (s ? -hi : hi);
+	return (hi);
 }
diff --git a/lib/msun/src/s_tanl.c b/lib/msun/src/s_tanl.c
index 2592dcd..3cf4a69 100644
--- a/lib/msun/src/s_tanl.c
+++ b/lib/msun/src/s_tanl.c
@@ -28,7 +28,6 @@
 __FBSDID("$FreeBSD$");
 
 /*
- * Compute tan(x) for x where x is reduced to y = x - k * pi / 2.
  * Limited testing on pseudorandom numbers drawn within [0:4e8] shows
  * an accuracy of <= 1.5 ULP where 247024 values of x out of 40 million
  * possibles resulted in tan(x) that exceeded 0.5 ULP (ie., 0.6%).
@@ -37,27 +36,22 @@ __FBSDID("$FreeBSD$");
 #include <float.h>
 
 #include "math.h"
+#define INLINE_REM_PIO2L
 #include "math_private.h"
-#include "fpmath.h"
-
 #if LDBL_MANT_DIG == 64
-#define	NX	3
-#define	PREC	2
+#include "../ld80/e_rem_pio2l.h"
 #elif LDBL_MANT_DIG == 113
-#define	NX	5
-#define	PREC	3
+#include "../ld128/e_rem_pio2l.h"
 #else
 #error "Unsupported long double format"
 #endif
 
-static const long double two24 = 1.67772160000000000000e+07L;
-
 long double
 tanl(long double x)
 {
 	union IEEEl2bits z;
-	int i, e0, s;
-	double xd[NX], yd[PREC];
+	int e0, s;
+	long double y[2];
 	long double hi, lo;
 
 	z.e = x;
@@ -78,26 +72,9 @@ tanl(long double x)
 		return (s ? -hi : hi);
 	}
 
-	/* Split z.e into a 24-bit representation. */
-	e0 = ilogbl(z.e) - 23;
-	z.e = scalbnl(z.e, -e0);
-	for (i = 0; i < NX; i++) {
-		xd[i] = (double)((int32_t)z.e);
-		z.e = (z.e - xd[i]) * two24;
-	}
-	
-	/* yd contains the pieces of xd rem pi/2 such that |yd| < pi/4. */
-	e0 = __kernel_rem_pio2(xd, yd, e0, NX, PREC);
-
-#if PREC == 2
-	hi = (long double)yd[0] + yd[1];
-	lo = yd[1] - (hi - yd[0]);
-#else /* PREC == 3 */
-	long double t;
-	t = (long double)yd[2] + yd[1];
-	hi = t + yd[0];
-	lo = yd[0] - (hi - t);
-#endif
+	e0 = __ieee754_rem_pio2l(x, y);
+	hi = y[0];
+	lo = y[1];
 
 	switch (e0 & 3) {
 	case 0:
@@ -110,5 +87,5 @@ tanl(long double x)
 	    break;
 	}
 
-	return (s ? -hi : hi);
+	return (hi);
 }