IMPLEMENTATION MODULE LMathLib;
(*------------------------------------------------------------------------*)
(* Stellt verschiedene mathematische Funktionen zur Verf"ugung *)
(* Provides different mathematical funktiones *)
(*------------------------------------------------------------------------*)
(* Letzte Bearbeitung: *)
(* *)
(* 93, MRi: Erstellen erste Version fuer JPI M2 *)
(* 16.04.15: MRi: Einfügen von MachEps *)
(* 31.05.15: MRi: Einfügen von MachEpsKonst in ExpM1 *)
(* 09.06.15: MRi: Einfügen von arctan2 *)
(* 31.07.15: MRi: Einfügen genauer Konstanten im Definitionmodul *)
(* 03.08.15: MRi: Einfügen von Pythag *)
(* 07.10.15, MRi: Verschieben von LnGamma,GammaU und Errf ins Module *)
(* SpezFunkt *)
(* 02.11.15, MRi: Pow10 hinzugefuegt (abgeleitet aus CardPot) *)
(* 13.12.15, MRi: GSchnitt eingefuegt *)
(* 04.02.16, MRi: MachEps mit <* NOOPTIMIZE + *> versehen. *)
(* 14.05.16, MRi: Prozedur Faktor hinzugefuegt *)
(* Prozeduren Bessel, Y01, BessJ0, H01 aus LMathLib her- *)
(* ausgenommen und Modul SpezFunkt3 erstellt *)
(* 16.06.16, MRi: lg in log10 umbenannt, weitere Konstanten eingefuegt *)
(* LTRUNC,LFLOAT eliminiert *)
(* 25.07.16, MRi: MachEpsConst in MachEps umbenannt, Prozedur MachEps *)
(* in CalcMachEps umbenannt *)
(* 12.08.16, MRi: MaxInt & MinInt eingefuegt *)
(* 08.09.16, MRi: IntPow hinzugefuegt (abgeleitet aus CardPot) *)
(* 17.11.16, MRi: Pow2 hinzugefuegt (abgeleitet aus CardPot) *)
(* 06.01.17, MRi: sinh & cosh fuer kleine x verbessert, EvalChebPoly *)
(* eingefuegt *)
(* 21.04.17, MRi: RationalApprox und gcd eingefuegt *)
(* 15.01.18, MRi: CalcPrimes eingefuegt *)
(* 03.03.18, MRi: Log2 eingefuegt (von P. Moylan) *)
(* 11.03.18, MRi: ceil & entier hierher verschoben *)
(* 26.07.18, MRi: SmallR4 eingefuegt *)
(*------------------------------------------------------------------------*)
(* Implementation : Michael Riedl *)
(* Licence : GNU Lesser General Public License (LGPL) *)
(*------------------------------------------------------------------------*)
(* $Id: LMathLib.mod,v 1.8 2018/03/11 11:20:00 mriedl Exp mriedl $ *)
IMPORT IEEE;
IMPORT TestReal;
FROM LongMath IMPORT exp,sqrt,ln,arctan,round;
FROM Errors IMPORT ErrOut;
PROCEDURE EvalChebPoly( x : LONGREAL;
VAR A : ARRAY OF LONGREAL;
n : CARDINAL ): LONGREAL;
(*----------------------------------------------------------------*)
(* Auswerten ein Chebyshev-Polynoms nach dem Clenshaw-Algorithmus *)
(*----------------------------------------------------------------*)
VAR b0,b1,b2 : LONGREAL;
i : CARDINAL;
BEGIN
b2 := 0.0;
b1 := 0.0;
b0 := 0.0;
x := 2.0*x;
FOR i:=n-1 TO 0 BY -1 DO
b2 := b1;
b1 := b0;
b0 := x*b1 - b2 + A[i];
END;
RETURN 0.5*(b0 - b2);
END EvalChebPoly;
PROCEDURE sqr(x : LONGREAL) : LONGREAL; BEGIN RETURN x*x; END sqr;
PROCEDURE Rest(x : LONGREAL) : LONGREAL;
BEGIN
IF (x < 0.0) THEN
RETURN x + VAL(LONGREAL,VAL(INTEGER,x));
ELSE
RETURN x - VAL(LONGREAL,VAL(INTEGER,x));
END;
END Rest;
PROCEDURE ceil(X : LONGREAL) : INTEGER;
VAR iret : INTEGER;
BEGIN
IF (X >= 0.0) THEN
iret := + VAL(INTEGER,TRUNC(ABS(X)));
ELSE
iret := - VAL(INTEGER,TRUNC(ABS(X)));
END;
(* iret := VAL(INTEGER,intpart(X)); *)
IF (LFLOAT(iret) < X) THEN INC(iret); END;
RETURN iret;
END ceil;
PROCEDURE entier(X : LONGREAL) : INTEGER;
VAR iret : INTEGER;
BEGIN
IF (X >= 0.0) THEN
iret := + VAL(INTEGER,TRUNC(ABS(X)));
ELSE
iret := - VAL(INTEGER,TRUNC(ABS(X)));
END;
(* iret := VAL(INTEGER,intpart(X)); *)
IF (LFLOAT(iret) > X) THEN DEC(iret); END;
RETURN iret;
END entier;
PROCEDURE CardPot(x : LONGREAL;
i : CARDINAL) : LONGREAL;
VAR z : LONGREAL;
BEGIN
z:=1.0;
LOOP
IF ODD(i) THEN z:=z*x END;
i:=i DIV 2;
IF (i = 0) THEN EXIT END;
x:=x*x;
END; (* LOOP *)
RETURN z;
END CardPot;
PROCEDURE IntPow(x : LONGREAL;
i : INTEGER) : LONGREAL;
VAR z : LONGREAL;
neg : BOOLEAN;
BEGIN
neg:= i < 0; i:=ABS(i);
z:=1.0;
LOOP
IF ODD(i) THEN z:=z*x END;
i:=i DIV 2;
IF (i = 0) THEN EXIT END;
x:=x*x;
END; (* LOOP *)
IF neg THEN z:= 1.0 / z; END;
RETURN z;
END IntPow;
PROCEDURE Pow10(pot : INTEGER) : LONGREAL;
VAR z,b : LONGREAL;
p : INTEGER;
BEGIN
b := 10.0; z := 1.0; p := ABS(pot);
LOOP
IF ODD(p) THEN z:=z*b END;
p:=p DIV 2;
IF (p = 0) THEN EXIT END;
b:=b*b;
END; (* LOOP *)
IF (pot < 0) THEN RETURN 1.0 / z; ELSE RETURN z; END;
END Pow10;
PROCEDURE log10(x : LONGREAL) : LONGREAL;
CONST ln10 = 2.302585092994045680;
BEGIN
RETURN ln(x)/ln10;
END log10;
PROCEDURE Pow2(pow : CARDINAL) : CARDINAL;
VAR z,b,p : CARDINAL;
BEGIN
b := 2; z := 1; p := pow;
LOOP
IF ODD(p) THEN z:=z*b END;
p:=p DIV 2;
IF (p = 0) THEN EXIT END;
b:=b*b;
END; (* LOOP *)
RETURN z;
END Pow2;
PROCEDURE Log2(N : CARDINAL) : CARDINAL;
(* This code is from Peter Moylan *)
VAR test, result: CARDINAL;
BEGIN
test := MAX(CARDINAL) DIV 2 + 1;
result := 31;
WHILE test > N DO
DEC (result); test := test DIV 2;
END;
RETURN result;
END Log2;
PROCEDURE arctan2(x,y : LONGREAL) : LONGREAL;
(* Maximale Abweichung von fdlibm.atan2 < MachEps, Testroutine *)
(* vorhanden. *)
CONST eps=0.0; (* Fusch *)
BEGIN
IF (y > 0.0) THEN
RETURN arctan(x / y);
ELSIF (y < 0.0) THEN
IF (ABS(x) > 0.0) THEN
RETURN 2.0*arctan((sqrt(x*x + y*y) - y) / x);
ELSE
RETURN arctan(x / y) - pi;
END;
ELSE (* y = 0.0 *)
IF (x > 0.0) THEN
RETURN pi / 2.0 + eps;
ELSIF (x < 0.0) THEN
RETURN -pi / 2.0 - eps;
ELSE
RETURN NAN;
END;
END;
END arctan2;
PROCEDURE sinh(x : LONGREAL) : LONGREAL;
TYPE ChebVek = ARRAY [0..8] OF LONGREAL;
CONST ChebKoeff = ChebVek{
0.173042194047179631675883846985E+00,
0.875942219227604771549002634544E-01,
0.107947777456713275024272706516E-02,
0.637484926075475048156855545659E-05,
0.220236640492305301591904955350E-07,
0.498794018041584931494262802066E-10,
0.797305355411573048133738827571E-13,
0.947315871307254445124531745434E-16,
0.869349205044812416919840906342E-19
};
VAR x2,y,sh : LONGREAL;
C : ChebVek;
BEGIN
IF (ABS(x) < 1.0) THEN
IF (ABS(x) <= 1.0E-08) THEN
x2:=x*x;
sh := x*(1.0 + x2/6.0 + x2*x2/120.0);
ELSE
C := ChebKoeff;
y := 2.0*x*x - 1.0;
sh := x + x*EvalChebPoly(y,C,9);
END;
ELSE
y:=exp(x);
sh := (y - 1.0 / y) / 2.0;
END;
RETURN sh;
END sinh;
PROCEDURE cosh(x : LONGREAL) : LONGREAL;
TYPE ChebVek = ARRAY [0..7] OF LONGREAL;
CONST ChebKoeff = ChebVek{
0.847409967933085972383472682702E-01,
0.706975331110557282636312722683E-03,
0.315232776534872632694769980081E-05,
0.874315531301413118696907944019E-08,
0.165355516210397415961496829471E-10,
0.226855895848535933261822520050E-13,
0.236057319781153495473072617928E-16,
0.192684134365813759102742156829E-19
};
VAR x2,y,ch : LONGREAL;
C : ChebVek;
BEGIN
IF (ABS(x) < 1.0) THEN
x2:=x*x;
IF (ABS(x) <= 1.0E-08) THEN
ch := 1.0 + 0.5*x2*(1.0 + x2/12.0);
ELSE
C := ChebKoeff;
y := 2.0*x2 - 1.0;
ch := 1.0 + 0.5*x2 + x2*x2*EvalChebPoly(y,C,8);
END;
ELSE
y := exp(ABS(x));
ch := (y + 1.0 / y) / 2.0;
END;
RETURN ch;
END cosh;
PROCEDURE tanh(x : LONGREAL) : LONGREAL;
BEGIN
RETURN sinh(x) / cosh(x);
END tanh;
PROCEDURE arcosh(x : LONGREAL) : LONGREAL;
BEGIN
RETURN ln(x + sqrt(x*x - 1.0));
END arcosh;
PROCEDURE arsinh(x : LONGREAL) : LONGREAL;
BEGIN
RETURN ln(x + sqrt(1.0 + x*x));
END arsinh;
PROCEDURE artanh(x : LONGREAL) : LONGREAL;
BEGIN
RETURN 0.5*ln((1.0 + x) / (1.0 - x));
END artanh;
PROCEDURE MinCard(i,j : CARDINAL) : CARDINAL;
BEGIN
IF (i <= j) THEN RETURN i ELSE RETURN j END;
END MinCard;
PROCEDURE MaxCard(i,j : CARDINAL) : CARDINAL;
BEGIN
IF (i >= j) THEN RETURN i ELSE RETURN j END;
END MaxCard;
PROCEDURE MinInt(i,j : INTEGER) : INTEGER;
BEGIN
IF (i <= j) THEN RETURN i ELSE RETURN j END;
END MinInt;
PROCEDURE MaxInt(i,j : INTEGER) : INTEGER;
BEGIN
IF (i >= j) THEN RETURN i ELSE RETURN j END;
END MaxInt;
(*
PROCEDURE Mantisse(x : LONGREAL) : LONGREAL;
VAR y : LONGREAL;
j : INTEGER;
BEGIN
y:=log10(x);
IF (y < 0.0) THEN
j := VAL(INTEGER,(ABS(y))
ELSE
j := - VAL(INTEGER,(ABS(y))
END;
RETURN x*IntPot(10.0,j);
END Mantisse;
*)
PROCEDURE Sign(x: LONGREAL): LONGREAL;
BEGIN
IF (x > 0.0) THEN
RETURN 1.0;
ELSIF (x < 0.0) THEN
RETURN -1.0;
ELSE
RETURN 0.0;
END;
END Sign;
PROCEDURE sign2(a,b : LONGREAL) : LONGREAL;
BEGIN
IF (b >= 0.0) THEN RETURN ABS(a); ELSE RETURN -ABS(a); END;
END sign2;
PROCEDURE fact(n : CARDINAL) : LONGREAL;
VAR p : CARDINAL;
Prod : LONGREAL;
BEGIN
CASE n OF
0 : RETURN 1.0000000000000000E+000;|
1 : RETURN 1.0000000000000000E+000;|
2 : RETURN 2.0000000000000000E+000;|
3 : RETURN 6.0000000000000000E+000;|
4 : RETURN 2.4000000000000000E+001;|
5 : RETURN 1.2000000000000000E+002;|
6 : RETURN 7.2000000000000000E+002;|
7 : RETURN 5.0400000000000000E+003;|
8 : RETURN 4.0320000000000000E+004;|
9 : RETURN 3.6288000000000000E+005;|
10 : RETURN 3.6288000000000000E+006;|
11 : RETURN 3.9916800000000000E+007;|
12 : RETURN 4.7900160000000000E+008;|
13 : RETURN 6.2270208000000000E+009;|
14 : RETURN 8.7178291200000000E+010;|
15 : RETURN 1.3076743680000000E+012;|
16 : RETURN 2.0922789888000000E+013;|
17 : RETURN 3.5568742809600000E+014;|
18 : RETURN 6.4023737057280000E+015;|
19 : RETURN 1.2164510040883200E+017;|
20 : RETURN 2.4329020081766400E+018;|
ELSE (* n > 20 *)
Prod:=2.4329020081766400E+018; (* 20! *)
FOR p:=21 TO n DO Prod:=Prod*VAL(LONGREAL,p); END;
RETURN Prod;
END;
END fact;
PROCEDURE BinKoeff(i,j : CARDINAL) : LONGREAL;
VAR k : CARDINAL;
Sum : LONGREAL;
BEGIN
IF (j > i) THEN
ErrOut('i > j (MathLib1.BinKoeff) !');
RETURN 1.0;
END;
Sum:=1.0;
FOR k:=i TO j+1 BY -1 DO Sum:=VAL(LONGREAL,k)*Sum; END;
RETURN Sum / fact(i-j);
END BinKoeff;
PROCEDURE Exp(x : LONGREAL) : LONGREAL;
VAR Sum,dy,fact,xi : LONGREAL;
i : CARDINAL;
BEGIN
Sum:=1.0 + x;
fact:=2.0; xi:=x*x; i:=2;
LOOP
dy := xi / fact;
Sum := Sum + dy;
IF (ABS(dy) < MachEps*ABS(Sum)) THEN EXIT END;
IF (i > 1200) THEN
ErrOut('Keine Konvergenz (MathLib1.exp) !');
EXIT;
END;
INC(i); fact:=fact*VAL(LONGREAL,i); xi:=xi*x;
END;
RETURN Sum;
END Exp;
PROCEDURE ExpM1(x : LONGREAL) : LONGREAL;
VAR xn,fact,S0,S1 : LONGREAL;
i : CARDINAL;
BEGIN
IF (x < 1.0E-20) THEN
S1:=x;
ELSE
S1:=x; fact:=2.0;
xn:=x*x; i:=2;
LOOP
S0:=S1;
S1:=S1 + xn/fact;
IF (ABS(S0 - S1) < MachEps*ABS(S1)) THEN EXIT END;
xn:=xn*x;
INC(i); fact:=fact*VAL(LONGREAL,i);
END;
END;
RETURN S1;
END ExpM1;
PROCEDURE Pythag(a,b : LONGREAL) : LONGREAL;
VAR p,r,s,t,u,su : LONGREAL;
BEGIN
IF ABS(a) >= ABS(b) THEN p:=ABS(a); ELSE p:=ABS(b); END;
IF (p # 0.0) THEN
IF ABS(a) >= ABS(b) THEN r:=ABS(b); ELSE r:=ABS(a); END;
r:=r / p; r:=r*r;
LOOP
t := 4.0 + r; IF (t = 4.0) THEN EXIT; END;
s := r / t;
u := 1.0 + 2.0*s;
p := u*p;
su := s / u; r := su*su * r;
END;
END;
RETURN p;
END Pythag;
PROCEDURE ggt(n,m : INTEGER) : INTEGER;
BEGIN
IF (n MOD m) = 0 THEN
RETURN m;
ELSE
RETURN ggt(m,(n MOD m));
END;
END ggt;
PROCEDURE kgv(i,j : INTEGER) : INTEGER;
BEGIN
RETURN i*j DIV ggt(i,j);
END kgv;
PROCEDURE gcd(m,n : INTEGER) : INTEGER;
VAR a,b,A : INTEGER;
BEGIN
a := ABS(m); b := ABS(n);
IF (a = 0) THEN
RETURN b;
ELSE
WHILE (b # 0) DO
A := b;
b := a REM b;
a := A;
END;
RETURN a;
END;
END gcd;
PROCEDURE RationalApprox( x : LONGREAL;
N : INTEGER;
VAR l,r : INTEGER);
(*---------------------------------------------------------------*)
(* If x = a1 + 1/(a2 + 1/(a3 + 1/(a4 + ...))) *)
(* then best approximation is found by truncating this series *)
(* (with some adjustments in the last term). *)
(* *)
(* Note the fraction can be recovered as the first column of the *)
(* matrix *)
(* *)
(* ( a1 1 ) ( a2 1 ) ( a3 1 ) ... *)
(* ( 1 0 ) ( 1 0 ) ( 1 0 ) *)
(* *)
(* Instead of keeping the sequence of continued fraction terms, *)
(* we just keep the last partial product of these matrices. *)
(* *)
(* David Eppstein / UC Irvine / 8 Aug 1993 *)
(* With corrections from Arno Formella, May 2008 *)
(*---------------------------------------------------------------*)
CONST HUGE = VAL(LONGREAL,MAX(REAL));
VAR t,ai,l2,r2 : INTEGER;
m00,m01,m10,m11 : INTEGER;
X0 : LONGREAL;
BEGIN
(* initialize matrix *)
m00:=1; m01:=0;
m10:=0; m11:=1;
X0:=x; ai := VAL(INTEGER,X0);
(* loop finding terms until denom gets too big *)
LOOP
IF ((m10*ai + m11) > N) THEN EXIT; END;
t := m00*ai + m01;
m01 := m00;
m00 := t;
t := m10*ai + m11;
m11 := m10;
m10 := t;
IF (x = LFLOAT(ai)) THEN EXIT; END; (* AF: division by zero *)
x := 1.0 / (x - LFLOAT(ai));
IF (x >= HUGE) THEN EXIT; END; (* AF: representation failure *)
ai := VAL(INTEGER,x);
END;
l:=m00; r:=m10;
(* now try other possibility *)
ai := (N - m11) DIV m10;
l2 := m00*ai + m01;
r2 := m10*ai + m11;
IF ABS(X0 - (LFLOAT(l ) / LFLOAT(r ))) >
ABS(X0 - (LFLOAT(l2) / LFLOAT(r2)))
THEN
l:=l2; r:=r2;
END;
r2:=gcd(l,r);
l:=l DIV r2;
r:=r DIV r2;
END RationalApprox;
PROCEDURE Faktor( N : INTEGER;
VAR Faktoren : ARRAY OF INTEGER;
VAR nFaktor : INTEGER);
VAR D,M,I : INTEGER;
IstNeg : BOOLEAN;
BEGIN
IstNeg := N < 0;
N:=ABS(N);
nFaktor:=0;
REPEAT (* Test if multiple of 2 *)
D := N MOD 2;
IF (D = 0) THEN
N:=N DIV 2;
Faktoren[nFaktor]:=2; INC(nFaktor);
END;
UNTIL (D > 0);
REPEAT (* Test if multiple of 3 *)
D := N MOD 3;
IF (D = 0) THEN
Faktoren[nFaktor]:=3; INC(nFaktor);
N:=N DIV 3;
END;
UNTIL (D > 0);
M:=round(sqrt(VAL(LONGREAL,N))) + 1;
I:=6;
WHILE (I < M + 1) DO
(* Test of divisors 6i-1 and 6i+1 up to sqrt(N) *)
(* Prime numbers are of the form 6i-1 or 6i+1 *)
REPEAT
D := N MOD (I - 1);
IF (D = 0) THEN
Faktoren[nFaktor]:=I - 1;
INC(nFaktor);
N:=N DIV (I - 1);
END;
UNTIL (D > 0);
REPEAT
D := N MOD (I + 1);
IF (D = 0) THEN
Faktoren[nFaktor]:=I + 1;
INC(nFaktor);
N:=N DIV (I + 1);
END;
UNTIL (D > 0);
I:=I + 6;
END;
IF (N > 1) THEN
Faktoren[nFaktor]:=N;
INC(nFaktor);
END;
IF IstNeg THEN Faktoren[0]:=-Faktoren[0]; END;
END Faktor;
PROCEDURE CalcPrimes(VAR Primes : ARRAY OF CARDINAL;
N : CARDINAL);
VAR i,j,p : CARDINAL;
BEGIN
IF (N-1 > HIGH(Primes)) THEN
ErrOut(" N zu gross (CalcPrimes)");
N := HIGH(Primes) + 1;
END;
Primes[0] := 1;
p := 1;
FOR i:=1 TO N-1 DO
REPEAT
INC(p);
j:=1;
WHILE (j < i) AND (p MOD Primes[j] # 0) DO INC(j); END;
UNTIL (j = i);
Primes[i] := p;
END;
END CalcPrimes;
(* <* NOOPTIMIZE + *> *)
PROCEDURE CalcMachEps() : LONGREAL;
VAR eps,EinsPlusEps : LONGREAL;
BEGIN
eps:=0.5; EinsPlusEps:=1.5;
WHILE (EinsPlusEps # 1.0) DO
eps:=0.5*eps;
EinsPlusEps:=1.0+eps;
END;
RETURN 2.0*eps;
END CalcMachEps;
PROCEDURE SmallR4() : LONGREAL;
VAR s,klein : REAL;
BEGIN
s:=1.0;
REPEAT
klein := s;
s := s / 16.0;
UNTIL (s = 0.0);
klein := 1.0*(klein / MachEpsR4);
RETURN VAL(LONGREAL,klein);
END SmallR4;
BEGIN
IF (MachEps # CalcMachEps()) THEN
ErrOut("Fehler in LMathLib CalcMachEps() # MachEps, bitte pruefen ");
END;
(*
INF := IEEE.INF;
NAN := IEEE.NAN;
*)
NAN := TestReal.Real8NaNquite();
INF := TestReal.Real8INF();
END LMathLib.