IMPLEMENTATION MODULE MyMathLibLong;
(*
  Created:   29.8.87
  Changed:   25.1.88/10.02.88/3.8.88/24.8.88/29.9.88 by 
             Stefan Salewski
             Stolper Weg 3
             2160 Stade   West-Germany
             Tel: 04141/61130
             
  Note: compiled with AMIGA Modula-2 System by AMSoft Verion from 5.5.88
   
  This Module may be freely copied. But please
  leave my name in. Thanks....Stefan 
*)    

  CONST
    PiH=        Pi*0.5;
    Euler=      2.71828182845904512; 
    einsLn10=   4.342944819032517E-1;    (*1.0/ln(10.0)*)
    MaxLongReal=MAX(LONGREAL);           (* 3. E+308 *)
    MinLongReal=MIN(LONGREAL);           (*-3. E+308 *)
    TwoPi=      6.28318530717958648;     (* 2.0*Pi *)
    OneOe=      1.0/Euler;               (* 0.368... *)
    EH10=       22026.4657948067434;     (* Euler^10 *)
    EH20=       485165195.409791529;     (* Euler^20 *)
    OneOeH10=   4.53999297624847755E-5;  (* OneOe^10 *)
    MaxInt2pi=32767.0*TwoPi;
    DegToRad=TwoPi/360.0;
    GonToRad=TwoPi/400.0;
    RadToDeg=360.0/TwoPi;
    RadToGon=400.0/TwoPi;    
(****************************************************************************)
  PROCEDURE MyUnit(w:LONGREAL):LONGREAL;
  (* rechnet Winkel in Grad oder Neugrad in Radiant um, wenn unit # rad     *)
  BEGIN
    IF unit=deg THEN
      RETURN w*DegToRad
    ELSIF unit=gon THEN
      RETURN w*GonToRad
    ELSE
      RETURN w
    END
  END MyUnit;
(****************************************************************************)
  PROCEDURE YourUnit(w:LONGREAL):LONGREAL;
  (* Rechnet Resultate von rad in die durch unit bestimmte Einheit um *)
  BEGIN
    IF unit=deg THEN
      RETURN w*RadToDeg
    ELSIF unit=gon THEN
      RETURN w*RadToGon
    ELSE
      RETURN w
    END
  END YourUnit;
(****************************************************************************)
  PROCEDURE ln(x:LONGREAL):LONGREAL;
  (* Natuerlicher Logarithmus zur Basis e=2.7182818 *)
    CONST Speichergroesse=30;
    VAR
      s,xm,xp,y:LONGREAL;
      speicher:ARRAY[0..Speichergroesse] OF LONGREAL;
      i,pos:CARDINAL;
  BEGIN
    IF x<=0.0 THEN (*Error*)
      errorNumber:=12;
      RETURN MinLongReal
    END;
    (*errorNumber:=0;*)
    IF x<OneOe THEN
      i:=0;
      WHILE x<1.0E-6 DO
        x:=x*EH10;
        INC(i,10)
      END;
      REPEAT
        x:=x*E;
        INC(i)
      UNTIL x>OneOe;
      RETURN ln(x)-LONGREAL(i)
    ELSIF x<=E THEN (* 0.368 < x < 2.72 *)
      i:=3;
      xm:=x-1.0;
      xp:=x+1.0;
      s:=xm/xp;
      y:=s*s;
      pos:=0;
      speicher[0]:=s+(y*s/3.0);
      s:=s*y;
      (*omi:=ABS(speicher[0])*epsilon;
        -0.1 < omi < 0.1  daher omi ~ 1 *)
      REPEAT
        INC(pos);
        INC(i,2);
        s:=s*y;
        speicher[pos]:=s/LONGREAL(i);
      (*UNTIL (ABS(speicher[pos])<=omi) OR (pos=Speichergroesse);*)
      UNTIL (ABS(speicher[pos])<=epsilon) OR (pos=Speichergroesse);
      s:=0.0;
      FOR i:=pos TO 0 BY (-1) DO
        s:=s+speicher[i]
      END;
      RETURN 2.0*s
    ELSE (* x>E *)
      i:=0;
      WHILE x>1.0E6 DO
        x:=x*OneOeH10;
        INC(i,10)
      END;
      REPEAT
        x:=x*OneOe;
        INC(i)
      UNTIL x<E;
      RETURN ln(x)+LONGREAL(i)
    END
  END ln;
(****************************************************************************)
  PROCEDURE exp(x:LONGREAL):LONGREAL;
  (* Exponentialfunktion *)
    CONST Speichergroesse=30;
    VAR
      s:LONGREAL;
      speicher:ARRAY[0..Speichergroesse] OF LONGREAL;
      i,pos:CARDINAL;
  BEGIN
    (*errorNumber:=0;*)
    IF x<=-710.0 THEN
      RETURN 0.0
    ELSIF x<0.0 THEN
      RETURN 1.0/exp(-x)
    ELSIF x<=1.0 THEN (* 0 < x < 1 *)
      i:=1;
      speicher[0]:=1.0;
      speicher[1]:=x;
      REPEAT
        INC(i);
        speicher[i]:=speicher[i-1]*x;
        speicher[i]:=speicher[i]/LONGREAL(i);
      UNTIL (speicher[i]<=epsilon) OR (i=Speichergroesse);
      s:=0.0;
      FOR pos:=i TO 0 BY (-1) DO
        s:=s+speicher[pos]
      END;
      RETURN s
    ELSIF x<=710.0 THEN
      pos:=CARDINAL(x); (* pos > 1 *)
      s:=exp(x-LONGREAL(pos));
      i:=pos;
      WHILE i>=20 DO
        s:=s*EH20;
        DEC(i,20);
      END;
      WHILE i>0 DO
        s:=s*E;
        DEC(i);
      END;
      RETURN s
    ELSE
      errorNumber:=73;
      RETURN MaxLongReal;
    END
  END exp;
(****************************************************************************)
  PROCEDURE sqrt(x:LONGREAL):LONGREAL;
    VAR
      ende:BOOLEAN;
      i:INTEGER;
      y,z,help:LONGREAL;
  BEGIN
    IF x<0.0 THEN
      errorNumber:=17;
      RETURN 0.0 (*ERROR*)
    ELSIF x=0.0 THEN
      (*errorNumber:=0;*)
      RETURN 0.0
    ELSE
      (*errorNumber:=0;*)
      IF x>=1.0 THEN
        IF x<1.3E154 THEN
          y:=2.0;
          REPEAT
            z:=y;
            y:=z*z;
          UNTIL y>=x;(*z > sqrt(x)*)
        ELSE
          z:=2.0E152
        END;
        y:=z;(* vergessen !!! *)
        REPEAT
          z:=y;
          y:=y*1.0E-4;
        UNTIL y*y<=x (*y < sqrt(x) *);
      ELSE
        IF x>7.0E-155 THEN
          y:=0.5;
          REPEAT
             z:=y;
             y:=z*z;
          UNTIL y<=x; (* z<= sqrt(x) *)
        ELSE
          z:=2.0E-152
        END;
        y:=z; (* vergessen !!! *)
        REPEAT
          z:=y;
          y:=y*1.0E+4
        UNTIL y*y>=x;(* y> sqrt(x)*)
      END;
      i:=0;
      ende:=FALSE;
      REPEAT
        INC(i);
        y:=z;
        z:=0.5*(y+(x/y));
        IF i>10 THEN
          help:=ABS(y-z);
          ende:=help <= (epsilon * z)
        END;
      UNTIL ende OR (i=99);
      RETURN z
    END
  END sqrt;
(****************************************************************************)
  PROCEDURE Mod2pi(VAR x:LONGREAL);
  BEGIN
    WHILE ABS(x)>MaxInt2pi DO
      IF x>0.0 THEN
        x:=x-MaxInt2pi
      ELSE
        x:=x+MaxInt2pi
      END;
    END;
    (*Fuer x > 0 und y > 0 gilt
      x=y (x DIV y) + x MOD y
      x MOD y = x - y (x DIV y)
              =x -y * CARDINAL(x / y)
    *)
    IF x>0.0 THEN
      x:=x-TwoPi*LONGREAL(CARDINAL(x/TwoPi))
    ELSE
      x:=x+TwoPi*LONGREAL(CARDINAL(-x/TwoPi))
    END
  END Mod2pi;
(****************************************************************************)
  PROCEDURE sinRad(x:LONGREAL):LONGREAL;
    CONST Speichergroesse=30;
    VAR
      coeff,s,sqrOfx,omi:LONGREAL;
      speicher:ARRAY[1..Speichergroesse] OF LONGREAL;
      i,j:INTEGER;
  BEGIN
    IF ABS(x) < 1.0E8 THEN
      (*errorNumber:=0;*)
      Mod2pi(x);
      (* -2*Pi < x < 2*Pi *)
      sqrOfx:=x*x;
      s:=x*sqrOfx;
      i:=1;
      coeff:=-1.0/6.0;
      speicher[1]:=x+s*coeff;
      omi:=ABS(speicher[1]*epsilon);
      REPEAT
        INC(i);
        coeff:=-coeff/LONGREAL((2*i)*(2*i+1));
        s:=s*sqrOfx;
        speicher[i]:=s*coeff;
      UNTIL (ABS(speicher[i]) <= omi) OR (i=Speichergroesse);
      s:=0.0;
      FOR j:= i TO 1 BY (-1) DO
        s:=s+speicher[j]
      END;
      RETURN s
    ELSE
      errorNumber:=18;
      RETURN 0.0
    END
  END sinRad;
  
  PROCEDURE sin(x:LONGREAL):LONGREAL;
  BEGIN
    RETURN sinRad(MyUnit(x))
  END sin;
(****************************************************************************)
  PROCEDURE cosRad(x:LONGREAL):LONGREAL;
    CONST Speichergroesse=30;
    VAR
      coeff,s,sqrOfx,omi:LONGREAL;
      speicher:ARRAY[1..Speichergroesse] OF LONGREAL;
      i,j:INTEGER;
  BEGIN
    IF ABS(x) < 1.0E8 THEN
      (*errorNumber:=0;*)
      Mod2pi(x);
      sqrOfx:=x*x;
      s:=sqrOfx;
      coeff:=-0.5;
      speicher[1]:=1.0+sqrOfx*coeff;
      omi:=ABS(speicher[1]*epsilon);
      i:=1;
      REPEAT
        INC(i);
        s:=s*sqrOfx;
        coeff:=-coeff/LONGREAL((2*i-1)*(2*i));
        speicher[i]:=s*coeff;
      UNTIL (ABS(speicher[i]) <= omi) OR (i=Speichergroesse);
      s:=0.0;
      FOR j:=i TO 1 BY (-1) DO 
        s:=s+speicher[j]
      END;
      RETURN s
    ELSE
      errorNumber:=18;
      RETURN 0.0
    END
  END cosRad;
  
  PROCEDURE cos(x:LONGREAL):LONGREAL;
  BEGIN
    RETURN cosRad(MyUnit(x))
  END cos;
(****************************************************************************)
  PROCEDURE arctanRad(x:LONGREAL):LONGREAL;
  (* Einheit des Resultats ist immer Bogenmass(Radiant) unabhaengig von unit *)
    CONST
      Speichergroesse=60;(* konvergiert langsam fuer x=1.5 *)
      PiMinusArctan2p1=2.0152155366959958; (* Pi - arctanRad(2.1*)
    VAR
      coeff,s,y,omi:LONGREAL;
      speicher:ARRAY[0..Speichergroesse] OF LONGREAL;
      i,j:CARDINAL;
      xNeg:BOOLEAN;
  BEGIN
    (*errorNumber:=0;*)
    xNeg:=x<0.0;
    x:=ABS(x);
    IF (x<1.5) AND (x>0.5) THEN
      s:=(x+2.1)/(1.0-(x*2.1));
      s:=arctanRad(s)+PiMinusArctan2p1;
      IF xNeg THEN
        RETURN -s
      ELSE
        RETURN s
      END;
    END;
    IF x <= 0.5 THEN (* 0 <x< 0.5 *)
      j:=0;
      i:=3;
      y:=x*x;
      s:=-x*y;
      speicher[0]:=x+s/3.0;
      omi:=ABS(speicher[0]*epsilon);
      REPEAT
        INC(j);
        INC(i,2);
        s:=-s*y;
        speicher[j]:=s/LONGREAL(i);
      UNTIL (ABS(speicher[j]) <= omi) OR (j=Speichergroesse)
    ELSE (* x>= 1.5*)
      s:=-1.0/x;
      y:=s*s;
      speicher[0]:=PiH+s;
      omi:=ABS(speicher[0]*epsilon);
      i:=1;
      j:=0;
      REPEAT
        INC(j);
        INC(i,2);
        s:=-s*y;
        speicher[j]:=s/LONGREAL(i);
      UNTIL (ABS(speicher[j]) <= omi)  OR (j=Speichergroesse)
    END;
    s:=0.0;
    FOR i:= j TO 0 BY (-1) DO
      s:=s+speicher[i]
    END;
    IF xNeg THEN
      RETURN -s
    ELSE
      RETURN s
    END
  END arctanRad;
  
  PROCEDURE arctan(x:LONGREAL):LONGREAL;
  (* Ist noetig da arctanRad rekursiv programmiert ist und daher das
     Resultat nicht in die gewuenschte Einheit umrechnen kann *)
  BEGIN
    RETURN YourUnit(arctanRad(x))
  END arctan;
(****************************************************************************)
  PROCEDURE neutraleFunc(x:LONGREAL):LONGREAL;
  BEGIN
    (*errorNumber:=0*)
    RETURN x
  END neutraleFunc;
(****************************************************************************)
  PROCEDURE abs(x:LONGREAL):LONGREAL;
  BEGIN
    (*errorNumber:=0*)
    RETURN ABS(x)
  END abs;
(****************************************************************************)
  PROCEDURE fac(x:LONGREAL):LONGREAL;
  (* Facultaet fuer ganze Zahlen 0 <= n <= 170 *)
    VAR
      j:[0..170];
      intx:INTEGER;
      z:LONGREAL;
      zuklein,zugross,istganz:BOOLEAN;
  BEGIN
    zugross:=x>170.0;
    zuklein:=x<0.0;
    IF (NOT zuklein) AND (NOT zugross) THEN
      intx:=INTEGER(x);
      istganz:=(x=LONGREAL(intx));
      IF istganz THEN
        (*errorNumber:=0*)
        z:=1.0;
        FOR j:=2 TO intx DO
          z:=z * LONGREAL(j)
        END;
        RETURN z
      ELSE
        errorNumber:=77;
        RETURN 0.0
      END
    ELSIF zugross THEN
      errorNumber:=75;
      RETURN 0.0
    ELSE
      errorNumber:=76;
      RETURN 0.0
    END
  END fac;
(****************************************************************************)
  PROCEDURE sqr(x:LONGREAL):LONGREAL;
  (* Quadrat *)
  BEGIN
    IF x<=1.0E154 THEN
      (*errorNumber:=0*)
      RETURN x*x
    ELSE errorNumber:=72;
      RETURN MaxLongReal
    END
  END sqr;
(****************************************************************************)
  PROCEDURE power(x,y:LONGREAL):LONGREAL;
  (*Raise x to the y th power  x^y *)
    CONST
      Epsilon=1.0E-16;
    VAR
      inty:INTEGER;
      j:CARDINAL;
      z:LONGREAL;
      expNegativ,ok:BOOLEAN;
  BEGIN
    (*errorNumber:=0*)
    ok:=(ABS(y)<20.0) AND (ABS(x)<=1.0E14);
    IF ok THEN
      IF y<0.0 THEN (* runden*)
        inty:=INTEGER(y-0.5)
      ELSE
        inty:=INTEGER(y+0.5)
      END;
    END;
    IF ok AND (ABS(y-LONGREAL(inty))<Epsilon) THEN
      expNegativ:= (inty<0);
      inty:=ABS(inty);
      z:=x;
      x:=1.0;
      FOR j:=1 TO inty DO
        x:=x*z
      END;
      IF expNegativ THEN
        IF x=0.0 THEN 
          errorNumber:=3
        ELSE 
          x:=1.0/x;
        END
      END
    ELSIF y=0.0 THEN
      x:=1.0
    ELSE
      IF x>0.0 THEN
        x:=exp(y*ln(x));
      ELSE
        x:=0.0;
        errorNumber:=4
      END
    END;
    RETURN x
  END power;
(****************************************************************************)
  PROCEDURE tan(x:LONGREAL):LONGREAL;
  (* tangens mit Fehlernummer *)
     VAR y:LONGREAL;
  BEGIN
    x:=MyUnit(x);
    y:=cosRad(x);
    IF y=0.0 THEN
      errorNumber:=5;
      RETURN MaxLongReal
    ELSE
      (*errorNumber:=0*)
      RETURN sinRad(x)/y
    END
  END tan;
(****************************************************************************)
  PROCEDURE cot(x:LONGREAL):LONGREAL;
  (* Kotangens  *)
    VAR y,z:LONGREAL;
  BEGIN
    x:=MyUnit(x);
    z:=cosRad(PiH-x);
    IF z=0.0 THEN
      errorNumber:=6;
      RETURN MaxLongReal
    ELSE 
      (*errorNumber:=0*)
      y:=sinRad(PiH-x)/z;
      RETURN y
    END
  END cot;
(****************************************************************************)
  PROCEDURE sec(x:LONGREAL):LONGREAL;
  (*Sekans = 1/cosRad(x) *)
    VAR y:LONGREAL;
  BEGIN
    x:=MyUnit(x);
    y:=cosRad(x);
    IF y=0.0 THEN
      errorNumber:=7;
      RETURN MaxLongReal
    ELSE 
      (*errorNumber:=0*)
      RETURN 1.0/y
    END
  END sec;
(****************************************************************************)
  PROCEDURE cosec(x:LONGREAL):LONGREAL;
  (* Kosekans =1/sinRad(x) *)
    VAR y:LONGREAL;
  BEGIN
    x:=MyUnit(x);
    y:=sinRad(x);
    IF y=0.0 THEN
      errorNumber:=8;
      RETURN MaxLongReal
    ELSE
      (*errorNumber:=0*)
      RETURN 1.0/sinRad(x)
    END
  END cosec;
(****************************************************************************)
  PROCEDURE arcsin(x:LONGREAL):LONGREAL;
  (* ArcusSinus= Umkehrfunktion des Sinus  -1<= x <= +1 *)
    VAR y:LONGREAL;
  BEGIN
    IF ABS(x)<=1.0 THEN
      (*errorNumber:=0*)
      IF x=1.0 THEN
        y:=PiH
      ELSIF x=-1.0 THEN
        y:=-PiH
      ELSE
        y:=arctanRad(x/sqrt(1.0-x*x))
      END;
      RETURN YourUnit(y)
    ELSE 
      errorNumber:=9;
      RETURN 0.0
    END
  END arcsin;
(****************************************************************************)
  PROCEDURE arccos(x:LONGREAL):LONGREAL;
  (* ArcusCosinus = Umkehrfunktion des Cosinus -1 <=x <= +1 *)
    VAR y:LONGREAL;
  BEGIN
    IF ABS(x)<=1.0 THEN
      (*errorNumber:=0*)
      y:=YourUnit(PiH)-arcsin(x);;
      RETURN y
    ELSE 
      errorNumber:=10;
      RETURN 0.0
    END
  END arccos;
(****************************************************************************)
  PROCEDURE arccot(x:LONGREAL):LONGREAL;
  (* ArcusKotangens = Umkehrfunktion des Kotangens *)
  BEGIN
    (*errorNumber:=0*)
    RETURN YourUnit(PiH-arctanRad(x))
  END arccot;
(****************************************************************************)
  PROCEDURE log(x:LONGREAL):LONGREAL;
  (* Logarithmus zur Basis 10 *)
  BEGIN
    IF x>0.0 THEN
      (*errorNumber:=0*)
      RETURN ln(x)*einsLn10
    ELSE
      errorNumber:=13;
      RETURN 0.0
    END
  END log;
(****************************************************************************)
  PROCEDURE sinh(x:LONGREAL):LONGREAL;
  (* Sinus Hyperbolicus  bzw. HyperbelSinus *)
    VAR y:LONGREAL;
  BEGIN
    (*errorNumber:=0*)
    y:=exp(x);
    y:=(y-1.0/y)*0.5;
    RETURN y
  END sinh;
(****************************************************************************)
  PROCEDURE cosh(x:LONGREAL):LONGREAL;
  (* Cosinus Hyperbolicus bzw. HyperbelCosinus *)
    VAR y:LONGREAL;
  BEGIN
    (*errorNumber:=0*)
    y:=exp(x);
    y:=y+1.0/y;
    RETURN y*0.5
  END cosh;
(****************************************************************************)
  PROCEDURE tanh(x:LONGREAL):LONGREAL;
  (* Tangens Hyperbolicus bzw. HyperbelTangens *)
    VAR 
      z,z1,y:LONGREAL;
  BEGIN
    (*errorNumber:=0*)
    z:=exp(x);
    IF z#0.0 THEN
      z1:=1.0/z;
      y:=(z-z1)/(z+z1);
      RETURN y
    ELSE
      RETURN 0.0
    END
  END tanh;
(****************************************************************************)
  PROCEDURE coth(x:LONGREAL):LONGREAL;
  (* Cotanges Hyperbolicus bzw. HyperbelCotangens *)
    VAR y,y1:LONGREAL;
  BEGIN
    IF x#0.0 THEN
      (*errorNumber:=0*)
      y:=exp(x);
      y1:=1.0/y;
      RETURN (y+y1)/(y-y1)
    ELSE
      errorNumber:=14;
      RETURN 0.0
    END
  END coth;
(****************************************************************************)
  PROCEDURE arsinh(x:LONGREAL):LONGREAL;
  (* AreaSinus = Umkehrfunktion von sinh(x) *)
    VAR y:LONGREAL;
  BEGIN
    (*errorNumber:=0*)
    y:=ln(x+sqrt(x*x+1.0));
    RETURN y
  END arsinh;
(****************************************************************************)
  PROCEDURE arcosh(x:LONGREAL):LONGREAL;
  (* AreaCosinus = Umkehrfunktion von cosh(x) *)
    VAR y:LONGREAL;
  BEGIN
    IF x>=1.0 THEN
      (*errorNumber:=0*)
      y:=ln(x+sqrt(x*x-1.0));
      RETURN y
    ELSE
      errorNumber:=15;
      RETURN 0.0
    END
  END arcosh;
(****************************************************************************)
  PROCEDURE artanh(x:LONGREAL):LONGREAL;
  (* AreaTangens = Umkehrfunktion tanh(x) *)
    VAR y:LONGREAL;
  BEGIN
    IF ABS(x)<1.0 THEN
      (*errorNumber:=0*)
      y:=0.5*ln((1.0+x)/(1.0-x));
      RETURN y
    ELSE
      errorNumber:=16;
      RETURN 0.0
    END
  END artanh;
(****************************************************************************)
  PROCEDURE arcoth(x:LONGREAL):LONGREAL;
  BEGIN
    IF ABS(x)>1.0 THEN
      (*errorNumber:=0*)
      RETURN 0.5*ln((x+1.0)/(x-1.0))
    ELSE
      errorNumber:=30;
      RETURN 0.0
    END
  END arcoth;
(****************************************************************************)
  PROCEDURE int(x:LONGREAL):LONGREAL;
  BEGIN
    IF ABS(x)<2147483648.0 THEN
      (*errorNumber:=0*)
      RETURN LONGREAL(LONGINT(x))
    ELSE
      errorNumber:=20;
      RETURN 0.0
    END
  END int;
(****************************************************************************)
BEGIN
(* 
 epsilon gibt die Genauigkeit an. Fuer eine groessere Geschwindigkeit
 mit weniger Genauigkeit kann epsilon bis ca. 1.0E-14 vergroessert werden
*)

  epsilon:=1.0E-20;
  unit:=rad;
  errorNumber:=0
END MyMathLibLong.mod