455 lines
9 KiB
Modula-2
455 lines
9 KiB
Modula-2
(*
|
|
(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;
|
|
|
|
(* From: Handbook of Mathematical Functions
|
|
Edited by M. Abramowitz and I.A. Stegun
|
|
National Bureau of Standards
|
|
Applied Mathematics Series 55
|
|
*)
|
|
|
|
CONST
|
|
OneRadianInDegrees = 57.295779513082320876798155D;
|
|
OneDegreeInRadians = 0.017453292519943295769237D;
|
|
|
|
(* 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 4 DO
|
|
temp := 0.5D*(temp + x/temp);
|
|
END;
|
|
RETURN temp;
|
|
END longsqrt;
|
|
|
|
PROCEDURE exp(x: REAL): REAL;
|
|
BEGIN
|
|
RETURN SHORT(longexp(LONG(x)));
|
|
END exp;
|
|
|
|
PROCEDURE longexp(x: LONGREAL): LONGREAL;
|
|
(*
|
|
* n = floor(x / ln2), d = x / ln2 - n
|
|
* exp(x) = exp((x / ln2) * ln2) = exp((n + d) * ln2) =
|
|
* exp(n * ln2) * exp(d * ln2) = 2 ** n * exp(d * ln2)
|
|
*)
|
|
CONST
|
|
a1 = -0.9999999995D;
|
|
a2 = 0.4999999206D;
|
|
a3 = -0.1666653019D;
|
|
a4 = 0.0416573475D;
|
|
a5 = -0.0083013598D;
|
|
a6 = 0.0013298820D;
|
|
a7 = -0.0001413161D;
|
|
VAR
|
|
neg: BOOLEAN;
|
|
polval: LONGREAL;
|
|
n: LONGREAL;
|
|
n1 : INTEGER;
|
|
BEGIN
|
|
neg := x < 0.0D;
|
|
IF neg THEN
|
|
x := -x;
|
|
END;
|
|
x := FIF(x, 1.0D/LONG(ln2), n) * LONG(ln2);
|
|
polval := 1.0D /(1.0D + x*(a1+x*(a2+x*(a3+x*(a4+x*(a5+x*(a6+x*a7)))))));
|
|
n1 := TRUNCD(n + 0.5D);
|
|
WHILE n1 >= 16 DO
|
|
polval := polval * 65536.0D;
|
|
n1 := n1 - 16;
|
|
END;
|
|
WHILE n1 > 0 DO
|
|
polval := polval * 2.0D;
|
|
DEC(n1);
|
|
END;
|
|
IF neg THEN RETURN 1.0D/polval; END;
|
|
RETURN polval;
|
|
END longexp;
|
|
|
|
PROCEDURE ln(x: REAL): REAL; (* natural log *)
|
|
BEGIN
|
|
RETURN SHORT(longln(LONG(x)));
|
|
END ln;
|
|
|
|
PROCEDURE longln(x: LONGREAL): LONGREAL; (* natural log *)
|
|
CONST
|
|
a1 = 0.9999964239D;
|
|
a2 = -0.4998741238D;
|
|
a3 = 0.3317990258D;
|
|
a4 = -0.2407338084D;
|
|
a5 = 0.1676540711D;
|
|
a6 = -0.0953293897D;
|
|
a7 = 0.0360884937D;
|
|
a8 = -0.0064535442D;
|
|
VAR
|
|
exp: INTEGER;
|
|
polval: LONGREAL;
|
|
|
|
BEGIN
|
|
IF x <= 0.0D THEN
|
|
Message("ln: argument <= 0");
|
|
HALT
|
|
END;
|
|
x := FEF(x, exp);
|
|
WHILE x < 1.0D DO
|
|
x := x + x;
|
|
DEC(exp);
|
|
END;
|
|
x := x - 1.0D;
|
|
polval := x*(a1+x*(a2+x*(a3+x*(a4+x*(a5+x*(a6+x*(a7+a8*x)))))));
|
|
RETURN polval + FLOATD(exp) * LONG(ln2);
|
|
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)/LONG(ln10);
|
|
END longlog;
|
|
|
|
(* trigonometric functions; arguments in radians *)
|
|
|
|
PROCEDURE sin(x: REAL): REAL;
|
|
BEGIN
|
|
RETURN SHORT(longsin(LONG(x)));
|
|
END sin;
|
|
|
|
PROCEDURE longsin(x: LONGREAL): LONGREAL;
|
|
CONST
|
|
a2 = -0.1666666664D;
|
|
a4 = 0.0083333315D;
|
|
a6 = -0.0001984090D;
|
|
a8 = 0.0000027526D;
|
|
a10 = -0.0000000239D;
|
|
VAR
|
|
xsqr: LONGREAL;
|
|
neg: BOOLEAN;
|
|
BEGIN
|
|
neg := FALSE;
|
|
IF x < 0.0D THEN
|
|
neg := TRUE;
|
|
x := -x;
|
|
END;
|
|
x := FIF(x, 1.0D / LONG(twicepi), (* dummy *) xsqr) * LONG(twicepi);
|
|
IF x >= LONG(pi) THEN
|
|
neg := NOT neg;
|
|
x := x - LONG(pi);
|
|
END;
|
|
IF x > LONG(halfpi) THEN
|
|
x := LONG(pi) - x;
|
|
END;
|
|
xsqr := x * x;
|
|
x := x * (1.0D + xsqr*(a2+xsqr*(a4+xsqr*(a6+xsqr*(a8+xsqr*a10)))));
|
|
IF neg THEN RETURN -x; END;
|
|
RETURN x;
|
|
END longsin;
|
|
|
|
PROCEDURE cos(x: REAL): REAL;
|
|
BEGIN
|
|
RETURN SHORT(longcos(LONG(x)));
|
|
END cos;
|
|
|
|
PROCEDURE longcos(x: LONGREAL): LONGREAL;
|
|
CONST
|
|
a2 = -0.4999999963D;
|
|
a4 = 0.0416666418D;
|
|
a6 = -0.0013888397D;
|
|
a8 = 0.0000247609D;
|
|
a10 = -0.0000002605D;
|
|
VAR
|
|
xsqr: LONGREAL;
|
|
neg: BOOLEAN;
|
|
BEGIN
|
|
neg := FALSE;
|
|
IF x < 0.0D THEN x := -x; END;
|
|
x := FIF(x, 1.0D / LONG(twicepi), (* dummy *) xsqr) * LONG(twicepi);
|
|
IF x >= LONG(pi) THEN
|
|
x := LONG(twicepi) - x;
|
|
END;
|
|
IF x > LONG(halfpi) THEN
|
|
neg := NOT neg;
|
|
x := LONG(pi) - x;
|
|
END;
|
|
xsqr := x * x;
|
|
x := 1.0D + xsqr*(a2+xsqr*(a4+xsqr*(a6+xsqr*(a8+xsqr*a10))));
|
|
IF neg THEN RETURN -x; END;
|
|
RETURN x;
|
|
END longcos;
|
|
|
|
PROCEDURE tan(x: REAL): REAL;
|
|
BEGIN
|
|
RETURN SHORT(longtan(LONG(x)));
|
|
END tan;
|
|
|
|
PROCEDURE longtan(x: LONGREAL): LONGREAL;
|
|
VAR cosinus: LONGREAL;
|
|
BEGIN
|
|
cosinus := longcos(x);
|
|
IF cosinus = 0.0D THEN
|
|
Message("tan: result does not exist");
|
|
HALT
|
|
END;
|
|
RETURN longsin(x)/cosinus;
|
|
END longtan;
|
|
|
|
PROCEDURE arcsin(x: REAL): REAL;
|
|
BEGIN
|
|
RETURN SHORT(longarcsin(LONG(x)));
|
|
END arcsin;
|
|
|
|
PROCEDURE longarcsin(x: LONGREAL): LONGREAL;
|
|
CONST
|
|
a0 = 1.5707963050D;
|
|
a1 = -0.2145988016D;
|
|
a2 = 0.0889789874D;
|
|
a3 = -0.0501743046D;
|
|
a4 = 0.0308918810D;
|
|
a5 = -0.0170881256D;
|
|
a6 = 0.0066700901D;
|
|
a7 = -0.0012624911D;
|
|
BEGIN
|
|
IF x < 0.0D THEN x := -x; END;
|
|
IF x > 1.0D THEN
|
|
Message("arcsin: argument > 1");
|
|
HALT
|
|
END;
|
|
RETURN LONG(halfpi) -
|
|
longsqrt(1.0D - x)*(a0+x*(a1+x*(a2+x*(a3+x*(a4+x*(a5+x*(a6+x*a7)))))));
|
|
END longarcsin;
|
|
|
|
PROCEDURE arccos(x: REAL): REAL;
|
|
BEGIN
|
|
RETURN SHORT(longarccos(LONG(x)));
|
|
END arccos;
|
|
|
|
PROCEDURE longarccos(x: LONGREAL): LONGREAL;
|
|
BEGIN
|
|
RETURN LONG(halfpi) - longarcsin(x);
|
|
END longarccos;
|
|
|
|
PROCEDURE arctan(x: REAL): REAL;
|
|
BEGIN
|
|
RETURN SHORT(longarctan(LONG(x)));
|
|
END arctan;
|
|
|
|
PROCEDURE longarctan(x: LONGREAL): LONGREAL;
|
|
CONST
|
|
a2 = -0.3333314528D;
|
|
a4 = 0.1999355085D;
|
|
a6 = -0.1420889944D;
|
|
a8 = 0.1065626393D;
|
|
a10 = -0.0752896400D;
|
|
a12 = 0.0429096318D;
|
|
a14 = -0.0161657367D;
|
|
a16 = 0.0028662257D;
|
|
VAR
|
|
xsqr: LONGREAL;
|
|
rev: BOOLEAN;
|
|
neg: BOOLEAN;
|
|
BEGIN
|
|
rev := FALSE;
|
|
neg := FALSE;
|
|
IF x < 0.0D THEN
|
|
neg := TRUE;
|
|
x := -x;
|
|
END;
|
|
IF x > 1.0D THEN
|
|
rev := TRUE;
|
|
x := 1.0D / x;
|
|
END;
|
|
xsqr := x * x;
|
|
x := x * (1.0D + xsqr*(a2+xsqr*(a4+xsqr*(a6+xsqr*(a8+xsqr*(a10+xsqr*(a12+xsqr*(a14+xsqr*a16))))))));
|
|
IF rev THEN
|
|
x := LONG(quartpi) - x;
|
|
END;
|
|
IF neg THEN RETURN -x; END;
|
|
RETURN x;
|
|
END longarctan;
|
|
|
|
(* hyperbolic functions *)
|
|
|
|
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;
|
|
|
|
END Mathlib.
|