--- a
+++ b/LMathLib.mod
@@ -0,0 +1,679 @@
+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                        *)
+  (*------------------------------------------------------------------------*)
+  (* 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;
+
+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.