/* C/80 3.0 Transcendental function library (6/6/83).  (c) 1983 Walt Bilofsky;
   (c) 1982 Knowledge Engineering Inc.	Reproduction of all or any part of
   this material in source or assembly language form is prohibited.
   Object programs containing executable code generated by this material
   may be distributed without restriction; credit to C/80 is appreciated. */

/* Source for the transcendental functions: "Computer Approximations,"
   by Hart, Cheney, et al., John Wiley & Sons, Inc., NY, 1968.	Part of
   the SIAM series in Applied Mathematics.  */

#define pi2    1.570796326795
#define pi     3.14159265359
#define sqrt10 3.162277660168
#define ln10   2.302585092994
#define log10e 0.43429448190
#define log2   0.30102999567

#ifneed cos
float cos(f) float f;
{   float sin(), fabs();
    return sin(fabs(f) + pi2); }
#endif

#ifneed sin
float sin(arg)
float arg;
{	/* Coefficients are #3370 from Hart & Cheney (18.80D). */
	static float p[5] = {
	  .1357884e8, -.49429081e7, .440103053e6, -.138472724e5, .145968841e3 };
	static float q[5] = {
	  .864455865e7, .408179225e6, .946309610e4, .132653491e3, 1.0 };
	static float y,ysq;
	static long quad,e;
	float F_poly();

	quad = (arg >= 0.0) ? 0 : (arg = -arg, 2);
	arg *= .63661977;		/* 2 / pi */
	y = arg - (e = arg);
	quad=(e + quad) % 4;
	if (quad & 1) y = 1.0 - y;
	if(quad & 2) y = -y;

	ysq = y * y;
	return y * F_poly(ysq,p,4)/F_poly(ysq,q,4);
}
#endif

#ifneed atan
float atan(f) float f;
/*	Eqns. 6.5.21, & 6.5.22	& Table 5097 */
{	static float p[6] = {
		33.05861847399, 58.655751569, 32.3909748562,
		5.853195211263, .1952374193623, -.002434603300441 };
	static float q[5] = {
		33.05861847399, 69.67529105952, 49.00434821822,
		12.97557886271, 1.0};
	static int sign, recip;
	float F_poly(),fabs();

	float f2, f2sq;
	if ((f2 = fabs(f)) == 0.0) return 0;
	if (recip = (f2 > 1.0)) f2 = 1.0 / f2;	/* 0 < f2 < 1 */
	f2sq = f2 * f2;
	f2 =* F_poly(f2sq, p, 5) / F_poly(f2sq, q, 4);
	if (recip) f2 = pi2 - f2;
	return f < 0.0 ? -f2 : f2;
}
#endif

#ifneed sqrt
float sqrt(f) float f;
/*	Eqns. 6.1.3 & 6.1.7; Initial guess from table 0231 */
{	static float p[] = {.58812e-2, .5267875, .5881229 };
	static float f2, temp, scale;
	static int t;
	float F_poly();

	if (f < 0.0) return -1.0;
	if (f == 0.0) return 0.0;
	scale = 1.0;
	while (f > 1.0) { f =* 0.01; scale =* 10.0; }
	while (f < 0.01) { f =* 100.0; scale =* 0.1; }
	f2 = F_poly(f, p, 2) / (f + .999998e-1);
	for (t = 10; t--; ) {
		temp = 0.5 * (f2 + f / f2);
		if (temp == f2) break;
		f2 = temp; }
	return f2 * scale;
}
#endif

#ifneed exp
float exp(f) float f;		/* = e^f */
{	float pow10();
	return pow10(log10e * f); }
#endif

#ifneed pow
float pow(f1,f2)	      /* f1 ^ f2 */
float f1,f2;
{	static int i;
	static float p,ans;
	float log(),fabs(),pow10();

	p = fabs(f2);		/* Note - doesn't check for neg power */
	ans = 1.0;
	if ((i = p) < 10 && (float)i == p) while (i--) ans =* f1;
	    else ans = pow10(p * log(f1));
	return f2 < 0.0 ?  1.0/ans : ans;
	}
#endif

#ifneed pow10
float pow10(f) float f; 	/* 10^f */
/*	Eq. 6.2.34 & Table 1403 */
{	static float c[7] = {
		1.0, .2302585092158E1, .2650949191501E1, .2034670312104E1,
		.1171493469046E1, .5358845480519, .2321569786604 };
	static float powtn[9] = {
		.1258925411794E+01, .1584893192461E+01, .1995262314969E+01,
		.2511886431510E+01, .3162277660168E+01, .3981071705535E+01,
		.5011872336273E+01, .6309573444802E+01, .7943282347243E+01 };
	static float powt[9] = {1e1,1e2,1e3,1e4,1e5,1e6,1e7,1e8,1e9};
	static float powtt[3] = {1e10, 1e20, 1e30};
	static int n1, n2, recip;
	static float f2;
	float F_poly();

	if (recip = (f < 0.0)) f = -f;
	if (f > 39.0) return 0.0;	/* Error - overflow */
	n1 = n2 = 0;
	if (f >= 0.1) {
		n1 = f; f =- n1;
		n2 = 10.0 * f; f =- 0.1 * n2; }
	f2 = f ? F_poly(f,c,6) : 1.0;
	if (n2) f2 =* powtn[n2-1];
	if (n1 >= 10) {	f2 =* powtt[n1/10 - 1];
			n1 =% 10; }
	if (n1) f2 =* powt[n1-1];
	return recip ? 1.0/f2 : f2; }
#endif

#ifneed ln
float ln(f) float f; { float log(); return log(f) * ln10; }
#endif

/* Natural and base 10 log routines */
#ifneed log
float log(f) float f;
/*	From Table 2325 & Eq. 6.3.28 */
{	static float p[4] = { -8.625170319686, 11.70031513942,
			      -3.932918863942, 0.1960631750250};
	static float q[4] = { -9.930094301066, 16.78051701465, -8.135425915213,
			       1.0 };
	static float n, z, zsq, num;
	static int pow2;
	float F_poly();
	union { float f_;
		char c_[4]; } ff;


	if (f <= 0) return 0;		/* Error */
	ff.f_ = f;
	pow2 = (0xFF & ff.c_[3]) - 0x80;/* Get exponent and set to 0 */
	ff.c_[3] = 0x80; f = ff.f_;	/* and put back in f. */
	if (f == 0.5) { z = 0.0; --pow2;}	/* Log 1 is 0 */
	else {	num = f * sqrt10;	/* Now compute log10 of fract. part */
		z = (num - 1.0) / (num + 1.0);
		zsq = z * z;
		z = (z * F_poly(zsq,p,3) / F_poly(zsq,q,3) - 0.5);
		}
	return z + pow2 * log2; 	/* Add exponent back in */
	}
#endif

#ifneed F_poly
float F_poly(x, p, n)  /* evaluate nth deg. polynomial, n >= 1 */
float x,p[];
{   static float result;  static int i;
    for (result = p[i = n]; i; ) result = p[--i] + x * result;
    return result; }
#endif

#ifneed fabs
float fabs(f) float f; { return f < 0 ? -f : f; }
#endif
; ) re