DEFINITION MODULE SpezFunkt2;

  (*------------------------------------------------------------------------*)
  (* Modul fuer einige spezielle Funktionen - Gammafunktion(en)             *)
  (* Module for the calculation of several forms of the Gamma function      *)
  (*------------------------------------------------------------------------*)
  (*------------------------------------------------------------------------*)
  (* Implementation : Michael Riedl                                         *)
  (* Licence        : GNU Lesser General Public License (LGPL)              *)
  (*------------------------------------------------------------------------*)


  (* $Id: SpezFunkt2.def,v 1.5 2018/05/25 13:35:18 mriedl Exp mriedl $ *)

CONST MaxGam =  171.624376956302;  (* Max. argument for Gamma   *)
      MaxLgm =    2.556348E+305;   (* Max. argument for LnGamma *)

PROCEDURE GammaIn(    X     : LONGREAL;
                      P     : LONGREAL;
                  VAR iFehl : INTEGER) : LONGREAL;
 
          (*----------------------------------------------------------------*)
          (* Berechnet die unvollständige Gammafunktion der Ordnung P       *)
          (* Computes the incomplete gamma ratio.                           *)
          (*                                                                *)
          (* A series expansion is used if P > X or X <= 1.  Otherwise, a   *)
          (* continued fraction approximation is used.                      *)
          (*                                                                *)
          (*   iFehl : 1 wenn P <= 0                                        *)
          (*           2 wenn X < 0                                         *)
          (*           4 wenn Fehler in lnGamma                             *)
          (*                                                                *)
          (* Ref.: Bhattacharjee, G; "Algorithm AS 32: The Incomplete Gamma *)
          (*       Integral", Applied Statistics, 19/3, 285-287 (1970)      *)
          (*                                                                *)
          (* Translation of a Fortran 90 version of Burkardt, John          *)
          (*                                                                *)
          (* This routine is less precise then GammaU (about 10-12 digits)  *)
          (* but considerably faster than its equivalent.                   *)
          (*----------------------------------------------------------------*)

PROCEDURE GammaU(a,x : LONGREAL) : LONGREAL;

          (*----------------------------------------------------------------*)
          (* Berechnet die regularisierte unvollständige Gammafunktion mit  *)
          (* der oberen Grenze x                                            *)
          (*                                                                *)
          (* P(a,x) = {1 \over \Gamma a} \int_0^x e^{-t} t^{a-1} dt         *)
          (*----------------------------------------------------------------*)

PROCEDURE IGamma(A,X : LONGREAL) : LONGREAL;

          (*----------------------------------------------------------------*)
          (* Incomplete Gamma function as defined by                        *)
          (*                                                                *)
          (*   IGamma(a,x) =  { 1 \over |a|} \int_0^x e^{-t} t^{a-1} dt     *)
          (*                                                                *)
          (* Based on dtmath86 / Cephes                                     *)
          (* This routine can be used as an alternative to GammaU           *)
          (*----------------------------------------------------------------*)

PROCEDURE JGamma(A,X : LONGREAL) : LONGREAL;

          (*----------------------------------------------------------------*)
          (* Complement of incomplete Gamma function as defined by          *)
          (*                                                                *)
          (*   JGamma(a,x) =  { 1 \over |a|} \int_x^\inf e^{-t} t^{a-1} dt  *)
          (*                                                                *)
          (* Based on dtmath86 / Cephes                                     *)
          (*----------------------------------------------------------------*)

PROCEDURE aLnGamma(    X     : LONGREAL;
                   VAR iFehl : INTEGER): LONGREAL;
 
          (*----------------------------------------------------------------*)
          (* Berechnet den natuerlichen Logarithmus der Gammafunktion       *)
          (*                                                                *)
          (*   iFehl : 1 wenn X <= 0                                        *)
          (*           2 wenn X > MaxX                                      *)
          (*                                                                *)
          (* Ref.: Macleod, Allan; "Algorithm AS 245, A Robust and Reliable *)
          (*       Algorithm for the Logarithm of the Gamma Function",      *)
          (*       Applied Statistics, 38/2, 397-402 (1989)                 *)
          (*                                                                *)
          (* Translation of a Fortran 90 version of Burkardt, John          *)
          (*                                                                *)
          (* This routine is less precise then dLnGamma (about              *)
          (* 10-12 digits) but considerably faster than its equivalent.     *)
          (*----------------------------------------------------------------*)
 
PROCEDURE LnGamma(z : LONGREAL) : LONGREAL;

          (*----------------------------------------------------------------*)
          (* Berechnet den Logarithmus der Gammafunktion mit                *)
          (* \Gamma z = \int_0^\infty t^{z-1} e^{-t} dt                     *)
          (* Genauigkeit  2.5E-13                                           *)
          (*----------------------------------------------------------------*)

PROCEDURE dLnGamma (X : LONGREAL) : LONGREAL;

          (*----------------------------------------------------------------*)
          (* This routine calculates the ln(Gamma) function for a           *)
          (* positive real argument X.                                      *)
          (* Computation is based on an algorithm outlined in               *)
          (* references 1 and 2. The program uses rational functions that   *)
          (* theoretically approximate ln(Gamma) to at least 18             *)
          (* significant decimal digits. The approximation for X > 12 is    *)
          (* from reference 3, while approximations for X < 12.0 are        *)
          (* similar to those in reference 1, but are unpublished. The      *)
          (* accuracy achieved depends on the arithmetic system, the        *)
          (* compiler, the intrinsic functions and proper selection of the  *)
          (* machine dependent constants in the implementation module.      *)
          (*                                                                *)
          (* Error returns:                                                 *)
          (*                                                                *)
          (*   The program returns the value INF when overflow would occur  *)
          (*   and NaN for negative arguments. The computation is believed  *)
          (*   to be free of underflow and overflow.                        *)
          (*                                                                *)
          (* References:                                                    *)
          (*                                                                *)
          (*  1) W. J. Cody and K. E. Hillstrom,                            *)
          (*     'Chebyshev Approximations for the Natural Logarithm of     *)
          (*     the Gamma Function,' Math. Comp. 21,1967, pp. 198-203.     *)
          (*  2) K. E. Hillstrom, ANL/AMD Program ANLC366S, DGAMMA/DLGAMA,  *)
          (*     May, 1969.                                                 *)
          (*  3) Hart, Et. Al., Computer Approximations, Wiley and sons,    *)
          (*     New York, 1968.                                            *)
          (*                                                                *)
          (*  Authors: Cody,W.J.; Stoltz,L. (Argonne National Laboratory)   *)
          (*----------------------------------------------------------------*)

PROCEDURE dGamma(x : LONGREAL) : LONGREAL;

          (*----------------------------------------------------------------*)
          (* Berechnet die vollst"andige Gammafunktion f"ur das Argument x  *)
          (*                                                                *)
          (* This routine calculates the Gamma function for a real          *)
          (* argument x.                                                    *)
          (*                                                                *)
          (* Lit: Broucke,R; Algorithm 446, CACM, 16, 254 (1973).           *)
          (*----------------------------------------------------------------*)

PROCEDURE Errf(x : LONGREAL) : LONGREAL;

          (*----------------------------------------------------------------*)
          (* Berechnet die Fehlerfunktion Errf(x) im Intervall [0,x)        *)
          (* Errf x = {2 \over \sqrt\pi} \int_0^x e^{-{t^2}} dt             *)
          (*----------------------------------------------------------------*)

PROCEDURE ErrorF(x : LONGREAL) : LONGREAL;

          (*----------------------------------------------------------------*)
          (* Berechnet die Fehlerfunktion durch Polynomapproximation        *)
          (*                                                                *)
          (* Lit.: Cody,W.J.; Waite,W.; "Software Manual for the Elementary *)
          (* Functions", Prentice-Hall, Englewood Cliffs, NJ (1980).        *)
          (*                                                                *)
          (* Basiert auf:                                                   *)
          (*                                                                *)
          (* http://history.dcs.ed.ac.uk/archive/os/emas/emas3/compilers/   *)
          (* ercs07/olddps.txt                                              *)
          (*                                                                *)
          (* Maximale Abweichung mit F77 Routine ist                        *)
          (*                                                                *)
          (* 4.4408921E-016 im Intervall (0 ,2 )                            *)
          (* 0.0            im Intervall (2 ,25)                            *)
          (*----------------------------------------------------------------*)

PROCEDURE CompErrorF(x : LONGREAL) : LONGREAL;

          (*----------------------------------------------------------------*)
          (* Berechnet die komplementaere Fehlerfunktion durch Polynom-     *)
          (* approximation.                                                 *)
          (*                                                                *)
          (* Lit. etc. siehe Kommentare für ApproxErf                       *)
          (*----------------------------------------------------------------*)

PROCEDURE CLnGamma(Z : LONGCOMPLEX) : LONGCOMPLEX;

          (*----------------------------------------------------------------*)
          (* Logarithm of the Gamma function for complex argument.          *)
          (* Returns CMPLX(INF,INF) in case of Z=0 or negative integer      *)
          (*                                                                *)
          (* Translation of CERN Libray function CLGAMA                     *)
          (*----------------------------------------------------------------*)

PROCEDURE CGamma(Z : LONGCOMPLEX) : LONGCOMPLEX;

          (*----------------------------------------------------------------*)
          (* Complex Gamma function                                         *)
          (* Returns CMPLX(INF,INF) in case of Z=0 or negative integer      *)
          (* or overflow                                                    *)
          (*----------------------------------------------------------------*)

PROCEDURE CdGamma(x : LONGCOMPLEX) : LONGCOMPLEX;

          (*----------------------------------------------------------------*)
          (* Complex gamma function in double precision                     *)
          (*----------------------------------------------------------------*)

PROCEDURE IBeta(A,B : LONGREAL;
                X   : LONGREAL) : LONGREAL;

          (*----------------------------------------------------------------*)
          (* Incomplete beta integral of the arguments, evaluated from zero *)
          (* to x.                                                          *)
          (*                                                                *)
          (*   IBeta(a,b,x) = { |a+b| \over |a| |b|}                        *)
          (*                  \int_0^x t^{a-1} (1-t)^{b-1} dt               *)
          (*                                                                *)
          (* Based on dtmath86 / Cephes                                     *)
          (*----------------------------------------------------------------*)

PROCEDURE SetGamitFast(opt : BOOLEAN);
 
          (*----------------------------------------------------------------*)
          (* Choose if aLnGamma / GammaIn are used in dGamit (speed) or     *)
          (* dLnGamma / GammaU (precision). If opt=true the faster variant  *)
          (* if opt=false the more precise subroutines are called.          *)
          (*----------------------------------------------------------------*)

PROCEDURE dGamit(A,X : LONGREAL): LONGREAL;

          (*----------------------------------------------------------------*)
          (* Calculate Tricomis form of the incomplete gamma function       *)
          (* defined by                                                     *)
          (*                                                                *)
          (*   gamit = { x^(-a) \over \Gamma a } \int_0^x e^{-t} t^{a - 1}  *)
          (*                                                                *)
          (* Code is a partial translation of W. Fullertons Fortran code    *)
          (* used in SLATEC.                                                *)
          (*                                                                *)
          (* Gautschi, W.; "A computational procedure for incomplete gamma  *)
          (* functions", ACM Transactions on mathematical software, 5,4     *)
          (* pp 466-481 (1979)                                              *)
          (* Gautschi, W.; "Incomplete gamma functions", Algorithm 542,     *)
          (* ACM Transactions on mathematical software 5,4; 482-489 (1979)  *)
          (* Code is a partial translation of W. Fullertons Fortran code    *)
          (* used in SLATEC.                                                *)
          (*----------------------------------------------------------------*)

END SpezFunkt2.
