(*
BSD 2-Clause License
Copyright (c) 2013-2014, 2018-2022 Anton Krotov
All rights reserved.
*)
MODULE Math;
IMPORT SYSTEM;
CONST
pi* = 3.141592653589793;
e* = 2.718281828459045;
PROCEDURE IsNan* (x: REAL): BOOLEAN;
VAR
h, l: SET;
BEGIN
SYSTEM.GET(SYSTEM.ADR(x), l);
SYSTEM.GET(SYSTEM.ADR(x) + 4, h)
RETURN (h * {20..30} = {20..30}) & ((h * {0..19} # {}) OR (l * {0..31} # {}))
END IsNan;
PROCEDURE IsInf* (x: REAL): BOOLEAN;
RETURN ABS(x) = SYSTEM.INF()
END IsInf;
PROCEDURE Max (a, b: REAL): REAL;
VAR
res: REAL;
BEGIN
IF a > b THEN
res := a
ELSE
res := b
END
RETURN res
END Max;
PROCEDURE Min (a, b: REAL): REAL;
VAR
res: REAL;
BEGIN
IF a < b THEN
res := a
ELSE
res := b
END
RETURN res
END Min;
PROCEDURE SameValue (a, b: REAL): BOOLEAN;
VAR
eps: REAL;
res: BOOLEAN;
BEGIN
eps := Max(Min(ABS(a), ABS(b)) * 1.0E-12, 1.0E-12);
IF a > b THEN
res := (a - b) <= eps
ELSE
res := (b - a) <= eps
END
RETURN res
END SameValue;
PROCEDURE IsZero (x: REAL): BOOLEAN;
RETURN ABS(x) <= 1.0E-12
END IsZero;
PROCEDURE [stdcall] sqrt* (x: REAL): REAL;
BEGIN
SYSTEM.CODE(
0DDH, 045H, 008H, (* fld qword [ebp + 08h] *)
0D9H, 0FAH, (* fsqrt *)
0C9H, (* leave *)
0C2H, 008H, 000H (* ret 08h *)
)
RETURN 0.0
END sqrt;
PROCEDURE [stdcall] sin* (x: REAL): REAL;
BEGIN
SYSTEM.CODE(
0DDH, 045H, 008H, (* fld qword [ebp + 08h] *)
0D9H, 0FEH, (* fsin *)
0C9H, (* leave *)
0C2H, 008H, 000H (* ret 08h *)
)
RETURN 0.0
END sin;
PROCEDURE [stdcall] cos* (x: REAL): REAL;
BEGIN
SYSTEM.CODE(
0DDH, 045H, 008H, (* fld qword [ebp + 08h] *)
0D9H, 0FFH, (* fcos *)
0C9H, (* leave *)
0C2H, 008H, 000H (* ret 08h *)
)
RETURN 0.0
END cos;
PROCEDURE [stdcall] tan* (x: REAL): REAL;
BEGIN
SYSTEM.CODE(
0DDH, 045H, 008H, (* fld qword [ebp + 08h] *)
0D9H, 0FBH, (* fsincos *)
0DEH, 0F9H, (* fdivp st1, st *)
0C9H, (* leave *)
0C2H, 008H, 000H (* ret 08h *)
)
RETURN 0.0
END tan;
PROCEDURE [stdcall] arctan2* (y, x: REAL): REAL;
BEGIN
SYSTEM.CODE(
0DDH, 045H, 008H, (* fld qword [ebp + 08h] *)
0DDH, 045H, 010H, (* fld qword [ebp + 10h] *)
0D9H, 0F3H, (* fpatan *)
0C9H, (* leave *)
0C2H, 010H, 000H (* ret 10h *)
)
RETURN 0.0
END arctan2;
PROCEDURE [stdcall] ln* (x: REAL): REAL;
BEGIN
SYSTEM.CODE(
0D9H, 0EDH, (* fldln2 *)
0DDH, 045H, 008H, (* fld qword [ebp + 08h] *)
0D9H, 0F1H, (* fyl2x *)
0C9H, (* leave *)
0C2H, 008H, 000H (* ret 08h *)
)
RETURN 0.0
END ln;
PROCEDURE [stdcall] log* (base, x: REAL): REAL;
BEGIN
SYSTEM.CODE(
0D9H, 0E8H, (* fld1 *)
0DDH, 045H, 010H, (* fld qword [ebp + 10h] *)
0D9H, 0F1H, (* fyl2x *)
0D9H, 0E8H, (* fld1 *)
0DDH, 045H, 008H, (* fld qword [ebp + 08h] *)
0D9H, 0F1H, (* fyl2x *)
0DEH, 0F9H, (* fdivp st1, st *)
0C9H, (* leave *)
0C2H, 010H, 000H (* ret 10h *)
)
RETURN 0.0
END log;
PROCEDURE [stdcall] exp* (x: REAL): REAL;
BEGIN
SYSTEM.CODE(
0DDH, 045H, 008H, (* fld qword [ebp + 08h] *)
0D9H, 0EAH, (* fldl2e *)
0DEH, 0C9H, 0D9H, 0C0H,
0D9H, 0FCH, 0DCH, 0E9H,
0D9H, 0C9H, 0D9H, 0F0H,
0D9H, 0E8H, 0DEH, 0C1H,
0D9H, 0FDH, 0DDH, 0D9H,
0C9H, (* leave *)
0C2H, 008H, 000H (* ret 08h *)
)
RETURN 0.0
END exp;
PROCEDURE [stdcall] round* (x: REAL): REAL;
BEGIN
SYSTEM.CODE(
0DDH, 045H, 008H, (* fld qword [ebp + 08h] *)
0D9H, 07DH, 0F4H, 0D9H,
07DH, 0F6H, 066H, 081H,
04DH, 0F6H, 000H, 003H,
0D9H, 06DH, 0F6H, 0D9H,
0FCH, 0D9H, 06DH, 0F4H,
0C9H, (* leave *)
0C2H, 008H, 000H (* ret 08h *)
)
RETURN 0.0
END round;
PROCEDURE [stdcall] frac* (x: REAL): REAL;
BEGIN
SYSTEM.CODE(
050H,
0DDH, 045H, 008H, (* fld qword [ebp + 08h] *)
0D9H, 0C0H, 0D9H, 03CH,
024H, 0D9H, 07CH, 024H,
002H, 066H, 081H, 04CH,
024H, 002H, 000H, 00FH,
0D9H, 06CH, 024H, 002H,
0D9H, 0FCH, 0D9H, 02CH,
024H, 0DEH, 0E9H,
0C9H, (* leave *)
0C2H, 008H, 000H (* ret 08h *)
)
RETURN 0.0
END frac;
PROCEDURE sqri* (x: INTEGER): INTEGER;
RETURN x * x
END sqri;
PROCEDURE sqrr* (x: REAL): REAL;
RETURN x * x
END sqrr;
PROCEDURE arcsin* (x: REAL): REAL;
RETURN arctan2(x, sqrt(1.0 - x * x))
END arcsin;
PROCEDURE arccos* (x: REAL): REAL;
RETURN arctan2(sqrt(1.0 - x * x), x)
END arccos;
PROCEDURE arctan* (x: REAL): REAL;
RETURN arctan2(x, 1.0)
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 := 1.0 - 2.0 / (exp(2.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;
RETURN ln(x + sqrt(x * x - 1.0))
END arcosh;
PROCEDURE artanh* (x: REAL): REAL;
VAR
res: REAL;
BEGIN
IF SameValue(x, 1.0) THEN
res := SYSTEM.INF()
ELSIF SameValue(x, -1.0) THEN
res := -SYSTEM.INF()
ELSE
res := 0.5 * ln((1.0 + x) / (1.0 - x))
END
RETURN res
END artanh;
PROCEDURE floor* (x: REAL): REAL;
VAR
f: REAL;
BEGIN
f := frac(x);
x := x - f;
IF f < 0.0 THEN
x := x - 1.0
END
RETURN x
END floor;
PROCEDURE ceil* (x: REAL): REAL;
VAR
f: REAL;
BEGIN
f := frac(x);
x := x - f;
IF f > 0.0 THEN
x := x + 1.0
END
RETURN x
END ceil;
PROCEDURE power* (base, exponent: REAL): REAL;
VAR
res: REAL;
BEGIN
IF exponent = 0.0 THEN
res := 1.0
ELSIF (base = 0.0) & (exponent > 0.0) THEN
res := 0.0
ELSE
res := exp(exponent * ln(base))
END
RETURN res
END power;
PROCEDURE ipower* (base: REAL; exponent: INTEGER): REAL;
VAR
i: INTEGER;
a: REAL;
BEGIN
a := 1.0;
IF base # 0.0 THEN
IF exponent # 0 THEN
IF exponent < 0 THEN
base := 1.0 / base
END;
i := ABS(exponent);
WHILE i > 0 DO
WHILE ~ODD(i) DO
i := LSR(i, 1);
base := sqrr(base)
END;
DEC(i);
a := a * base
END
ELSE
a := 1.0
END
ELSE
ASSERT(exponent > 0);
a := 0.0
END
RETURN a
END ipower;
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 DegToRad* (x: REAL): REAL;
RETURN x * (pi / 180.0)
END DegToRad;
PROCEDURE RadToDeg* (x: REAL): REAL;
RETURN x * (180.0 / pi)
END RadToDeg;
(* Return hypotenuse of triangle *)
PROCEDURE hypot* (x, y: REAL): REAL;
VAR
a: REAL;
BEGIN
x := ABS(x);
y := ABS(y);
IF x > y THEN
a := x * sqrt(1.0 + sqrr(y / x))
ELSE
IF x > 0.0 THEN
a := y * sqrt(1.0 + sqrr(x / y))
ELSE
a := y
END
END
RETURN a
END hypot;
END Math.