(*  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.  See the copyright notice in the ACK home directory, in the file "Copyright".*)(*$R-*)IMPLEMENTATION MODULE Mathlib;(*  Module:	Mathematical functions  Author:	Ceriel J.H. Jacobs  Version:	$Header: /cvsup/minix/src/lib/ack/libm2/Mathlib.mod,v 1.1 2005/10/10 15:27:46 beng Exp $*)  FROM	EM IMPORT	FIF, FEF;  FROM	Traps IMPORT	Message;  CONST	OneRadianInDegrees	= 57.295779513082320876798155D;	OneDegreeInRadians	=  0.017453292519943295769237D;	OneOverSqrt2		= 0.70710678118654752440084436210484904D;  (* basic functions *)  PROCEDURE pow(x: REAL; i: INTEGER): REAL;  BEGIN	RETURN SHORT(longpow(LONG(x), i));  END pow;  PROCEDURE longpow(x: LONGREAL; i: INTEGER): LONGREAL;    VAR	val: LONGREAL;	ri: LONGREAL;  BEGIN	ri := FLOATD(i);	IF x < 0.0D THEN		val := longexp(longln(-x) * ri);		IF ODD(i) THEN RETURN -val;		ELSE RETURN val;		END;	ELSIF x = 0.0D THEN		RETURN 0.0D;	ELSE		RETURN longexp(longln(x) * ri);	END;  END longpow;  PROCEDURE sqrt(x: REAL): REAL;  BEGIN	RETURN SHORT(longsqrt(LONG(x)));  END sqrt;  PROCEDURE longsqrt(x: LONGREAL): LONGREAL;    VAR	temp: LONGREAL;	exp, i: INTEGER;  BEGIN	IF x <= 0.0D THEN		IF x < 0.0D THEN			Message("sqrt: negative argument");			HALT		END;		RETURN 0.0D;	END;	temp := FEF(x,exp);	(*	 * NOTE	 * this wont work on 1's comp	 *)	IF ODD(exp) THEN		temp := 2.0D * temp;		DEC(exp);	END;	temp := 0.5D*(1.0D + temp);	WHILE exp > 28 DO		temp := temp * 16384.0D;		exp := exp - 28;	END;	WHILE exp < -28 DO		temp := temp / 16384.0D;		exp := exp + 28;	END;	WHILE exp >= 2 DO		temp := temp * 2.0D;		exp := exp - 2;	END;	WHILE exp <= -2 DO		temp := temp / 2.0D;		exp := exp + 2;	END;	FOR i := 0 TO 5 DO		temp := 0.5D*(temp + x/temp);	END;	RETURN temp;  END longsqrt;  PROCEDURE ldexp(x:LONGREAL; n: INTEGER): LONGREAL;  BEGIN	WHILE n >= 16 DO		x := x * 65536.0D;		n := n - 16;	END;	WHILE n > 0 DO		x := x * 2.0D;		DEC(n);	END;	WHILE n <= -16 DO		x := x / 65536.0D;		n := n + 16;	END;	WHILE n < 0 DO		x := x / 2.0D;		INC(n);	END;	RETURN x;  END ldexp;  PROCEDURE exp(x: REAL): REAL;  BEGIN	RETURN SHORT(longexp(LONG(x)));  END exp;  PROCEDURE longexp(x: LONGREAL): LONGREAL;  (*	Algorithm and coefficients from:		"Software manual for the elementary functions"		by W.J. Cody and W. Waite, Prentice-Hall, 1980  *)    CONST	p0 = 0.25000000000000000000D+00;	p1 = 0.75753180159422776666D-02;	p2 = 0.31555192765684646356D-04;	q0 = 0.50000000000000000000D+00;	q1 = 0.56817302698551221787D-01;	q2 = 0.63121894374398503557D-03;	q3 = 0.75104028399870046114D-06;    VAR	neg: BOOLEAN;	n: INTEGER;	xn, g, x1, x2: LONGREAL;  BEGIN	neg := x < 0.0D;	IF neg THEN		x := -x;	END;	n := TRUNC(x/longln2 + 0.5D);	xn := FLOATD(n);	x1 := FLOATD(TRUNCD(x));	x2 := x - x1;	g := ((x1 - xn * 0.693359375D)+x2) - xn * (-2.1219444005469058277D-4);	IF neg THEN		g := -g;		n := -n;	END;	xn := g*g;	x := g*((p2*xn+p1)*xn+p0);	INC(n);	RETURN ldexp(0.5D + x/((((q3*xn+q2)*xn+q1)*xn+q0) - x), n);  END longexp;  PROCEDURE ln(x: REAL): REAL;	(* natural log *)  BEGIN	RETURN SHORT(longln(LONG(x)));  END ln;  PROCEDURE longln(x: LONGREAL): LONGREAL;	(* natural log *)  (*	Algorithm and coefficients from:		"Software manual for the elementary functions"		by W.J. Cody and W. Waite, Prentice-Hall, 1980   *)    CONST	p0 = -0.64124943423745581147D+02;	p1 =  0.16383943563021534222D+02;	p2 = -0.78956112887491257267D+00;	q0 = -0.76949932108494879777D+03;	q1 =  0.31203222091924532844D+03;	q2 = -0.35667977739034646171D+02;	q3 =  1.0D;    VAR	exp: INTEGER;	z, znum, zden, w: LONGREAL;  BEGIN	IF x <= 0.0D THEN		Message("ln: argument <= 0");		HALT	END;	x := FEF(x, exp);	IF x > OneOverSqrt2 THEN		znum := (x - 0.5D) - 0.5D;		zden := x * 0.5D + 0.5D;	ELSE		znum := x - 0.5D;		zden := znum * 0.5D + 0.5D;		DEC(exp);	END;	z := znum / zden;	w := z * z;	x := z + z * w * (((p2*w+p1)*w+p0)/(((q3*w+q2)*w+q1)*w+q0));	z := FLOATD(exp);	x := x + z * (-2.121944400546905827679D-4);	RETURN x + z * 0.693359375D;  END longln;  PROCEDURE log(x: REAL): REAL;	(* log with base 10 *)  BEGIN	RETURN SHORT(longlog(LONG(x)));  END log;  PROCEDURE longlog(x: LONGREAL): LONGREAL;	(* log with base 10 *)  BEGIN	RETURN longln(x)/longln10;  END longlog;  (* trigonometric functions; arguments in radians *)  PROCEDURE sin(x: REAL): REAL;  BEGIN	RETURN SHORT(longsin(LONG(x)));  END sin;  PROCEDURE sinus(x: LONGREAL; cosflag: BOOLEAN) : LONGREAL;  (*	Algorithm and coefficients from:		"Software manual for the elementary functions"		by W.J. Cody and W. Waite, Prentice-Hall, 1980  *)    CONST	r0 = -0.16666666666666665052D+00;	r1 =  0.83333333333331650314D-02;	r2 = -0.19841269841201840457D-03;	r3 =  0.27557319210152756119D-05;	r4 = -0.25052106798274584544D-07;	r5 =  0.16058936490371589114D-09;	r6 = -0.76429178068910467734D-12;	r7 =  0.27204790957888846175D-14;	A1 =  3.1416015625D;	A2 = -8.908910206761537356617D-6;    VAR	x1, x2, y : LONGREAL;	neg : BOOLEAN;  BEGIN	IF x < 0.0D THEN		neg := TRUE;		x := -x	ELSE	neg := FALSE	END;	IF cosflag THEN		neg := FALSE;		y := longhalfpi + x	ELSE		y := x	END;	y := y / longpi + 0.5D;	IF FIF(y, 1.0D, y) < 0.0D THEN ; END;	IF FIF(y, 0.5D, x1) # 0.0D THEN neg := NOT neg END;	IF cosflag THEN y := y - 0.5D END;	x2 := FIF(x, 1.0, x1);	x := x1 - y * A1;	x := x + x2;	x := x - y * A2;	IF x < 0.0D THEN		neg := NOT neg;		x := -x	END;	y := x * x;	x := x + x * y * (((((((r7*y+r6)*y+r5)*y+r4)*y+r3)*y+r2)*y+r1)*y+r0);	IF neg THEN RETURN -x END;	RETURN x;  END sinus;  PROCEDURE longsin(x: LONGREAL): LONGREAL;  BEGIN	RETURN sinus(x, FALSE);  END longsin;  PROCEDURE cos(x: REAL): REAL;  BEGIN	RETURN SHORT(longcos(LONG(x)));  END cos;  PROCEDURE longcos(x: LONGREAL): LONGREAL;  BEGIN	IF x < 0.0D THEN x := -x; END;	RETURN sinus(x, TRUE);	  END longcos;  PROCEDURE tan(x: REAL): REAL;  BEGIN	RETURN SHORT(longtan(LONG(x)));  END tan;  PROCEDURE longtan(x: LONGREAL): LONGREAL;  (*	Algorithm and coefficients from:		"Software manual for the elementary functions"		by W.J. Cody and W. Waite, Prentice-Hall, 1980  *)    CONST	p1 = -0.13338350006421960681D+00;	p2 =  0.34248878235890589960D-02;	p3 = -0.17861707342254426711D-04;	q0 =  1.0D;	q1 = -0.46671683339755294240D+00;	q2 =  0.25663832289440112864D-01;	q3 = -0.31181531907010027307D-03;	q4 =  0.49819433993786512270D-06;	A1 =  1.57080078125D;	A2 = -4.454455103380768678308D-06;    VAR y, x1, x2: LONGREAL;	negative: BOOLEAN;	invert: BOOLEAN;  BEGIN	negative := x < 0.0D;	y := x / longhalfpi + 0.5D;        (*      Use extended precision to calculate reduced argument.                Here we used 12 bits of the mantissa for a1.                Also split x in integer part x1 and fraction part x2.        *)	IF FIF(y, 1.0D, y) < 0.0D THEN ; END;	invert := FIF(y, 0.5D, x1) # 0.0D;	x2 := FIF(x, 1.0D, x1);	x := x1 - y * A1;	x := x + x2;	x := x - y * A2;	y := x * x;	x := x + x * y * ((p3*y+p2)*y+p1);	y := (((q4*y+q3)*y+q2)*y+q1)*y+q0;	IF negative THEN x := -x END;	IF invert THEN RETURN -y/x END;	RETURN x/y;  END longtan;  PROCEDURE arcsin(x: REAL): REAL;  BEGIN	RETURN SHORT(longarcsin(LONG(x)));  END arcsin;  PROCEDURE arcsincos(x: LONGREAL; cosfl: BOOLEAN): LONGREAL;    CONST	p0 = -0.27368494524164255994D+02;	p1 =  0.57208227877891731407D+02;	p2 = -0.39688862997540877339D+02;	p3 =  0.10152522233806463645D+02;	p4 = -0.69674573447350646411D+00;	q0 = -0.16421096714498560795D+03;	q1 =  0.41714430248260412556D+03;	q2 = -0.38186303361750149284D+03;	q3 =  0.15095270841030604719D+03;	q4 = -0.23823859153670238830D+02;	q5 =  1.0D;    VAR	negative : BOOLEAN;	big: BOOLEAN;	g: LONGREAL;  BEGIN	negative := x < 0.0D;	IF negative THEN x := -x; END;	IF x > 0.5D THEN		big := TRUE;		IF x > 1.0D THEN			Message("arcsin or arccos: argument > 1");			HALT		END;		g := 0.5D - 0.5D * x;		x := -longsqrt(g);		x := x + x;	ELSE		big := FALSE;		g := x * x;	END;	x := x + x * g *	  ((((p4*g+p3)*g+p2)*g+p1)*g+p0)/(((((q5*g+q4)*g+q3)*g+q2)*g+q1)*g+q0);	IF cosfl AND NOT negative THEN x := -x END;	IF cosfl = NOT big THEN		x := (x + longquartpi) + longquartpi;	ELSIF cosfl AND negative AND big THEN		x := (x + longhalfpi) + longhalfpi;	END;	IF negative AND NOT cosfl THEN x := -x END;	RETURN x;  END arcsincos;	  PROCEDURE longarcsin(x: LONGREAL): LONGREAL;  BEGIN	RETURN arcsincos(x, FALSE);  END longarcsin;  PROCEDURE arccos(x: REAL): REAL;  BEGIN	RETURN SHORT(longarccos(LONG(x)));  END arccos;  PROCEDURE longarccos(x: LONGREAL): LONGREAL;  BEGIN	RETURN arcsincos(x, TRUE);  END longarccos;  PROCEDURE arctan(x: REAL): REAL;  BEGIN	RETURN SHORT(longarctan(LONG(x)));  END arctan;  VAR A: ARRAY[0..3] OF LONGREAL;      arctaninit: BOOLEAN;  PROCEDURE longarctan(x: LONGREAL): LONGREAL;  (*	Algorithm and coefficients from:		"Software manual for the elementary functions"		by W.J. Cody and W. Waite, Prentice-Hall, 1980  *)    CONST	p0 = -0.13688768894191926929D+02;	p1 = -0.20505855195861651981D+02;	p2 = -0.84946240351320683534D+01;	p3 = -0.83758299368150059274D+00;	q0 =  0.41066306682575781263D+02;	q1 =  0.86157349597130242515D+02;	q2 =  0.59578436142597344465D+02;	q3 =  0.15024001160028576121D+02;	q4 =  1.0D;    VAR	g: LONGREAL;	neg: BOOLEAN;	n: INTEGER;  BEGIN	IF NOT arctaninit THEN		arctaninit := TRUE;		A[0] := 0.0D;		A[1] := 0.52359877559829887307710723554658381D;	(* p1/6 *)		A[2] := longhalfpi;		A[3] := 1.04719755119659774615421446109316763D; (* pi/3 *)	END;	neg := FALSE;	IF x < 0.0D THEN		neg := TRUE;		x := -x;	END;	IF x > 1.0D THEN		x := 1.0D/x;		n := 2	ELSE		n := 0	END;	IF x > 0.26794919243112270647D (* 2-sqrt(3) *) THEN		INC(n);		x := (((0.73205080756887729353D*x-0.5D)-0.5D)+x)/			(1.73205080756887729353D + x);	END;	g := x*x;	x := x + x * g * (((p3*g+p2)*g+p1)*g+p0) / ((((q4*g+q3)*g+q2)*g+q1)*g+q0);	IF n > 1 THEN x := -x END;	x := x + A[n];	IF neg THEN RETURN -x; END;	RETURN x;  END longarctan;  (* hyperbolic functions *)  (* The C math library has better implementations for some of these, but     they depend on some properties of the floating point implementation,     and, for now, we don't want that in the Modula-2 system.  *)  PROCEDURE sinh(x: REAL): REAL;  BEGIN	RETURN SHORT(longsinh(LONG(x)));  END sinh;  PROCEDURE longsinh(x: LONGREAL): LONGREAL;    VAR expx: LONGREAL;  BEGIN	expx := longexp(x);	RETURN (expx - 1.0D/expx)/2.0D;  END longsinh;  PROCEDURE cosh(x: REAL): REAL;  BEGIN	RETURN SHORT(longcosh(LONG(x)));  END cosh;  PROCEDURE longcosh(x: LONGREAL): LONGREAL;    VAR expx: LONGREAL;  BEGIN	expx := longexp(x);	RETURN (expx + 1.0D/expx)/2.0D;  END longcosh;  PROCEDURE tanh(x: REAL): REAL;  BEGIN	RETURN SHORT(longtanh(LONG(x)));  END tanh;  PROCEDURE longtanh(x: LONGREAL): LONGREAL;    VAR expx: LONGREAL;  BEGIN	expx := longexp(x);	RETURN (expx - 1.0D/expx) / (expx + 1.0D/expx);  END longtanh;  PROCEDURE arcsinh(x: REAL): REAL;  BEGIN	RETURN SHORT(longarcsinh(LONG(x)));  END arcsinh;  PROCEDURE longarcsinh(x: LONGREAL): LONGREAL;    VAR neg: BOOLEAN;  BEGIN	neg := FALSE;	IF x < 0.0D THEN		neg := TRUE;		x := -x;	END;	x := longln(x + longsqrt(x*x+1.0D));	IF neg THEN RETURN -x; END;	RETURN x;  END longarcsinh;  PROCEDURE arccosh(x: REAL): REAL;  BEGIN	RETURN SHORT(longarccosh(LONG(x)));  END arccosh;  PROCEDURE longarccosh(x: LONGREAL): LONGREAL;  BEGIN	IF x < 1.0D THEN		Message("arccosh: argument < 1");		HALT	END;	RETURN longln(x + longsqrt(x*x - 1.0D));  END longarccosh;  PROCEDURE arctanh(x: REAL): REAL;  BEGIN	RETURN SHORT(longarctanh(LONG(x)));  END arctanh;  PROCEDURE longarctanh(x: LONGREAL): LONGREAL;  BEGIN	IF (x <= -1.0D) OR (x >= 1.0D) THEN		Message("arctanh: ABS(argument) >= 1");		HALT	END;	RETURN longln((1.0D + x)/(1.0D - x)) / 2.0D;  END longarctanh;  (* conversions *)  PROCEDURE RadianToDegree(x: REAL): REAL;  BEGIN	RETURN SHORT(longRadianToDegree(LONG(x)));  END RadianToDegree;  PROCEDURE longRadianToDegree(x: LONGREAL): LONGREAL;  BEGIN	RETURN x * OneRadianInDegrees;  END longRadianToDegree;  PROCEDURE DegreeToRadian(x: REAL): REAL;  BEGIN	RETURN SHORT(longDegreeToRadian(LONG(x)));  END DegreeToRadian;  PROCEDURE longDegreeToRadian(x: LONGREAL): LONGREAL;  BEGIN	RETURN x * OneDegreeInRadians;  END longDegreeToRadian;BEGIN	arctaninit := FALSE;END Mathlib.