DEFINITION MODULE LMathLib;

  (*------------------------------------------------------------------------*)
  (* Stellt verschiedene mathematische Funktionen zur Verf"ugung            *)
  (* Provides different mathematical funktiones                             *)
  (*                                                                        *)
  (* Konstanten zum Teil aus:                                               *)
  (*                                                                        *)
  (*   Computer Approximations                                              *)
  (*   Hart, Cheney, e.a.                                                   *)
  (*   The SIAM Series in Applied Mathematics                               *)
  (*   John Wiley & Sons, INC. New York London Sydney, 1968                 *)
  (*------------------------------------------------------------------------*)
  (* Implementation : Michael Riedl                                         *)
  (* Licence        : GNU Lesser General Public License (LGPL)              *)
  (*------------------------------------------------------------------------*)

  (* $Id: LMathLib.def,v 1.11 2018/03/11 11:19:56 mriedl Exp mriedl $ *)

IMPORT IEEE;

CONST C         = 0.57721566490153286060651209008240243; (* Eulerkonstante. *)
      sqrtpi    = 1.772453850905516027298167483341145;   (* \sqrt{\pi} *)
      rezsqrtpi = 0.564189583547756286948079451560773;   (* \sqrt{1\over\pi} *)
      sqrt2     = 1.41421356237309504880168872420969808; (* 2^(1/2) *)
      sqrt05    = 0.70710678118654752440084436210484904; (* 2^(1/2)/2 *)
      sqrt3     = 1.73205080756887729352744634150587237; (* 3^(1/2) *)
      sqrt5     = 2.23606797749978969640917366873127624;

CONST pi        = 3.14159265358979323846264338327950288;
      twopi     = 6.28318530717958647692528676655900576; (* 2 \pi        *)
      halfpi    = 1.57079632679489661923132169163975144; (* \pi \over 2  *)
      quartpi   = 0.78539816339744830961566084581987572; (* \pi \over 4  *)
      sqrt2pi   = 2.50662827463100050241576528481105;    (* \sqrt{2\pi}  *)
      lnpi      = 1.144729885849400174143427351353059;   (* ln(pi)       *)
      sqrtpi05  = 0.88622693778351338858761284182448424; (* 0.5*sqrt(pi) *)
      e         = 2.71828182845904523536028747135266250;
      ln10      = 2.30258509299404568401799145468436421;
      ln2       = 0.69314718055994530941723212145817657;
      log2e     = 1.44269504088896340735992468100189213; (* log2(e)  *)

      GoldSec   = (1.0 + sqrt5) / 2.0; (* Goldener Schnitt, gold section *)
      CGold     = 2.0 - GoldSec;       (* (3 - sqrt(5))/2 = 0,381966...  *)

CONST MachEpsR8    =  2.22044604925031308E-16;
      MachEpsR4    =  1.192092895507813E-07;   (* 2^(-23) *)
      MachEps      =  MachEpsR8;
      Small        =  1.0 / VAL(LONGREAL,MAX(REAL));
      MinReal8     =  2.225073858507202E-308; (* 2^(-1022) *)
      MinReal4     =  1.175494350822288E-38;  (* 2^(-126)  *)
      MaxLog       =  709.7827128933840; (* Max. argument for Exp = Ln(MaxNum) *) 
      MinLog       = -708.3964185322641; (* Min. argument for Exp = Ln(MinNum) *) 
      MinExp       = -745.1332191019417; (* Number for which exp(MinExp) = 0   *) 
      MaxFac       =  170;               (* Max. argument for Factorial        *)

(*
 * CONST NAN          = IEEE.NAN;
 *       INF          = IEEE.INF;
 *)

VAR   INF          : LONGREAL;
      NAN          : LONGREAL;

  (*------------------------------------------------------------------------*)
  (* Funktion1 und FunktionN werden in Ausgleich & SpezFunk1 genutzt und    *)
  (* hier an zentraler Stelle definiert                                     *)
  (*------------------------------------------------------------------------*)

TYPE  Funktion1  = PROCEDURE(LONGREAL) : LONGREAL;
      FunktionN  = PROCEDURE(VAR ARRAY OF LONGREAL,
                                 CARDINAL) : LONGREAL;

(*--------------------------------------------------------------------------*)

PROCEDURE sqr(x : LONGREAL) : LONGREAL;

          (*----------------------------------------------------------------*)
          (* Analog zu Pascal sqr                                           *)
          (*----------------------------------------------------------------*)

PROCEDURE Rest(x : LONGREAL) : LONGREAL;

          (*----------------------------------------------------------------*)
          (* Berechet die Nachkommastellen von x                            *)
          (*----------------------------------------------------------------*)

PROCEDURE ceil(X : LONGREAL) : INTEGER;

          (*----------------------------------------------------------------*)
          (* Kleinsten ganzzahligen Wert, der nicht kleiner als x ist.      *)
          (* Bei negativen Zahlen werden also nur die Nachkommastellen      *)
          (* abgeschnitten, bei positiven Zahlen wird auf die naechst-      *)
          (* groessere ganze Zahl aufgerundet. Entspricht Rundung gegen     *)
          (* Plus Unendlich (Aufrundungsfunktion).                          *)
          (*                                                                *)
          (* Returns the least integer greater than or equal to X.          *)
          (*                                                                *)
          (* Beispiele / examples:                                          *)
          (*                                                                *)
          (*   ceil(-1.5) = -1.0, ceil(1.5) = 2.0                           *)
          (*   ceil(-1.0) = -1.0, ceil(1.0) = 1.0                           *)
          (*----------------------------------------------------------------*)

PROCEDURE entier(X : LONGREAL) : INTEGER;

          (*----------------------------------------------------------------*)
          (* Groessten ganzzahligen Wert, der nicht groesser als X ist.     *)
          (* Bei negativen Zahlen wird also auf die naechstkleinere ganze   *)
          (* Zahl abgerundet, bei positiven Zahlen werden nur die Nach-     *)
          (* kommastellen abgeschnitten.  Entspricht Rundung gegen Minus    *)
          (* Unendlich (Abrundungsfunktion).                                *)
          (*                                                                *)
          (* Returns the largest integer not grater than the value of X     *)
          (* ('rounding toward zero'). In other languages this is referred  *)
          (* to as the floor function.                                      *)
          (*                                                                *)
          (*   X - 1 < entier(X) <= X                                       *)
          (*                                                                *)
          (* Beispiele / examples:                                          *)
          (*                                                                *)
          (*   entier(-1.5) = -2.0, entier(1.5) = 1.0                       *)
          (*   entier(-1.0) = -1.0, entier(1.0) = 1.0                       *)
          (*                                                                *)
          (* This function is equal to the entier function defined in the   *)
          (* "Revised report on the algorithmic language Algol 60",         *)
          (* sec. 3.2.5.                                                    *)
          (*----------------------------------------------------------------*)

PROCEDURE CardPot(x : LONGREAL;
                  i : CARDINAL) : LONGREAL;

          (*----------------------------------------------------------------*)
          (* Berechet i-te Potenz von x / x to the power of i               *)
          (*----------------------------------------------------------------*)

PROCEDURE IntPow(x : LONGREAL;
                 i : INTEGER) : LONGREAL;

          (*----------------------------------------------------------------*)
          (* Berechet i-te Potenz von x / x to the integer power of i       *)
          (*----------------------------------------------------------------*)

PROCEDURE Pow10(pot : INTEGER) : LONGREAL;

          (*----------------------------------------------------------------*)
          (* Berechnet $ 10^{pot} $ = 10**pot / power of 10**pot            *)
          (*----------------------------------------------------------------*)

PROCEDURE log10(x : LONGREAL) : LONGREAL;

          (*----------------------------------------------------------------*)
          (* Berechnet den dekadischen Logaritmus von x (ln base 10 of x)   *)
          (* Calculates the logarithm of x to the base of 10                *)
          (*----------------------------------------------------------------*)

PROCEDURE Pow2(pow : CARDINAL) : CARDINAL;
          
          (*----------------------------------------------------------------*)
          (* Berechnet die "pow"-sche Potenz von 2 = 2**pow                 *)
          (* Two to the power of "pow"                                      *)
          (*----------------------------------------------------------------*)

PROCEDURE Log2(N : CARDINAL) : CARDINAL;

          (*----------------------------------------------------------------*)
          (* Logarithm to base 2 (rounded down to the next-lower integral   *)
          (* value, in case N is not a power of 2.)  The result is in the   *)
          (* range 0..31.                                                   *)
          (*                                                                *)
          (* This code is from Peter Moylan                                 *)
          (*----------------------------------------------------------------*)

PROCEDURE arctan2(x,y : LONGREAL) : LONGREAL;

          (*----------------------------------------------------------------*)
          (* Berechnet den Arcus Tangens von x / y.                         *)
          (*                                                                *)
          (* Das Ergebnis liegt im Interval [-pi,pii] und ist, je nach      *)
          (* nach Quadrant, vorzeichenbehaftet.                             *)
          (*                                                                *)
          (* Inverse tangent of x / y. Result is in range -pi to pi.        *)
          (*----------------------------------------------------------------*)

PROCEDURE sinh(x : LONGREAL) : LONGREAL;

          (*----------------------------------------------------------------*)
          (* Berechnet den Sinushyperbolikus von x im ersten und vierten    *)
          (* Quadranten, je nach Vorzeichen von x                           *)
          (*----------------------------------------------------------------*)

PROCEDURE cosh(x : LONGREAL) : LONGREAL;

          (*----------------------------------------------------------------*)
          (* Berechnet den Cosinushyperbolikus von x                        *)
          (*----------------------------------------------------------------*)

PROCEDURE tanh(x : LONGREAL) : LONGREAL;

          (*----------------------------------------------------------------*)
          (* Berechnet den Tangenshyperbolikus von x                        *)
          (*----------------------------------------------------------------*)

PROCEDURE arcosh(x : LONGREAL) : LONGREAL;  

          (*----------------------------------------------------------------*)
          (* Berechnet den Areacosinushyperbolikus von x                    *)
          (*----------------------------------------------------------------*)

PROCEDURE arsinh(x : LONGREAL) : LONGREAL;  

          (*----------------------------------------------------------------*)
          (* Berechnet den Areasinushyperbolikus von x                      *)
          (*----------------------------------------------------------------*)

PROCEDURE artanh(x : LONGREAL) : LONGREAL;  

          (*----------------------------------------------------------------*)
          (* Berechnet den Areatangenshyperbolikus von x                    *)
          (*----------------------------------------------------------------*)

PROCEDURE MinCard(i,j : CARDINAL) : CARDINAL;

          (*----------------------------------------------------------------*)
          (* Minimum von i und j                                            *)
          (*----------------------------------------------------------------*)

PROCEDURE MaxCard(i,j : CARDINAL) : CARDINAL;

          (*----------------------------------------------------------------*)
          (* Maximum von i und j                                            *)
          (*----------------------------------------------------------------*)

PROCEDURE MinInt(i,j : INTEGER) : INTEGER;

          (*----------------------------------------------------------------*)
          (* Minimum von i und j                                            *)
          (*----------------------------------------------------------------*)

PROCEDURE MaxInt(i,j : INTEGER) : INTEGER;

          (*----------------------------------------------------------------*)
          (* Maximum von i und j                                            *)
          (*----------------------------------------------------------------*)

PROCEDURE Sign(x: LONGREAL): LONGREAL;

          (*----------------------------------------------------------------*)
          (* Entsprechung der FORTRAN SIGN-Funktion                         *)
          (*----------------------------------------------------------------*)

PROCEDURE sign2(a,b : LONGREAL) : LONGREAL;

          (*----------------------------------------------------------------*)
          (* Gibt a mit dem Vorzeichen von b zurueck.                       *)
          (* Entspricht FORTRAN SIGN(a,b)                                   *)
          (*                                                                *)
          (* Returns a with the sign of b                                   *)
          (*----------------------------------------------------------------*)

PROCEDURE fact(n : CARDINAL) : LONGREAL;

          (*----------------------------------------------------------------*)
          (* Berechet die Fakult"at x! von x, n \in [0..MaxFac]             *)
          (*----------------------------------------------------------------*)

PROCEDURE BinKoeff(i,j : CARDINAL) : LONGREAL;

          (*----------------------------------------------------------------*)
          (* Berechnet i "uber j                                            *)
          (*----------------------------------------------------------------*)

PROCEDURE Exp(x : LONGREAL) : LONGREAL; 

          (*----------------------------------------------------------------*)
          (* Nur zu Testzecken exp(x)                                       *)
          (*----------------------------------------------------------------*)

PROCEDURE ExpM1(x : LONGREAL) : LONGREAL; 

          (*----------------------------------------------------------------*)
          (* Berechnet  exp(x) - 1.0  f"ur sehr kleine x                    *)
          (*----------------------------------------------------------------*)

PROCEDURE Pythag(a,b : LONGREAL) : LONGREAL;

          (*---------------------------------------------------------------*)
          (* Finds sqrt(a**2 + b**2) without overflow or destructive       *)
          (* underflow                                                     *)
          (*                                                               *)
          (* Modula-2 version of Eispack pythag subroutine                 *)
          (*---------------------------------------------------------------*)

PROCEDURE gcd(m,n : INTEGER) : INTEGER;

          (*---------------------------------------------------------------*)
          (* Iterative Greatest Common Divisor routine                     *)
          (*---------------------------------------------------------------*)

PROCEDURE ggt(n,m : INTEGER) : INTEGER;

          (*----------------------------------------------------------------*)
          (* Berechnet den gr"o\3ten gemeinsamen Teiler von n und m         *)
          (* Greatest common divisor of n and m                             *)
          (*----------------------------------------------------------------*)

PROCEDURE kgv(i,j : INTEGER) : INTEGER;

          (*----------------------------------------------------------------*)
          (* Berechnet das kleinstes gemeinsame Vielfache von i und j       *)
          (* Smallest common multiple of i and j                            *)
          (*----------------------------------------------------------------*)

PROCEDURE RationalApprox(    x   : LONGREAL;
                             N   : INTEGER;
                         VAR l,r : INTEGER);

          (*---------------------------------------------------------------*)
          (* Berechnet die rationale Approximation der Zahl x mit l,r < N  *)
          (*                                                               *)
          (*   x = \approx(l / r), l,r <=N                                 *)
          (*                                                               *)
          (* Find rational approximation to a given real number x based    *)
          (* on the theory of continued fractions                          *)
          (*---------------------------------------------------------------*)

PROCEDURE Faktor(    N        : INTEGER;
                 VAR Faktoren : ARRAY OF INTEGER;
                 VAR nFaktor  : INTEGER);

          (*----------------------------------------------------------------*)
          (* Faktorisiert die Zahl N und gibt die Faktoren im Feld Faktoren *)
          (* zurueck so dass gilt N = \PI_{i=0]^{nFaktor-1} Faktoren_i      *)
          (*                                                                *)
          (* Factors the numner N and returns the factors in array Faktoren *)
          (* such that N = \PI_{i=0]^{nFaktor-1} Faktoren_i                 *)
          (*----------------------------------------------------------------*)

PROCEDURE CalcPrimes(VAR Primes : ARRAY OF CARDINAL;
                         N      : CARDINAL);

          (*----------------------------------------------------------------*)
          (* Calculate the first N prime numbers using the algorithm of     *)
          (* Eratosthenes. Do not use for large numbers of N.               *)
          (*----------------------------------------------------------------*)

PROCEDURE CalcMachEps() : LONGREAL;

          (*----------------------------------------------------------------*)
          (* Berechne d. Maschinengenauigkeit / calculate machine precision *)
          (*----------------------------------------------------------------*)

PROCEDURE SmallR4() : LONGREAL;

          (*----------------------------------------------------------------*)
          (* Berechne eine kleine REAL4 - Zahl / calculate small real value *)
          (*----------------------------------------------------------------*)

END LMathLib.
