Subversion Repositories Kolibri OS

(*

    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.