(*
Copyright 2016 Anton Krotov
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*)
MODULE Math;
IMPORT sys := SYSTEM;
CONST pi* = 3.141592653589793D+00;
e* = 2.718281828459045D+00;
VAR Inf*, nInf*: LONGREAL;
PROCEDURE IsNan*(x: LONGREAL): BOOLEAN;
VAR h, l: SET;
BEGIN
sys.GET(sys.ADR(x), l);
sys.GET(sys.ADR(x) + 4, h);
RETURN (h * {20..30} = {20..30}) & ((h * {0..19} # {}) OR (l * {0..31} # {}))
END IsNan;
PROCEDURE IsInf*(x: LONGREAL): BOOLEAN;
RETURN ABS(x) = sys.INF(LONGREAL)
END IsInf;
PROCEDURE Max(A, B: LONGREAL): LONGREAL;
VAR Res: LONGREAL;
BEGIN
IF A > B THEN
Res := A
ELSE
Res := B
END
RETURN Res
END Max;
PROCEDURE Min(A, B: LONGREAL): LONGREAL;
VAR Res: LONGREAL;
BEGIN
IF A < B THEN
Res := A
ELSE
Res := B
END
RETURN Res
END Min;
PROCEDURE SameValue(A, B: LONGREAL): BOOLEAN;
VAR Epsilon: LONGREAL; Res: BOOLEAN;
BEGIN
Epsilon := Max(Min(ABS(A), ABS(B)) * 1.0D-12, 1.0D-12);
IF A > B THEN
Res := (A - B) <= Epsilon
ELSE
Res := (B - A) <= Epsilon
END
RETURN Res
END SameValue;
PROCEDURE IsZero(x: LONGREAL): BOOLEAN;
RETURN ABS(x) <= 1.0D-12
END IsZero;
PROCEDURE [stdcall] sqrt*(x: LONGREAL): LONGREAL;
BEGIN
sys.CODE("DD4508D9FAC9C20800")
RETURN 0.0D0
END sqrt;
PROCEDURE [stdcall] sin*(x: LONGREAL): LONGREAL;
BEGIN
sys.CODE("DD4508D9FEC9C20800")
RETURN 0.0D0
END sin;
PROCEDURE [stdcall] cos*(x: LONGREAL): LONGREAL;
BEGIN
sys.CODE("DD4508D9FFC9C20800")
RETURN 0.0D0
END cos;
PROCEDURE [stdcall] tan*(x: LONGREAL): LONGREAL;
BEGIN
sys.CODE("DD4508D9F2DEC9C9C20800")
RETURN 0.0D0
END tan;
PROCEDURE [stdcall] arctan2*(y, x: LONGREAL): LONGREAL;
BEGIN
sys.CODE("DD4508DD4510D9F3C9C21000")
RETURN 0.0D0
END arctan2;
PROCEDURE [stdcall] ln*(x: LONGREAL): LONGREAL;
BEGIN
sys.CODE("D9EDDD4508D9F1C9C20800")
RETURN 0.0D0
END ln;
PROCEDURE [stdcall] log*(base, x: LONGREAL): LONGREAL;
BEGIN
sys.CODE("D9E8DD4510D9F1D9E8DD4508D9F1DEF9C9C21000")
RETURN 0.0D0
END log;
PROCEDURE [stdcall] exp*(x: LONGREAL): LONGREAL;
BEGIN
sys.CODE("DD4508D9EADEC9D9C0D9FCDCE9D9C9D9F0D9E8DEC1D9FDDDD9C9C20800")
RETURN 0.0D0
END exp;
PROCEDURE [stdcall] round*(x: LONGREAL): LONGREAL;
BEGIN
sys.CODE("DD4508D97DF4D97DF666814DF60003D96DF6D9FCD96DF4C9C20800")
RETURN 0.0D0
END round;
PROCEDURE [stdcall] frac*(x: LONGREAL): LONGREAL;
BEGIN
sys.CODE("50DD4508D9C0D93C24D97C240266814C2402000FD96C2402D9FCD92C24DEE9C9C20800")
RETURN 0.0D0
END frac;
PROCEDURE arcsin*(x: LONGREAL): LONGREAL;
RETURN arctan2(x, sqrt(1.0D0 - x * x))
END arcsin;
PROCEDURE arccos*(x: LONGREAL): LONGREAL;
RETURN arctan2(sqrt(1.0D0 - x * x), x)
END arccos;
PROCEDURE arctan*(x: LONGREAL): LONGREAL;
RETURN arctan2(x, 1.0D0)
END arctan;
PROCEDURE sinh*(x: LONGREAL): LONGREAL;
VAR Res: LONGREAL;
BEGIN
IF IsZero(x) THEN
Res := 0.0D0
ELSE
Res := (exp(x) - exp(-x)) / 2.0D0
END
RETURN Res
END sinh;
PROCEDURE cosh*(x: LONGREAL): LONGREAL;
VAR Res: LONGREAL;
BEGIN
IF IsZero(x) THEN
Res := 1.0D0
ELSE
Res := (exp(x) + exp(-x)) / 2.0D0
END
RETURN Res
END cosh;
PROCEDURE tanh*(x: LONGREAL): LONGREAL;
VAR Res: LONGREAL;
BEGIN
IF IsZero(x) THEN
Res := 0.0D0
ELSE
Res := sinh(x) / cosh(x)
END
RETURN Res
END tanh;
PROCEDURE arcsinh*(x: LONGREAL): LONGREAL;
RETURN ln(x + sqrt((x * x) + 1.0D0))
END arcsinh;
PROCEDURE arccosh*(x: LONGREAL): LONGREAL;
RETURN ln(x + sqrt((x - 1.0D0) / (x + 1.0D0)) * (x + 1.0D0))
END arccosh;
PROCEDURE arctanh*(x: LONGREAL): LONGREAL;
VAR Res: LONGREAL;
BEGIN
IF SameValue(x, 1.0D0) THEN
Res := Inf
ELSIF SameValue(x, -1.0D0) THEN
Res := nInf
ELSE
Res := 0.5D0 * ln((1.0D0 + x) / (1.0D0 - x))
END
RETURN Res
END arctanh;
PROCEDURE floor*(x: LONGREAL): LONGREAL;
VAR f: LONGREAL;
BEGIN
f := frac(x);
x := x - f;
IF f < 0.0D0 THEN
x := x - 1.0D0
END
RETURN x
END floor;
PROCEDURE ceil*(x: LONGREAL): LONGREAL;
VAR f: LONGREAL;
BEGIN
f := frac(x);
x := x - f;
IF f > 0.0D0 THEN
x := x + 1.0D0
END
RETURN x
END ceil;
PROCEDURE power*(base, exponent: LONGREAL): LONGREAL;
VAR Res: LONGREAL;
BEGIN
IF exponent = 0.0D0 THEN
Res := 1.0D0
ELSIF (base = 0.0D0) & (exponent > 0.0D0) THEN
Res := 0.0D0
ELSE
Res := exp(exponent * ln(base))
END
RETURN Res
END power;
PROCEDURE sgn*(x: LONGREAL): INTEGER;
VAR Res: INTEGER;
BEGIN
IF x > 0.0D0 THEN
Res := 1
ELSIF x < 0.0D0 THEN
Res := -1
ELSE
Res := 0
END
RETURN Res
END sgn;
BEGIN
Inf := sys.INF(LONGREAL);
nInf := -sys.INF(LONGREAL)
END Math.