450 lines
9.0 KiB
Plaintext
450 lines
9.0 KiB
Plaintext
(*
|
|
BSD 2-Clause License
|
|
|
|
Copyright (c) 2013-2014, 2018-2020 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 := 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;
|
|
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. |