(*
BSD 2-Clause License
Copyright (c) 2019, Anton Krotov
All rights reserved.
*)
MODULE Math;
IMPORT SYSTEM;
CONST
e *= 2.71828182845904523;
pi *= 3.14159265358979324;
ln2 *= 0.693147180559945309;
eps = 1.0E-16;
MaxCosArg = 1000000.0 * pi;
VAR
Exp: ARRAY 710 OF REAL;
PROCEDURE [stdcall64] sqrt* (x: REAL): REAL;
BEGIN
ASSERT(x >= 0.0);
SYSTEM.CODE(
0F2H, 0FH, 51H, 45H, 10H, (* sqrtsd xmm0, qword[rbp + 10h] *)
05DH, (* pop rbp *)
0C2H, 08H, 00H (* ret 8 *)
)
RETURN 0.0
END sqrt;
PROCEDURE exp* (x: REAL): REAL;
CONST
e25 = 1.284025416687741484; (* exp(0.25) *)
VAR
a, s, res: REAL;
neg: BOOLEAN;
n: INTEGER;
BEGIN
neg := x < 0.0;
IF neg THEN
x := -x
END;
IF x < FLT(LEN(Exp)) THEN
res := Exp[FLOOR(x)];
x := x - FLT(FLOOR(x));
WHILE x >= 0.25 DO
res := res * e25;
x := x - 0.25
END
ELSE
res := SYSTEM.INF();
x := 0.0
END;
n := 0;
a := 1.0;
s := 1.0;
REPEAT
INC(n);
a := a * x / FLT(n);
s := s + a
UNTIL a < eps;
IF neg THEN
res := 1.0 / (res * s)
ELSE
res := res * s
END
RETURN res
END exp;
PROCEDURE ln* (x: REAL): REAL;
VAR
a, x2, res: REAL;
n: INTEGER;
BEGIN
ASSERT(x > 0.0);
UNPK(x, n);
x := (x - 1.0) / (x + 1.0);
x2 := x * x;
res := x + FLT(n) * (ln2 * 0.5);
n := 1;
REPEAT
INC(n, 2);
x := x * x2;
a := x / FLT(n);
res := res + a
UNTIL a < eps
RETURN res * 2.0
END ln;
PROCEDURE power* (base, exponent: REAL): REAL;
BEGIN
ASSERT(base > 0.0)
RETURN exp(exponent * ln(base))
END power;
PROCEDURE log* (base, x: REAL): REAL;
BEGIN
ASSERT(base > 0.0);
ASSERT(x > 0.0)
RETURN ln(x) / ln(base)
END log;
PROCEDURE cos* (x: REAL): REAL;
VAR
a, res: REAL;
n: INTEGER;
BEGIN
x := ABS(x);
ASSERT(x <= MaxCosArg);
x := x - FLT( FLOOR(x / (2.0 * pi)) ) * (2.0 * pi);
x := x * x;
res := 0.0;
a := 1.0;
n := -1;
REPEAT
INC(n, 2);
res := res + a;
a := -a * x / FLT(n*n + n)
UNTIL ABS(a) < eps
RETURN res
END cos;
PROCEDURE sin* (x: REAL): REAL;
BEGIN
ASSERT(ABS(x) <= MaxCosArg);
x := cos(x)
RETURN sqrt(1.0 - x * x)
END sin;
PROCEDURE tan* (x: REAL): REAL;
BEGIN
ASSERT(ABS(x) <= MaxCosArg);
x := cos(x)
RETURN sqrt(1.0 - x * x) / x
END tan;
PROCEDURE arcsin* (x: REAL): REAL;
PROCEDURE arctan (x: REAL): REAL;
VAR
z, p, k: REAL;
BEGIN
p := x / (x * x + 1.0);
z := p * x;
x := 0.0;
k := 0.0;
REPEAT
k := k + 2.0;
x := x + p;
p := p * k * z / (k + 1.0)
UNTIL p < eps
RETURN x
END arctan;
BEGIN
ASSERT(ABS(x) <= 1.0);
IF ABS(x) >= 0.707 THEN
x := 0.5 * pi - arctan(sqrt(1.0 - x * x) / x)
ELSE
x := arctan(x / sqrt(1.0 - x * x))
END
RETURN x
END arcsin;
PROCEDURE arccos* (x: REAL): REAL;
BEGIN
ASSERT(ABS(x) <= 1.0)
RETURN 0.5 * pi - arcsin(x)
END arccos;
PROCEDURE arctan* (x: REAL): REAL;
RETURN arcsin(x / sqrt(1.0 + x * x))
END arctan;
PROCEDURE sinh* (x: REAL): REAL;
BEGIN
x := exp(x)
RETURN (x - 1.0 / x) * 0.5
END sinh;
PROCEDURE cosh* (x: REAL): REAL;
BEGIN
x := exp(x)
RETURN (x + 1.0 / x) * 0.5
END cosh;
PROCEDURE tanh* (x: REAL): REAL;
BEGIN
IF x > 15.0 THEN
x := 1.0
ELSIF x < -15.0 THEN
x := -1.0
ELSE
x := exp(2.0 * x);
x := (x - 1.0) / (x + 1.0)
END
RETURN x
END tanh;
PROCEDURE arsinh* (x: REAL): REAL;
RETURN ln(x + sqrt(x * x + 1.0))
END arsinh;
PROCEDURE arcosh* (x: REAL): REAL;
BEGIN
ASSERT(x >= 1.0)
RETURN ln(x + sqrt(x * x - 1.0))
END arcosh;
PROCEDURE artanh* (x: REAL): REAL;
BEGIN
ASSERT(ABS(x) < 1.0)
RETURN 0.5 * ln((1.0 + x) / (1.0 - x))
END artanh;
PROCEDURE sgn* (x: REAL): INTEGER;
VAR
res: INTEGER;
BEGIN
IF x > 0.0 THEN
res := 1
ELSIF x < 0.0 THEN
res := -1
ELSE
res := 0
END
RETURN res
END sgn;
PROCEDURE fact* (n: INTEGER): REAL;
VAR
res: REAL;
BEGIN
res := 1.0;
WHILE n > 1 DO
res := res * FLT(n);
DEC(n)
END
RETURN res
END fact;
PROCEDURE init;
VAR
i: INTEGER;
BEGIN
Exp[0] := 1.0;
FOR i := 1 TO LEN(Exp) - 1 DO
Exp[i] := Exp[i - 1] * e
END
END init;
BEGIN
init
END Math.