From bd06cb56ab7f07dd1d6e23daa8b00a6b75f2b37c Mon Sep 17 00:00:00 2001 From: bde Date: Wed, 6 Feb 2008 06:35:21 +0000 Subject: 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. --- lib/msun/src/e_expf.c | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) (limited to 'lib') 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) { -- cgit v1.1