573 lines
13 KiB
Modula-2
573 lines
13 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;
|
|
|
|
CONST
|
|
OneRadianInDegrees = 57.295779513082320876798155D;
|
|
OneDegreeInRadians = 0.017453292519943295769237D;
|
|
Sqrt2 = 1.41421356237309504880168872420969808D;
|
|
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 exp(x: REAL): REAL;
|
|
BEGIN
|
|
RETURN SHORT(longexp(LONG(x)));
|
|
END exp;
|
|
|
|
PROCEDURE longexp(x: LONGREAL): LONGREAL;
|
|
(* 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 *)
|
|
CONST
|
|
p0 = 0.2080384346694663001443843411D+07;
|
|
p1 = 0.3028697169744036299076048876D+05;
|
|
p2 = 0.6061485330061080841615584556D+02;
|
|
q0 = 0.6002720360238832528230907598D+07;
|
|
q1 = 0.3277251518082914423057964422D+06;
|
|
q2 = 0.1749287689093076403844945335D+04;
|
|
q3 = 0.1000000000000000000000000000D+01;
|
|
|
|
VAR
|
|
neg: BOOLEAN;
|
|
xPxx, Qxx: LONGREAL;
|
|
n: LONGREAL;
|
|
n1 : INTEGER;
|
|
xsq : LONGREAL;
|
|
large: BOOLEAN;
|
|
BEGIN
|
|
neg := x < 0.0D;
|
|
IF neg THEN
|
|
x := -x;
|
|
END;
|
|
x := FIF(x/longln2, 1.0D, n);
|
|
large := x > 0.5D;
|
|
IF large THEN x := x - 0.5D; END;
|
|
xsq := x*x;
|
|
xPxx := x*((p2*xsq+p1)*xsq+p0);
|
|
Qxx := ((q3*xsq+q2)*xsq+q1)*xsq+q0;
|
|
x := (Qxx + xPxx)/(Qxx - xPxx);
|
|
IF large THEN
|
|
x := x * Sqrt2;
|
|
END;
|
|
n1 := TRUNCD(n + 0.5D);
|
|
WHILE n1 >= 16 DO
|
|
x := x * 65536.0D;
|
|
n1 := n1 - 16;
|
|
END;
|
|
WHILE n1 > 0 DO
|
|
x := x * 2.0D;
|
|
DEC(n1);
|
|
END;
|
|
IF neg THEN RETURN 1.0D/x; END;
|
|
RETURN x;
|
|
END longexp;
|
|
|
|
PROCEDURE ln(x: REAL): REAL; (* natural log *)
|
|
BEGIN
|
|
RETURN SHORT(longln(LONG(x)));
|
|
END ln;
|
|
|
|
PROCEDURE longln(x: LONGREAL): LONGREAL; (* natural log *)
|
|
(* log(x) = z*P(z*z)/Q(z*z), z = (x-1)/(x+1), x in [1/sqrt(2), sqrt(2)]
|
|
Hart & Cheney #2707
|
|
*)
|
|
CONST
|
|
p0 = 0.7504094990777122217455611007D+02;
|
|
p1 = -0.1345669115050430235318253537D+03;
|
|
p2 = 0.7413719213248602512779336470D+02;
|
|
p3 = -0.1277249755012330819984385000D+02;
|
|
p4 = 0.3327108381087686938144000000D+00;
|
|
q0 = 0.3752047495388561108727775374D+02;
|
|
q1 = -0.7979028073715004879439951583D+02;
|
|
q2 = 0.5616126132118257292058560360D+02;
|
|
q3 = -0.1450868091858082685362325000D+02;
|
|
q4 = 0.1000000000000000000000000000D+01;
|
|
VAR
|
|
exp: INTEGER;
|
|
z, zsq: LONGREAL;
|
|
|
|
BEGIN
|
|
IF x <= 0.0D THEN
|
|
Message("ln: argument <= 0");
|
|
HALT
|
|
END;
|
|
x := FEF(x, exp);
|
|
WHILE x < OneOverSqrt2 DO
|
|
x := x + x;
|
|
DEC(exp);
|
|
END;
|
|
z := (x - 1.0D) / (x + 1.0D);
|
|
zsq := z*z;
|
|
RETURN z * ((((p4*zsq+p3)*zsq+p2)*zsq+p1)*zsq+p0) /
|
|
((((q4*zsq+q3)*zsq+q2)*zsq+q1)*zsq+q0) +
|
|
FLOATD(exp) * longln2;
|
|
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; quadrant: INTEGER) : LONGREAL;
|
|
(* sin(0.5*pi*x) = x * P(x*x)/Q(x*x) for x in [0,1]
|
|
Hart & Cheney # 3374
|
|
*)
|
|
CONST
|
|
p0 = 0.4857791909822798473837058825D+10;
|
|
p1 = -0.1808816670894030772075877725D+10;
|
|
p2 = 0.1724314784722489597789244188D+09;
|
|
p3 = -0.6351331748520454245913645971D+07;
|
|
p4 = 0.1002087631419532326179108883D+06;
|
|
p5 = -0.5830988897678192576148973679D+03;
|
|
q0 = 0.3092566379840468199410228418D+10;
|
|
q1 = 0.1202384907680254190870913060D+09;
|
|
q2 = 0.2321427631602460953669856368D+07;
|
|
q3 = 0.2848331644063908832127222835D+05;
|
|
q4 = 0.2287602116741682420054505174D+03;
|
|
q5 = 0.1000000000000000000000000000D+01;
|
|
A1 = 6.2822265625D;
|
|
A2 = 0.00095874467958647692528676655900576D;
|
|
VAR
|
|
xsq, x1, x2, n : LONGREAL;
|
|
t : INTEGER;
|
|
BEGIN
|
|
IF x < 0.0D THEN
|
|
INC(quadrant, 2);
|
|
x := -x;
|
|
END;
|
|
IF longhalfpi - x = longhalfpi THEN
|
|
CASE quadrant OF
|
|
| 0,2:
|
|
RETURN 0.0D;
|
|
| 1:
|
|
RETURN 1.0D;
|
|
| 3:
|
|
RETURN -1.0D;
|
|
END;
|
|
END;
|
|
IF x >= longtwicepi THEN
|
|
IF x <= FLOATD(MAX(LONGINT)) THEN
|
|
(* Use extended precision to calculate reduced argument.
|
|
Split 2pi in 2 parts a1 and a2, of which the first only
|
|
uses some bits of the mantissa, so that n * a1 is
|
|
exactly representable, where n is the integer part of
|
|
x/pi.
|
|
Here we used 12 bits of the mantissa for a1.
|
|
Also split x in integer part x1 and fraction part x2.
|
|
We then compute x-n*2pi as ((x1 - n*a1) + x2) - n*a2.
|
|
*)
|
|
n := FLOATD(TRUNCD(x/longtwicepi));
|
|
x1 := FLOATD(TRUNCD(x));
|
|
x2 := x - x1;
|
|
x := ((x1 - n * A1) + x2) - n * A2;
|
|
ELSE
|
|
x := FIF(x/longtwicepi, 1.0D, x1) * longtwicepi;
|
|
END
|
|
END;
|
|
x := x / longhalfpi;
|
|
t := TRUNC(x);
|
|
x := x - FLOATD(t);
|
|
quadrant := (quadrant + t MOD 4) MOD 4;
|
|
IF ODD(quadrant) THEN
|
|
x := 1.0D - x;
|
|
END;
|
|
IF quadrant > 1 THEN
|
|
x := -x;
|
|
END;
|
|
xsq := x * x;
|
|
RETURN x * (((((p5*xsq+p4)*xsq+p3)*xsq+p2)*xsq+p1)*xsq+p0) /
|
|
(((((q5*xsq+q4)*xsq+q3)*xsq+q2)*xsq+q1)*xsq+q0);
|
|
END sinus;
|
|
|
|
PROCEDURE longsin(x: LONGREAL): LONGREAL;
|
|
BEGIN
|
|
RETURN sinus(x, 0);
|
|
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, 1);
|
|
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 arcsincos(x: LONGREAL; cosfl: BOOLEAN): LONGREAL;
|
|
VAR
|
|
negative : BOOLEAN;
|
|
BEGIN
|
|
negative := x <= 0.0D;
|
|
IF negative THEN x := -x; END;
|
|
IF x > 1.0D THEN
|
|
Message("arcsin or arccos: argument > 1");
|
|
HALT
|
|
END;
|
|
IF x = 1.0D THEN
|
|
x := longhalfpi;
|
|
ELSE
|
|
x := longarctan(x/longsqrt(1.0D - x*x));
|
|
END;
|
|
IF negative THEN x := -x; END;
|
|
IF cosfl THEN
|
|
RETURN longhalfpi - 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;
|
|
|
|
TYPE
|
|
precomputed = RECORD
|
|
X: LONGREAL; (* partition point *)
|
|
arctan: LONGREAL; (* arctan of evaluation node *)
|
|
OneOverXn: LONGREAL; (* 1/xn *)
|
|
OneOverXnSquarePlusone: LONGREAL; (* ... *)
|
|
END;
|
|
|
|
VAR arctaninit: BOOLEAN;
|
|
precomp : ARRAY[0..4] OF precomputed;
|
|
|
|
PROCEDURE longarctan(x: LONGREAL): LONGREAL;
|
|
(* The interval [0, infinity) is treated as follows:
|
|
Define partition points Xi
|
|
X0 = 0
|
|
X1 = tan(pi/16)
|
|
X2 = tan(3pi/16)
|
|
X3 = tan(5pi/16)
|
|
X4 = tan(7pi/16)
|
|
X5 = infinity
|
|
and evaluation nodes xi
|
|
x2 = tan(2pi/16)
|
|
x3 = tan(4pi/16)
|
|
x4 = tan(6pi/16)
|
|
x5 = infinity
|
|
An argument x in [Xn-1, Xn] is now reduced to an argument
|
|
t in [-X1, X1] by the following formulas:
|
|
|
|
t = 1/xn - (1/(xn*xn) + 1)/((1/xn) + x)
|
|
|
|
arctan(x) = arctan(xi) + arctan(t)
|
|
|
|
For the interval [0, p/16] an approximation is used:
|
|
arctan(x) = x * P(x*x)/Q(x*x)
|
|
*)
|
|
(* Hart & Cheney # 5037 *)
|
|
CONST
|
|
p0 = 0.7698297257888171026986294745D+03;
|
|
p1 = 0.1557282793158363491416585283D+04;
|
|
p2 = 0.1033384651675161628243434662D+04;
|
|
p3 = 0.2485841954911840502660889866D+03;
|
|
p4 = 0.1566564964979791769948970100D+02;
|
|
q0 = 0.7698297257888171026986294911D+03;
|
|
q1 = 0.1813892701754635858982709369D+04;
|
|
q2 = 0.1484049607102276827437401170D+04;
|
|
q3 = 0.4904645326203706217748848797D+03;
|
|
q4 = 0.5593479839280348664778328000D+02;
|
|
q5 = 0.1000000000000000000000000000D+01;
|
|
VAR
|
|
xsqr: LONGREAL;
|
|
neg: BOOLEAN;
|
|
i: INTEGER;
|
|
BEGIN
|
|
IF NOT arctaninit THEN
|
|
arctaninit := TRUE;
|
|
WITH precomp[0] DO
|
|
X := 0.19891236737965800691159762264467622;
|
|
arctan := 0.0D;
|
|
OneOverXn := 0.0D;
|
|
OneOverXnSquarePlusone := 0.0D;
|
|
END;
|
|
WITH precomp[1] DO
|
|
X := 0.66817863791929891999775768652308076;
|
|
arctan := longpi/8.0D;
|
|
OneOverXn := 2.41421356237309504880168872420969808;
|
|
OneOverXnSquarePlusone := 6.82842712474619009760337744841939616;
|
|
END;
|
|
WITH precomp[2] DO
|
|
X := 1.49660576266548901760113513494247691;
|
|
arctan := longquartpi;
|
|
OneOverXn := 1.0;
|
|
OneOverXnSquarePlusone := 2.0;
|
|
END;
|
|
WITH precomp[3] DO
|
|
X := 5.02733949212584810451497507106407238;
|
|
arctan := 3.0D*longpi/8.0D;
|
|
OneOverXn := 0.41421356237309504880168872420969808;
|
|
OneOverXnSquarePlusone := 1.17157287525380998659662255158060384;
|
|
END;
|
|
WITH precomp[4] DO
|
|
X := 0.0;
|
|
arctan := longhalfpi;
|
|
OneOverXn := 0.0;
|
|
OneOverXnSquarePlusone := 1.0;
|
|
END;
|
|
END;
|
|
neg := FALSE;
|
|
IF x < 0.0D THEN
|
|
neg := TRUE;
|
|
x := -x;
|
|
END;
|
|
i := 0;
|
|
WHILE (i <= 3) AND (x <= precomp[i].X) DO
|
|
INC(i);
|
|
END;
|
|
IF (i # 0) THEN
|
|
WITH precomp[i] DO
|
|
x := arctan + longarctan(OneOverXn-OneOverXnSquarePlusone/(OneOverXn+x));
|
|
END
|
|
ELSE
|
|
xsqr := x * x;
|
|
x := x * ((((p4*xsqr+p3)*xsqr+p2)*xsqr+p1)*xsqr+p0) /
|
|
(((((q5*xsqr+q4)*xsqr+q3)*xsqr+q2)*xsqr+q1)*xsqr+q0);
|
|
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;
|
|
|
|
BEGIN
|
|
arctaninit := FALSE;
|
|
END Mathlib.
|