As for the float trig functions and logf, use a minimax polynomial

that is specialized for float precision.  The new polynomial has degree
5 instead of 11, and a maximum error of 2**-27.74 ulps instead
of 2**-30.64.  This doesn't affect the final error significantly; the
maximum error was and is about 0.9101 ulps on amd64 -01 and the number
of cases with an error of > 0.5 ulps is actually reduced by epsilon
despite the larger error in the polynomial.

This is about 15% faster on amd64 (A64), i386 (A64) and ia64.  The asm
version is still used instead of this on i386 since it is faster and
more accurate.

1 files changed, 7 insertions, 6 deletions

diff --git a/lib/msun/src/e_expf.c b/lib/msun/src/e_expf.c
index 50b91da..949fa54 100644
--- a/lib/msun/src/e_expf.c
+++ b/lib/msun/src/e_expf.c
@@ -31,11 +31,12 @@ ln2HI[2]   ={ 6.9314575195e-01,		/* 0x3f317200 */
 ln2LO[2]   ={ 1.4286067653e-06,  	/* 0x35bfbe8e */
 	     -1.4286067653e-06,},	/* 0xb5bfbe8e */
 invln2 =  1.4426950216e+00, 		/* 0x3fb8aa3b */
-P1   =  1.6666667163e-01, /* 0x3e2aaaab */
-P2   = -2.7777778450e-03, /* 0xbb360b61 */
-P3   =  6.6137559770e-05, /* 0x388ab355 */
-P4   = -1.6533901999e-06, /* 0xb5ddea0e */
-P5   =  4.1381369442e-08; /* 0x3331bb4c */
+/*
+ * Domain [-0.34568, 0.34568], range ~[-4.278e-9, 4.447e-9]:
+ * |x*(exp(x)+1)/(exp(x)-1) - p(x)| < 2**-27.74
+ */
+P1 =  1.6666625440e-1,		/*  0xaaaa8f.0p-26 */
+P2 = -2.7667332906e-3;		/* -0xb55215.0p-32 */
 
 static volatile float twom100 = 7.8886090522e-31;      /* 2**-100=0x0d800000 */
 
@@ -79,7 +80,7 @@ __ieee754_expf(float x)	/* default IEEE double exp */
 
     /* x is now in primary range */
 	t  = x*x;
-	c  = x - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
+	c  = x - t*(P1+t*P2);
 	if(k==0) 	return one-((x*c)/(c-(float)2.0)-x);
 	else 		y = one-((lo-(x*c)/((float)2.0-c))-hi);
 	if(k >= -125) {