use new math routines

lang/basic/lib/exp.c

@@ -10,19 +10,6 @@
 #define __NO_DEFS
 #include <math.h>
 
-static double
-floor(x)
-	double x;
-{
-	extern double _fif();
-	double val;
-
-	return _fif(x, 1.0, &val) < 0 ? val - 1.0 : val ;
-	/*	this also works if _fif always returns a positive
-		fractional part
-	*/
-}
-
 static double
 ldexp(fl,exp)
 	double fl;
@@ -57,52 +44,54 @@ ldexp(fl,exp)
 
 double
 _exp(x)
-	double x;
+	double	x;
 {
-	/*	2**x = (Q(x*x)+x*P(x*x))/(Q(x*x)-x*P(x*x)) for x in [0,0.5] */
-	/*	Hart & Cheney #1069 */
+	/*	Algorithm and coefficients from:
+			"Software manual for the elementary functions"
+			by W.J. Cody and W. Waite, Prentice-Hall, 1980
+	*/
 
-	static double p[3] = {
-		 0.2080384346694663001443843411e+07,
-		 0.3028697169744036299076048876e+05,
-		 0.6061485330061080841615584556e+02
+	static double p[] = {
+	        0.25000000000000000000e+0,
+	        0.75753180159422776666e-2,
+		0.31555192765684646356e-4
 	};
 
-	static double q[4] = {
-		 0.6002720360238832528230907598e+07,
-		 0.3277251518082914423057964422e+06,
-		 0.1749287689093076403844945335e+04,
-		 0.1000000000000000000000000000e+01
+	static double q[] = {
+	        0.50000000000000000000e+0,
+	        0.56817302698551221787e-1,
+	        0.63121894374398503557e-3,
+		0.75104028399870046114e-6
 	};
+	double	xn, g;
+	int	n;
+	int	negative = x < 0;
 
-	int negative = x < 0;
-	int ipart, large = 0;
-	double xsqr, xPxx, Qxx;
-
-	if (x < M_LN_MIN_D) {
+	if (x <= M_LN_MIN_D) {
 		return M_MIN_D;
 	}
 	if (x >= M_LN_MAX_D) {
 		if (x > M_LN_MAX_D) error(3);
 		return M_MAX_D;
 	}
+	if (negative) x = -x;
 
+	/* ??? avoid underflow ??? */
+
+	n = x * M_LOG2E + 0.5;	/* 1/ln(2) = log2(e), 0.5 added for rounding */
+	xn = n;
+	{
+		double	x1 = (long) x;
+		double	x2 = x - x1;
+
+		g = ((x1-xn*0.693359375)+x2) - xn*(-2.1219444005469058277e-4);
+	}
 	if (negative) {
-		x = -x;
+		g = -g;
+		n = -n;
 	}
-	x /= M_LN2;
-	ipart = floor(x);
-	x -= ipart;
-	if (x > 0.5) {
-		large = 1;
-		x -= 0.5;
-	}
-	xsqr = x * x;
-	xPxx = x * POLYNOM2(xsqr, p);
-	Qxx = POLYNOM3(xsqr, q);
-	x = (Qxx + xPxx) / (Qxx - xPxx);
-	if (large) x *= M_SQRT2;
-	x = ldexp(x, ipart);
-	if (negative) return 1.0/x;
-	return x;
+	xn = g * g;
+	x = g * POLYNOM2(xn, p);
+	n += 1;
+	return (ldexp(0.5 + x/(POLYNOM3(xn, q) - x), n));
 }