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.