minix/lib/libm2/Mathlib.mod
2005-04-21 14:53:53 +00:00

576 lines
12 KiB
Modula-2
Executable file

(*
(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$
*)
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.