338 lines
6.3 KiB
Modula-2
338 lines
6.3 KiB
Modula-2
IMPLEMENTATION MODULE MathLib0;
|
|
(* Rewritten in Modula-2.
|
|
The originals came from the Pascal runtime library.
|
|
*)
|
|
|
|
FROM FIFFEF IMPORT FIF, FEF;
|
|
|
|
CONST
|
|
HUGE = 1.701411733192644270E38;
|
|
|
|
PROCEDURE sinus(arg: REAL; quad: INTEGER): REAL;
|
|
|
|
(* Coefficients for sin/cos are #3370 from Hart & Cheney (18.80D).
|
|
*)
|
|
CONST
|
|
twoopi = 0.63661977236758134308;
|
|
p0 = 0.1357884097877375669092680E8;
|
|
p1 = -0.4942908100902844161158627E7;
|
|
p2 = 0.4401030535375266501944918E6;
|
|
p3 = -0.1384727249982452873054457E5;
|
|
p4 = 0.1459688406665768722226959E3;
|
|
q0 = 0.8644558652922534429915149E7;
|
|
q1 = 0.4081792252343299749395779E6;
|
|
q2 = 0.9463096101538208180571257E4;
|
|
q3 = 0.1326534908786136358911494E3;
|
|
VAR
|
|
e, f: REAL;
|
|
ysq: REAL;
|
|
x,y: REAL;
|
|
k: INTEGER;
|
|
temp1, temp2: REAL;
|
|
BEGIN
|
|
x := arg;
|
|
IF x < 0.0 THEN
|
|
x := -x;
|
|
quad := quad + 2;
|
|
END;
|
|
x := x*twoopi; (*underflow?*)
|
|
IF x>32764.0 THEN
|
|
y := FIF(x, 10.0, e);
|
|
e := e + FLOAT(quad);
|
|
temp1 := FIF(0.25, e, f);
|
|
quad := TRUNC(e - 4.0*f);
|
|
ELSE
|
|
k := TRUNC(x);
|
|
y := x - FLOAT(k);
|
|
quad := (quad + k) MOD 4;
|
|
END;
|
|
IF ODD(quad) THEN
|
|
y := 1.0-y;
|
|
END;
|
|
IF quad > 1 THEN
|
|
y := -y;
|
|
END;
|
|
|
|
ysq := y*y;
|
|
temp1 := ((((p4*ysq+p3)*ysq+p2)*ysq+p1)*ysq+p0)*y;
|
|
temp2 := ((((ysq+q3)*ysq+q2)*ysq+q1)*ysq+q0);
|
|
RETURN temp1/temp2;
|
|
END sinus;
|
|
|
|
PROCEDURE cos(arg: REAL): REAL;
|
|
BEGIN
|
|
IF arg < 0.0 THEN
|
|
arg := -arg;
|
|
END;
|
|
RETURN sinus(arg, 1);
|
|
END cos;
|
|
|
|
PROCEDURE sin(arg: REAL): REAL;
|
|
BEGIN
|
|
RETURN sinus(arg, 0);
|
|
END sin;
|
|
|
|
(*
|
|
floating-point arctangent
|
|
|
|
arctan returns the value of the arctangent of its
|
|
argument in the range [-pi/2,pi/2].
|
|
|
|
coefficients are #5077 from Hart & Cheney. (19.56D)
|
|
*)
|
|
|
|
CONST
|
|
sq2p1 = 2.414213562373095048802E0;
|
|
sq2m1 = 0.414213562373095048802E0;
|
|
pio2 = 1.570796326794896619231E0;
|
|
pio4 = 0.785398163397448309615E0;
|
|
p4 = 0.161536412982230228262E2;
|
|
p3 = 0.26842548195503973794141E3;
|
|
p2 = 0.11530293515404850115428136E4;
|
|
p1 = 0.178040631643319697105464587E4;
|
|
p0 = 0.89678597403663861959987488E3;
|
|
q4 = 0.5895697050844462222791E2;
|
|
q3 = 0.536265374031215315104235E3;
|
|
q2 = 0.16667838148816337184521798E4;
|
|
q1 = 0.207933497444540981287275926E4;
|
|
q0 = 0.89678597403663861962481162E3;
|
|
|
|
(*
|
|
xatan evaluates a series valid in the
|
|
range [-0.414...,+0.414...].
|
|
*)
|
|
|
|
PROCEDURE xatan(arg: REAL) : REAL;
|
|
VAR
|
|
argsq, value: REAL;
|
|
BEGIN
|
|
argsq := arg*arg;
|
|
value := ((((p4*argsq + p3)*argsq + p2)*argsq + p1)*argsq + p0);
|
|
value := value/(((((argsq + q4)*argsq + q3)*argsq + q2)*argsq + q1)*argsq + q0);
|
|
RETURN value*arg;
|
|
END xatan;
|
|
|
|
PROCEDURE satan(arg: REAL): REAL;
|
|
BEGIN
|
|
IF arg < sq2m1 THEN
|
|
RETURN xatan(arg);
|
|
ELSIF arg > sq2p1 THEN
|
|
RETURN pio2 - xatan(1.0/arg);
|
|
ELSE
|
|
RETURN pio4 + xatan((arg-1.0)/(arg+1.0));
|
|
END;
|
|
END satan;
|
|
|
|
(*
|
|
atan makes its argument positive and
|
|
calls the inner routine satan.
|
|
*)
|
|
|
|
PROCEDURE arctan(arg: REAL): REAL;
|
|
BEGIN
|
|
IF arg>0.0 THEN
|
|
RETURN satan(arg);
|
|
ELSE
|
|
RETURN -satan(-arg);
|
|
END;
|
|
END arctan;
|
|
|
|
(*
|
|
sqrt returns the square root of its floating
|
|
point argument. Newton's method.
|
|
*)
|
|
|
|
PROCEDURE sqrt(arg: REAL): REAL;
|
|
VAR
|
|
x, temp: REAL;
|
|
exp, i: INTEGER;
|
|
BEGIN
|
|
IF arg <= 0.0 THEN
|
|
IF arg < 0.0 THEN
|
|
(* ??? *)
|
|
;
|
|
END;
|
|
RETURN 0.0;
|
|
END;
|
|
x := FEF(arg,exp);
|
|
(*
|
|
* NOTE
|
|
* this wont work on 1's comp
|
|
*)
|
|
IF ODD(exp) THEN
|
|
x := 2.0 * x;
|
|
DEC(exp);
|
|
END;
|
|
temp := 0.5*(1.0 + x);
|
|
|
|
WHILE exp > 28 DO
|
|
temp := temp * 16384.0;
|
|
exp := exp - 28;
|
|
END;
|
|
WHILE exp < -28 DO
|
|
temp := temp / 16384.0;
|
|
exp := exp + 28;
|
|
END;
|
|
WHILE exp >= 2 DO
|
|
temp := temp * 2.0;
|
|
exp := exp - 2;
|
|
END;
|
|
WHILE exp <= -2 DO
|
|
temp := temp / 2.0;
|
|
exp := exp + 2;
|
|
END;
|
|
FOR i := 0 TO 4 DO
|
|
temp := 0.5*(temp + arg/temp);
|
|
END;
|
|
RETURN temp;
|
|
END sqrt;
|
|
|
|
(*
|
|
ln returns the natural logarithm of its floating
|
|
point argument.
|
|
|
|
The coefficients are #2705 from Hart & Cheney. (19.38D)
|
|
*)
|
|
PROCEDURE ln(arg: REAL): REAL;
|
|
CONST
|
|
log2 = 0.693147180559945309E0;
|
|
sqrto2 = 0.707106781186547524E0;
|
|
p0 = -0.240139179559210510E2;
|
|
p1 = 0.309572928215376501E2;
|
|
p2 = -0.963769093368686593E1;
|
|
p3 = 0.421087371217979714E0;
|
|
q0 = -0.120069589779605255E2;
|
|
q1 = 0.194809660700889731E2;
|
|
q2 = -0.891110902798312337E1;
|
|
VAR
|
|
x,z, zsq, temp: REAL;
|
|
exp: INTEGER;
|
|
BEGIN
|
|
IF arg <= 0.0 THEN
|
|
(* ??? *)
|
|
RETURN -HUGE;
|
|
END;
|
|
x := FEF(arg,exp);
|
|
IF x<sqrto2 THEN
|
|
x := x + x;
|
|
DEC(exp);
|
|
END;
|
|
|
|
z := (x-1.0)/(x+1.0);
|
|
zsq := z*z;
|
|
|
|
temp := ((p3*zsq + p2)*zsq + p1)*zsq + p0;
|
|
temp := temp/(((zsq + q2)*zsq + q1)*zsq + q0);
|
|
temp := temp*z + FLOAT(exp)*log2;
|
|
RETURN temp;
|
|
END ln;
|
|
|
|
(*
|
|
exp returns the exponential function of its
|
|
floating-point argument.
|
|
|
|
The coefficients are #1069 from Hart and Cheney. (22.35D)
|
|
*)
|
|
|
|
PROCEDURE floor(d: REAL): REAL;
|
|
BEGIN
|
|
IF d < 0.0 THEN
|
|
d := -d;
|
|
IF FIF(d, 1.0, d) # 0.0 THEN
|
|
d := d + 1.0;
|
|
END;
|
|
d := -d;
|
|
ELSE
|
|
IF FIF(d, 1.0, d) # 0.0 THEN
|
|
(* just ignore result of FIF *)
|
|
;
|
|
END;
|
|
END;
|
|
RETURN d;
|
|
END floor;
|
|
|
|
PROCEDURE ldexp(fr: REAL; exp: INTEGER): REAL;
|
|
VAR
|
|
neg,i: INTEGER;
|
|
BEGIN
|
|
neg := 1;
|
|
IF fr < 0.0 THEN
|
|
fr := -fr;
|
|
neg := -1;
|
|
END;
|
|
fr := FEF(fr, i);
|
|
exp := exp + i;
|
|
IF exp > 127 THEN
|
|
(* Too large. ??? *)
|
|
RETURN FLOAT(neg) * HUGE;
|
|
END;
|
|
IF exp < -127 THEN
|
|
RETURN 0.0;
|
|
END;
|
|
WHILE exp > 14 DO
|
|
fr := fr * 16384.0;
|
|
exp := exp - 14;
|
|
END;
|
|
WHILE exp < -14 DO
|
|
fr := fr / 16384.0;
|
|
exp := exp + 14;
|
|
END;
|
|
WHILE exp > 0 DO
|
|
fr := fr + fr;
|
|
DEC(exp);
|
|
END;
|
|
WHILE exp < 0 DO
|
|
fr := fr / 2.0;
|
|
INC(exp);
|
|
END;
|
|
RETURN FLOAT(neg) * fr;
|
|
END ldexp;
|
|
|
|
PROCEDURE exp(arg: REAL): REAL;
|
|
CONST
|
|
p0 = 0.2080384346694663001443843411E7;
|
|
p1 = 0.3028697169744036299076048876E5;
|
|
p2 = 0.6061485330061080841615584556E2;
|
|
q0 = 0.6002720360238832528230907598E7;
|
|
q1 = 0.3277251518082914423057964422E6;
|
|
q2 = 0.1749287689093076403844945335E4;
|
|
log2e = 1.4426950408889634073599247;
|
|
sqrt2 = 1.4142135623730950488016887;
|
|
maxf = 10000.0;
|
|
VAR
|
|
fract: REAL;
|
|
temp1, temp2, xsq: REAL;
|
|
ent: INTEGER;
|
|
BEGIN
|
|
IF arg = 0.0 THEN
|
|
RETURN 1.0;
|
|
END;
|
|
IF arg < -maxf THEN
|
|
RETURN 0.0;
|
|
END;
|
|
IF arg > maxf THEN
|
|
(* result too large ??? *)
|
|
RETURN HUGE;
|
|
END;
|
|
arg := arg * log2e;
|
|
ent := TRUNC(floor(arg));
|
|
fract := (arg-FLOAT(ent)) - 0.5;
|
|
xsq := fract*fract;
|
|
temp1 := ((p2*xsq+p1)*xsq+p0)*fract;
|
|
temp2 := ((xsq+q2)*xsq+q1)*xsq + q0;
|
|
RETURN ldexp(sqrt2*(temp2+temp1)/(temp2-temp1), ent);
|
|
END exp;
|
|
|
|
PROCEDURE entier(x: REAL): INTEGER;
|
|
BEGIN
|
|
RETURN TRUNC(x); (* ??? *)
|
|
END entier;
|
|
|
|
PROCEDURE real(x: INTEGER): REAL;
|
|
BEGIN
|
|
RETURN FLOAT(x); (* ??? *)
|
|
END real;
|
|
|
|
BEGIN
|
|
END MathLib0.
|