SpezFunkt2.mod.m2pp
2063 lines (1826 with data), 72.2 kB
IMPLEMENTATION MODULE SpezFunkt2;
(*========================================================================*)
(* HINWEIS : Bitte nur die Datei SpezFunkt2.mod.m2pp editieren *)
(*========================================================================*)
(* Es sind 2 Versionen enthalten die mit *)
(* *)
(* m2pp -D __{Parameter}__ < SpezFunkt2.mod.m2pp > SpezFunkt2.mod *)
(* *)
(* mit Parameter = {XDS|GM2} erzeugt werden koennen. *)
(* *)
(* XDS : Der Export von Objekten aus lokalen Modulen wird mit *)
(* "IMPORT" angezeigt - das scheint nicht ganz "regelkonform" *)
(* GM2 : Normales "EXPORT" wird genutzt (ELSE-Fall). *)
(* *)
(* ansonsten gibt es keine Aenderungen am Quellcode *)
(*------------------------------------------------------------------------*)
(* Modul fuer einige Formen der Gammafunktion *)
(* Module for the calculation of several forms of the Gamma function *)
(*------------------------------------------------------------------------*)
(* IGamma,JGamm and IBeta are based on work found in dtmath086 Pascal *)
(* library with is in parts based on Cephes library (www.moshier.net) *)
(* Code for aLnGamma,GammaIn and dGamit is a partial translation of *)
(* W. Fullertons Fortran code used in SLATEC. *)
(* CLnGamma is a translation of CERN Library function GLGAMA *)
(*------------------------------------------------------------------------*)
(* Letzte Bearbeitung: *)
(* *)
(* 06.10.15, MRi: ErrorF und CompErrorF eingefuegt (erste Version) *)
(* 07.10.15, MRi: Verschieben von LnGamma,GammaU und Errf aus dem Modul *)
(* LMathLib *)
(* 08.10.15, MRi: dLnGamma & dGamma eingefuegt *)
(* 27.05.16, MRi: Erstellen der ersten Verison von GamSmall aud Basis von *)
(* Konstanten aus dtmath086 Bibliothek (Pascal). *)
(* 29.05.16, MRi: Hilfsprozedure D1Mach in Module F77func verschoben *)
(* GamSmall in dGamma eingebaut. *)
(* 4 Compilerwarnungen eleminiert durch 1. Ersetzen von *)
(* s1 in ErrorF durch B[17] resp. A[18] in Zeilen die mit *)
(* s = BJP2 bzw s = BJ gekennzeichent sind und 2. durch *)
(* Umwandlung der LOOP durch REPEAT-Schleifen in D9GamL *)
(* 27.08.16, MRi: Argument von dGamma ist nicht mehr "VAR" *)
(* 27.08.16, MRi: Erstellen der ersten Verionsn von IGamma,JGamma,IBeta *)
(* 16.02.17, MRi: Erstellen der ersten Version (aLnGamma,GamamaIn,dGamit) *)
(* 14.02.18, MRi: Hinzufuegen von CLnGamma und CGamma *)
(* 14.02.18, MRi: Erstellen der ersten Version von CdGamma *)
(* 23.05.18, MRi: Erweiterung von CdGamma und aufnahme in SpezFunkt2 *)
(*------------------------------------------------------------------------*)
(* Testroutinen vorhanden, IBeta gegen unabhaengigen Fortran-Code ge- *)
(* testet, IGamma & JGamma nur gegen Pascal-Aequivalent *)
(* Testcode (parially against Fortran versions) in TstSpzFun2 *)
(* dGamit wird in BoysInt2Lib intensiv genutzt *)
(*------------------------------------------------------------------------*)
(* Offene Punkte *)
(* *)
(* - Modulstruktur durchsehen und eventuell bereinigen. *)
(* - besserer Test fuer CLnGamma, CGamma *)
(* - Zusammenfuehren von SpezFunkt2 und SpezFunkt2b *)
(*------------------------------------------------------------------------*)
(* Implementation : Michael Riedl *)
(* Licence : GNU Lesser General Public License (LGPL) *)
(*------------------------------------------------------------------------*)
(* $Id: SpezFunkt2.mod.m2pp,v 1.4 2018/02/17 09:17:19 mriedl Exp $ *)
FROM Errors IMPORT FatalError,ErrOut,Fehler,Fehlerflag;
IMPORT LowLong;
FROM TestReal IMPORT Real8NaNquite,IsNearInteger;
IMPORT TestReal;
FROM LongMath IMPORT pi,sqrt,power,ln,exp,sin,cos,round;
FROM LMathLib IMPORT INF,NAN,MachEps,MinLog,MaxLog,sqrtpi,
sinh,cosh,tanh,arctan2;
IMPORT LongComplexMath;
FROM LongComplexMath IMPORT zero,one,conj,scalarMult;
FROM F77func IMPORT D1Mach;
VAR fast : BOOLEAN;
Sqrt2Pi : LONGREAL; (* sqrt(2 pi) *)
bot : LONGREAL; (* wird in d9Gamit genutzt *)
PROCEDURE sqr(x : LONGREAL) : LONGREAL; BEGIN RETURN x*x; END sqr;
PROCEDURE DSIGN(a,b : LONGREAL) : LONGREAL;
BEGIN
IF (b >= 0.0) THEN
RETURN ABS(a);
ELSE
RETURN - ABS(a);
END;
END DSIGN;
PROCEDURE GammaIn( X : LONGREAL; (* SLATEC *)
P : LONGREAL;
VAR iFehl : INTEGER) : LONGREAL;
CONST eps = 1.0E-08;
Big = 1.0E+37;
Small = 1.0E-37;
VAR a,an,arg,b,dif,factor : LONGREAL;
g,gin,rn,term,GInc : LONGREAL;
i : INTEGER;
pn : ARRAY [1..6] OF LONGREAL;
BEGIN
(* check the input. *)
IF (P <= 0.0) THEN
iFehl := 1;
RETURN 0.0;
ELSIF (X < 0.0) THEN
iFehl:=2;
RETURN 0.0;
ELSIF (X = 0.0) THEN
iFehl:=0;
RETURN 0.0;
END;
g := aLnGamma(P,iFehl);
IF (iFehl # 0) THEN
iFehl:=4;
RETURN 0.0;
END;
arg := P*ln(X) - X - g;
IF (arg < ln(Small)) THEN
iFehl:=3;
RETURN 0.0;
END;
iFehl:=0;
factor := exp(arg);
IF (X <= 1.0) OR (X < P) THEN (* calculation by series expansion. *)
gin := 1.0;
term := 1.0;
rn := P;
LOOP
rn := rn + 1.0;
term := term*X / rn;
gin := gin + term;
IF (term <= eps) THEN
EXIT;
END;
END; (* LOOP *)
GInc := gin*factor / P;
ELSE (* calculation by continued fraction. *)
a := 1.0 - P;
b := a + X + 1.0;
term := 0.0;
pn[1] := 1.0;
pn[2] := X;
pn[3] := X + 1.0;
pn[4] := X*b;
gin := pn[3] / pn[4];
LOOP
a:=a + 1.0;
b:=b + 2.0;
term:=term + 1.0;
an := a*term;
FOR i:=1 TO 2 DO
pn[i+4] := b*pn[i+2] - an*pn[i];
END;
IF (pn[6] # 0.0) THEN
rn := pn[5] / pn[6];
dif := ABS(gin - rn);
IF (dif <= eps) THEN (* absolute error tolerance satisfied *)
IF (dif <= eps*rn) THEN (* relative error tolerance satisfied *)
GInc := 1.0 - factor*gin;
EXIT;
END;
END;
gin := rn;
END;
FOR i:=1 TO 4 DO
pn[i] := pn[i+2];
END;
IF (ABS(pn[5]) >= Big) THEN
FOR i:=1 TO 4 DO
pn[i] := pn[i] / Big;
END;
END;
END; (* LOOP *)
END; (* IF *)
RETURN GInc;
END GammaIn;
PROCEDURE GammaU(a,x : LONGREAL) : LONGREAL;
(*----------------------------------------------------------------*)
(* Ben"otigt die PROCEDURE LnGamma *)
(*----------------------------------------------------------------*)
CONST MaxIter = 100;
eps = 2.5E-16;
VAR Iter : CARDINAL;
Sum,Delta,Fakt,I,Wert : LONGREAL;
g0,g1,b0,b1,a0,a1,ana,anf,ap : LONGREAL;
BEGIN
IF (x < 0.0) OR (a <= 0.0) THEN
FatalError('Argument (x < 0) bzw (a <= 0) (Funkt.GammaU) !)');
END;
IF (x = 0.0) THEN RETURN 0.0; END; (* !!! *)
Iter:=0;
IF (x < (a + 1.0)) THEN (* Reihenentwicklung *)
ap := a; Sum:=1.0 / a; Delta:=Sum;
REPEAT
INC(Iter);
ap := ap + 1.0; Delta:=Delta*x / ap; Sum:=Sum + Delta;
UNTIL (ABS(Delta) < eps*ABS(Sum)) OR (Iter > MaxIter);
Wert := Sum*exp(-x + a*ln(x) - LnGamma(a));
ELSE (* *)
g0:=0.0; a0:=1.0; a1:=x; b0:=0.0; b1:=1.0; Fakt:=1.0;
LOOP
INC(Iter); I:=LFLOAT(Iter);
ana := I - a;
a0 := (a1 + a0*ana)*Fakt;
b0 := (b1 + b0*ana)*Fakt;
anf:=I*Fakt;
a1:=x*a0 + anf*a1;
b1:=x*b0 + anf*b1;
IF (a1 # 0.0) THEN
Fakt := 1.0 / a1;
g1 := b1*Fakt;
IF (ABS((g1- g0) / g1) < eps) OR (Iter > MaxIter) THEN EXIT END;
g0 := g1;
END;
END; (* LOOP *)
Wert := 1.0 - g1*exp(-x + a*ln(x) - LnGamma(a));
END; (* IF x *)
IF (Iter > MaxIter) THEN
Fehler:=TRUE; Fehlerflag:='MaxIter "uberschritten (Funkt.GammaU) !';
ErrOut(Fehlerflag);
END;
RETURN Wert;
END GammaU;
PROCEDURE IGamma(A,X : LONGREAL) : LONGREAL;
CONST MinLn = 1.0 / MAX(LONGREAL);
VAR Ans,Ax,C,R : LONGREAL;
BEGIN
IF (X <= MinLn) OR (A <= 0.0) THEN
RETURN 0.0;
END;
IF (X > 1.0) AND (X > A) THEN
RETURN 1.0 - JGamma(A,X);
END;
(* Ax := A*ln(X) - X - LnGamma(A); *)
Ax := A*ln(X) - X - dLnGamma(A);
IF (Ax < MinLog) THEN
RETURN 0.0;
END;
Ax := exp(Ax);
R := A; C := 1.0; Ans := 1.0;
REPEAT (* Power series *)
R := R + 1.0;
C := C * X / R;
Ans := Ans + C;
UNTIL C / Ans <= MachEps;
RETURN Ans*Ax / A;
END IGamma;
PROCEDURE JGamma(A,X : LONGREAL) : LONGREAL;
CONST Big = 1.0 / MachEps;
VAR Ans,C,Yc,Ax,Y,Z,R,T : LONGREAL;
Pk,Pkm1,Pkm2,Qk,Qkm1,Qkm2 : LONGREAL;
BEGIN
IF (X <= 0.0) OR (A <= 0.0) THEN
RETURN 1.0;
END;
IF (X < 1.0) OR (X < A) THEN
RETURN 1.0 - IGamma(A,X);
END;
Ax := A*ln(X) - X - dLnGamma(A);
IF (Ax < MinLog) THEN
RETURN 0.0;
END;
Ax := exp(Ax);
(* Continued fraction *)
Y := 1.0 - A;
Z := X + Y + 1.0;
C := 0.0;
Pkm2 := 1.0;
Qkm2 := X;
Pkm1 := X + 1.0;
Qkm1 := Z * X;
Ans := Pkm1 / Qkm1;
REPEAT
C := C + 1.0;
Y := Y + 1.0;
Z := Z + 2.0;
Yc := Y * C;
Pk := Pkm1*Z - Pkm2*Yc;
Qk := Qkm1*Z - Qkm2*Yc;
IF (Qk # 0.0) THEN
R := Pk / Qk;
T := ABS((Ans - R) / R);
Ans := R;
ELSE
T := 1.0;
END;
Pkm2 := Pkm1;
Pkm1 := Pk;
Qkm2 := Qkm1;
Qkm1 := Qk;
IF (ABS(Pk) > Big) THEN
Pkm2:=Pkm2 * MachEps;
Pkm1:=Pkm1 * MachEps;
Qkm2:=Qkm2 * MachEps;
Qkm1:=Qkm1 * MachEps;
END;
UNTIL (T <= MachEps);
RETURN Ans*Ax;
END JGamma;
PROCEDURE aLnGamma( X : LONGREAL;
VAR iFehl : INTEGER): LONGREAL; (* SLATEC *)
CONST Hln2pi = 0.918938533204673; (* 0.5*ln(2.0*pi) *)
Xlow = 5.10E+05;
Xlarge = 1.0E+30;
TYPE Vec9 = ARRAY [1..9] OF LONGREAL;
Vec5 = ARRAY [1..5] OF LONGREAL;
CONST r1 = Vec9{
-2.66685511495E+00, -2.44387534237E+01,
-2.19698958928E+01, 1.11667541262E+01,
3.13060547623E+00, 6.07771387771E-01,
1.19400905721E+01, 3.14690115749E+01,
1.52346874070E+01
};
r2 = Vec9{
-7.83359299449E+01, -1.42046296688E+02,
1.37519416416E+02, 7.86994924154E+01,
4.16438922228E+00, 4.70668766060E+01,
3.13399215894E+02, 2.63505074721E+02,
4.33400022514E+01
};
r3 = Vec9{
-2.12159572323E+05, 2.30661510616E+05,
2.74647644705E+04, -4.02621119975E+04,
-2.29660729780E+03, -1.16328495004E+05,
-1.46025937511E+05, -2.42357409629E+04,
-5.70691009324E+02
};
r4 = Vec5{
0.2791953179185250, 0.4917317610505968,
0.0692910599291889, 3.3503438150223040,
6.0124592597641030
};
VAR x1,x2,y,LnG : LONGREAL;
BEGIN
IF (X <= 0.0) THEN (* check the input. *)
iFehl:=1;
RETURN 0.0;
ELSIF (X >= Xlarge) THEN
iFehl:=2;
RETURN 0.0;
END;
iFehl:=0;
(* calculation for 0 < x < 0.5 and 0.5 <= x < 1.5 combined. *)
IF (X < 1.5) THEN
IF (X < 0.5) THEN
LnG := - ln(X);;
y := X + 1.0;
(* test whether x < machine epsilon. *)
IF (y = 1.0) THEN
RETURN LnG;
END;
ELSE
LnG := 0.0;
y := X;
X := (X-0.5) - 0.5;
END;
LnG:=LnG +
X*((((r1[5]*y + r1[4])*y + r1[3])*y + r1[2])*y + r1[1]) /
((((y + r1[9])*y + r1[8])*y + r1[7])*y + r1[6]);
RETURN LnG;
END;
IF (X < 4.0) THEN
(* calculation for 1.5 <= x < 4.0. *)
y := X - 2.0;
LnG := y*((((r2[5]*X + r2[4])*X + r2[3])*X + r2[2])*X + r2[1]) /
((((X + r2[9])*X + r2[8])*X + r2[7])*X + r2[6]);
ELSIF (X < 12.0) THEN
(* calculation for 4.0 <= x < 12.0. *)
LnG := ((((r3[5]*X + r3[4])*X + r3[3])*X + r3[2])*X + r3[1]) /
(((( X + r3[9])*X + r3[8])*X + r3[7])*X + r3[6]);
ELSE
(* calculation for 12.0 <= x. *)
y := ln(X);
LnG := X*(y-1.0) - 0.5*y + Hln2pi;
IF (X <= Xlow) THEN
x1 := 1.0 / X;
x2 := x1*x1;
LnG := LnG + x1*((r4[3]*x2 + r4[2])*x2 + r4[1]) /
((x2 + r4[5])*x2 + r4[4]);
END;
END;
RETURN LnG;
END aLnGamma;
MODULE GAMMA; (* Stellt LnGamma zur Verfuegung *)
IMPORT pi,sqrt,ln,sin,Sqrt2Pi;
<* IF (__XDS__) THEN *>
IMPORT LnGamma; (* EXPORT *)
<* ELSE *>
EXPORT LnGamma;
<* END *>
VAR GamKoeff : ARRAY [0..15] OF LONGREAL;
PROCEDURE LnGamma(z : LONGREAL) : LONGREAL;
(*----------------------------------------------------------------*)
(* Lit.: Luke,Y.L.; "Special Functions and their Approximations", *)
(* Academic Press, New York, USA 1969. *)
(*----------------------------------------------------------------*)
VAR k : CARDINAL;
x,tmp,Sum,Z,N,Hk : LONGREAL;
BEGIN
IF (z = 1.0) THEN RETURN 0.0; END;
IF (z <= 1.0) THEN
x:=1.0 - z; tmp:=pi*x;
RETURN ln(tmp / (sin(tmp))) - LnGamma(x + 1.0); (* !!! *)
END;
x := z - 1.0;
tmp:=x + 5.5; tmp:= (x + 0.5)*ln(tmp) - tmp;
Sum:=0.0; Hk:=1.0; Z:=x; N:=z;
FOR k:=0 TO 15 DO
Sum:=Sum + GamKoeff[k] * Hk;
Hk := (Z / N) * Hk;
Z:=Z - 1.0; N:=N + 1.0;
END;
RETURN tmp + ln(Sqrt2Pi*Sum);
END LnGamma;
BEGIN
GamKoeff[ 0] := 41.624436916439068; GamKoeff[ 1] := -51.224241022374774;
GamKoeff[ 2] := 11.338755813488977; GamKoeff[ 3] := -0.747732687772388;
GamKoeff[ 4] := 0.008782877493061; GamKoeff[ 5] := -0.000001899030264;
GamKoeff[ 6] := 0.000000001946335; GamKoeff[ 7] := -0.000000000199345;
GamKoeff[ 8] := 0.000000000008433; GamKoeff[ 9] := 0.000000000001486;
GamKoeff[10] := -0.000000000000806; GamKoeff[11] := 0.000000000000293;
GamKoeff[12] := -0.000000000000102; GamKoeff[13] := 0.000000000000037;
GamKoeff[14] := -0.000000000000014; GamKoeff[15] := 0.000000000000006;
END GAMMA;
MODULE LibdLnGamma; (* Stellt dLnGamma zur Verfuegung *)
IMPORT TestReal; (* Modul *)
IMPORT ln,sqrt;
IMPORT MachEps;
<* IF (__XDS__) THEN *>
IMPORT dLnGamma; (* EXPORT *)
<* ELSE *>
EXPORT dLnGamma;
<* END *>
(*-----------------------------------------------------------------------*)
(* See also lot of hints given at http://dlmf.nist.gov/5 *)
(*-----------------------------------------------------------------------*)
CONST pnt68 = 0.6796875;
sqrtpi = 0.9189385332046727417803297; (* not sqrt(pi) !!! *)
(* Machine dependent parameters *)
xinf = MAX(LONGREAL);
xbig = 2.0*VAL(LONGREAL,MAX(REAL)); (* ??? *)
VAR frtbig : LONGREAL;
TYPE KoeffVek8 = ARRAY[1..8] OF LONGREAL;
KoeffVek7 = ARRAY[1..7] OF LONGREAL;
(*-------------------------------------------------------------*)
(* Numerator and denominator coefficients for rational minimax *)
(* approximation over (0.5,1.5). *)
(*-------------------------------------------------------------*)
CONST d1 = -5.772156649015328605195174E-01;
P1 = KoeffVek8{
4.945235359296727046734888E+00, 2.018112620856775083915565E+02,
2.290838373831346393026739E+03, 1.131967205903380828685045E+04,
2.855724635671635335736389E+04, 3.848496228443793359990269E+04,
2.637748787624195437963534E+04, 7.225813979700288197698961E+03
};
Q1 = KoeffVek8{
6.748212550303777196073036E+01, 1.113332393857199323513008E+03,
7.738757056935398733233834E+03, 2.763987074403340708898585E+04,
5.499310206226157329794414E+04, 6.161122180066002127833352E+04,
3.635127591501940507276287E+04, 8.785536302431013170870835E+03
};
(*-------------------------------------------------------------*)
(* Numerator and denominator coefficients for rational minimax *)
(* approximation over (1.5,4.0). *)
(*-------------------------------------------------------------*)
CONST d2 = 4.227843350984671393993777E-01;
P2 = KoeffVek8{
4.974607845568932035012064E+00, 5.424138599891070494101986E+02,
1.550693864978364947665077E+04, 1.847932904445632425417223E+05,
1.088204769468828767498470E+06, 3.338152967987029735917223E+06,
5.106661678927352456275255E+06, 3.074109054850539556250927E+06
};
Q2 = KoeffVek8{
1.830328399370592604055942E+02, 7.765049321445005871323047E+03,
1.331903827966074194402448E+05, 1.136705821321969608938755E+06,
5.267964117437946917577538E+06, 1.346701454311101692290052E+07,
1.782736530353274213975932E+07, 9.533095591844353613395747E+06
};
(*-------------------------------------------------------------*)
(* Numerator and denominator coefficients for rational minimax *)
(* approximation over (4.0,12.0). *)
(*-------------------------------------------------------------*)
CONST d4 = 1.791759469228055000094023;
P4 = KoeffVek8{
1.474502166059939948905062E+04, 2.426813369486704502836312E+06,
1.214755574045093227939592E+08, 2.663432449630976949898078E+09,
2.940378956634553899906876E+10, 1.702665737765398868392998E+11,
4.926125793377430887588120E+11, 5.606251856223951465078242E+11
};
Q4 = KoeffVek8{
2.690530175870899333379843E+03, 6.393885654300092398984238E+05,
4.135599930241388052042842E+07, 1.120872109616147941376570E+09,
1.488613728678813811542398E+10, 1.016803586272438228077304E+11,
3.417476345507377132798597E+11, 4.463158187419713286462081E+11
};
(*---------------------------------------------------------*)
(* Coefficients for minimax approximation over (12, INF). *)
(*---------------------------------------------------------*)
CONST C = KoeffVek7 {
-1.910444077728E-03, 8.4171387781295E-04,
-5.952379913043012E-04, 7.93650793500350248E-04,
-2.777777777777681622553E-03, 8.333333333333333331554247E-02,
5.7083835261E-03
};
PROCEDURE dLnGamma (X : LONGREAL) : LONGREAL;
(*----------------------------------------------------------------*)
(* Explanation of machine-dependent constants *)
(* *)
(* beta - radix for the floating-point representation *)
(* maxexp - the smallest positive power of beta that overflows *)
(* XBig - largest argument for which LN(GAMMA(X)) is *)
(* representable in the machine, i.e., the solution *)
(* to the equation ln(GAMMA(XBIG)) = beta**maxexp *)
(* XINF - largest machine representable floating-point *)
(* number; approximately beta**maxexp. *)
(* eps - The smallest positive floating-point number such *)
(* that 1.0 + eps > 1.0 *)
(* FRTBIG - Rough estimate of the fourth root of XBIG *)
(* *)
(* Latest modification: June 16, 1988 *)
(* *)
(* Adapted to Modula-2: Riedl,Michael 03.05.2015 *)
(*----------------------------------------------------------------*)
CONST eps = MachEps;
VAR i : CARDINAL;
res,corr : LONGREAL;
xden,xm1,xm2,xm4,xnum,Xsq : LONGREAL;
BEGIN
IF (X > 0.0) AND (X <= xbig) THEN
IF (X <= eps) THEN
res := -ln(X);
ELSIF (X <= 1.5) THEN (* eps < X <= 1.5 *)
IF (X < pnt68) THEN
corr:= -ln(X);
xm1 := X;
ELSE
corr := 0.0;
xm1 := (X - 0.5) - 0.5; (* Wieso nicht X - 1.0 ??? *)
END;
IF (X <= 0.5) OR (X >= pnt68) THEN
xden := 1.0;
xnum := 0.0;
FOR i:=1 TO 8 DO
xnum:=xnum*xm1 + P1[i];
xden:=xden*xm1 + Q1[i];
(*100:*) END;
res := corr + (xm1*(d1 + xm1*(xnum/xden)));
ELSE
xm2 := (X - 0.5) - 0.5;
xden := 1.0;
xnum := 0.0;
FOR i:=1 TO 8 DO
xnum:=xnum*xm2 + P2[i];
xden:=xden*xm2 + Q2[i];
(*200:*) END;
res := corr+xm2*(d2+xm2*(xnum/xden));
END;
ELSIF (X <= 4.0) THEN (* 1.5 < X <= 4.0 *)
xm2 := X - 2.0;
xden := 1.0;
xnum := 0.0;
FOR i:=1 TO 8 DO
xnum:=xnum*xm2 + P2[i];
xden:=xden*xm2 + Q2[i];
(*300:*) END;
res := xm2*(d2+xm2*(xnum/xden))
ELSIF (X <= 12.0) THEN (* 4 < X <= 12 *)
xm4 := X - 4.0;
xden := -1.0;
xnum := 0.0;
FOR i:=1 TO 8 DO
xnum:=xnum*xm4 + P4[i];
xden:=xden*xm4 + Q4[i];
(*400:*) END;
res := d4+xm4*(xnum/xden)
ELSE (* X >= 12.0 *)
res:=0.0;
IF (X <= frtbig) THEN
res := C[7];
Xsq := X*X;
FOR i:=1 TO 6 DO
res:=res / Xsq + C[i];
END;
END;
res := res / X;
corr := ln(X);
res:=res + sqrtpi - 0.5*corr;
res:=res + X*(corr - 1.0);
END;
ELSIF (X <= 0.0) THEN (* Return for bad arguments *)
res:=TestReal.Real8NaNquite();
ELSE (* this means that X > XBIG and so that log(gamma) overflows *)
res:=TestReal.Real8INF();
END;
RETURN res;
END dLnGamma;
BEGIN
frtbig := sqrt(sqrt(xbig));
END LibdLnGamma;
PROCEDURE GamSmall(X : LONGREAL) : LONGREAL;
(*----------------------------------------------------------------*)
(* Gamma function for small ABS(X) in the Range (0,0.03125] *)
(* *)
(* Ref.: Moshier, Stephen; "Methods and Programs for Mathematical *)
(* Functions", Prentice-Hall, Englewood Cliffs, NY (US) (1989) *)
(*----------------------------------------------------------------*)
VAR P : LONGREAL;
BEGIN
(* IF (X = 0.0) THEN
* (* Singulaer *)
* RETURN MAX(LONGREAL);
* ELSIF (ABS(X) > 0.03125) THEN
* (* Auserhalb des Wertebereichs *)
* RETURN MAX(LONGREAL);
* ELS*)
IF (X > 0.0) THEN (* peak relative error 4.2E-23 *)
P := ((((((((-1.193945051381510095614E-03 *X +
7.220599478036909672331E-03)*X -
9.622023360406271645744E-03)*X -
4.219773360705915470089E-02)*X +
1.665386113720805206758E-01)*X -
4.200263503403344054473E-02)*X -
6.558780715202540684668E-01)*X +
5.772156649015328608253E-01)*X + 1.0);
ELSE
X := - X; (* Peak relative error 5.16E-23 *)
P := (((((((( 1.133374167243894382010E-03 *X +
7.220837261893170325704E-03)*X +
9.621911155035976733706E-03)*X -
4.219773343731191721664E-02)*X -
1.665386113944413519335E-01)*X -
4.200263503402112910504E-02)*X +
6.558780715202536547116E-01)*X +
5.772156649015328608727E-01)*X - 1.0);
END;
RETURN 1.0 / (X * P);
END GamSmall;
MODULE LibdGamma; (* Stellt dGamma zur Verfuegung *)
IMPORT LowLong;
IMPORT TestReal;
IMPORT MachEps;
IMPORT D1Mach;
IMPORT pi,ln,exp,sqrt,sin,power,round;
IMPORT sqr,INF;
IMPORT ErrOut,FatalError;
IMPORT GamSmall;
<* IF (__XDS__) THEN *>
IMPORT dGamma; (* EXPORT *)
<* ELSE *>
EXPORT dGamma;
<* END *>
PROCEDURE DMAX1(x,y:LONGREAL):LONGREAL;
BEGIN
IF (x > y) THEN RETURN x ELSE RETURN y; END;
END DMAX1;
PROCEDURE DMIN(x,y:LONGREAL):LONGREAL;
BEGIN
IF (x > y) THEN RETURN y; ELSE RETURN x; END;
END DMIN;
PROCEDURE DINT(x : LONGREAL) : LONGREAL;
BEGIN
RETURN LowLong.intpart(x);
END DINT;
CONST sq2pil = 0.91893853320467274178032973640562; (* ln(sqrt(2*pi)) *)
TYPE Vektor15 = ARRAY [1..15] OF LONGREAL;
VAR Algmcs : Vektor15;
CONST algmcs = Vektor15{
+0.1666389480451863247205729650822E+00,
-0.1384948176067563840732986059135E-04,
+0.9810825646924729426157171547487E-08,
-0.1809129475572494194263306266719E-10,
+0.6221098041892605227126015543416E-13,
-0.3399615005417721944303330599666E-15,
+0.2683181998482698748957538846666E-17,
-0.2868042435334643284144622399999E-19,
+0.3962837061046434803679306666666E-21,
-0.6831888753985766870111999999999E-23,
+0.1429227355942498147573333333333E-24,
-0.3547598158101070547199999999999E-26,
+0.1025680058010470912000000000000E-27,
-0.3401102254316748799999999999999E-29,
+0.1276642195630062933333333333333E-30
};
TYPE Vektor42 = ARRAY [1..42] OF LONGREAL;
VAR Gamcs : Vektor42;
CONST gamcs = Vektor42{
+0.8571195590989331421920062399942E-02,
+0.4415381324841006757191315771652E-02,
+0.5685043681599363378632664588789E-01,
-0.4219835396418560501012500186624E-02,
+0.1326808181212460220584006796352E-02,
-0.1893024529798880432523947023886E-03,
+0.3606925327441245256578082217225E-04,
-0.6056761904460864218485548290365E-05,
+0.1055829546302283344731823509093E-05,
-0.1811967365542384048291855891166E-06,
+0.3117724964715322277790254593169E-07,
-0.5354219639019687140874081024347E-08,
+0.9193275519859588946887786825940E-09,
-0.1577941280288339761767423273953E-09,
+0.2707980622934954543266540433089E-10,
-0.4646818653825730144081661058933E-11,
+0.7973350192007419656460767175359E-12,
-0.1368078209830916025799499172309E-12,
+0.2347319486563800657233471771688E-13,
-0.4027432614949066932766570534699E-14,
+0.6910051747372100912138336975257E-15,
-0.1185584500221992907052387126192E-15,
+0.2034148542496373955201026051932E-16,
-0.3490054341717405849274012949108E-17,
+0.5987993856485305567135051066026E-18,
-0.1027378057872228074490069778431E-18,
+0.1762702816060529824942759660748E-19,
-0.3024320653735306260958772112042E-20,
+0.5188914660218397839717833550506E-21,
-0.8902770842456576692449251601066E-22,
+0.1527474068493342602274596891306E-22,
-0.2620731256187362900257328332799E-23,
+0.4496464047830538670331046570666E-24,
-0.7714712731336877911703901525333E-25,
+0.1323635453126044036486572714666E-25,
-0.2270999412942928816702313813333E-26,
+0.3896418998003991449320816639999E-27,
-0.6685198115125953327792127999999E-28,
+0.1146998663140024384347613866666E-28,
-0.1967938586345134677295103999999E-29,
+0.3376448816585338090334890666666E-30,
-0.5793070335782135784625493333333E-31
};
VAR ngam : INTEGER; (* Global für dGamma *)
dxrel,xmax,xmin,xsml : LONGREAL;
VAR nalgm : INTEGER; (* Global für D9lnMC *)
xbig,D9lnMCxmax : LONGREAL;
PROCEDURE dGamma(X: LONGREAL): LONGREAL;
PROCEDURE Dcsevl( x : LONGREAL;
VAR A : ARRAY OF LONGREAL;
n : INTEGER) : LONGREAL;
(*-------------------------------------------------------------*)
(* Getestet im Interval (-1,1) fuer beide Koeff.-vektoren *)
(* *)
(* Evaluate the n-term Chebyshev series a at x. *)
(* *)
(* Lit: Broucke,R; Algorithm 446, c.a.c.m., 16, 254 (1973). *)
(* *)
(* Input arguments *)
(* *)
(* x value at which the series is to be evaluated. *)
(* A array of n terms of a chebyshev series. in *)
(* evaluating a, only half the first coef is summed. *)
(* n number of terms in array A. *)
(*-------------------------------------------------------------*)
CONST eps = 1.0E-10;
VAR twox,b0,b1,b2 : LONGREAL;
i : INTEGER;
BEGIN
IF (n < 1) THEN
FatalError('SpezFunkt2.dGamma.Dcsevl: number of terms <= 0');
ELSIF (n > 1000) THEN
FatalError('SpezFunkt2.sGamma.Dcsevl: number of terms > 1000');
END;
IF (x < -(1.0+eps)) OR (x > (1.0+eps)) THEN
ErrOut('SpezFunkt2.dGamma.Dcsevl: x outside (-1,+1)');
END;
twox := 2.0*x;
b1:=0.0; b0:=0.0;
b2:=0.0; (* Wg. Compilerwarnung *)
FOR i:=1 TO n DO
b2 := b1;
b1 := b0;
b0 := twox*b1 - b2 + A[n-i];
END;
RETURN 0.5*(b0-b2);
END Dcsevl;
PROCEDURE D9lnMC(VAR x: LONGREAL): LONGREAL;
(*-------------------------------------------------------------*)
(* Compute the log gamma correction factor for x > 10 so that *)
(* ln (dGamma(x)) = ln(sqrt(2*pi)) + (x-.5)*ln(x) *)
(* - x + D9lnMC(x) *)
(* *)
(* Series for ALGM on the interval 0. to 1.00000E-02 *)
(* with weighted error 1.28E-31 *)
(* ln weighted error 30.89 *)
(* significant figures required 29.81 *)
(* decimal places required 31.48 *)
(*-------------------------------------------------------------*)
BEGIN
IF (x < 10.0) THEN
ErrOut('SpezFunkt2.dGamma.D9lnMC: x must be > 10');
END;
IF (x >= xmax) THEN
ErrOut('SpezFunkt2.dGamma.D9lnMC: x so big D9lnMC underflows');
RETURN 0.0;
END;
IF (x < xbig) THEN
RETURN Dcsevl(2.0*sqr(10.0/x) - 1.0,Algmcs,nalgm) / x;
ELSE
RETURN 1.0 / (12.0*x);
END;
END D9lnMC;
VAR i,n : INTEGER;
sinpiy,y : LONGREAL;
result : LONGREAL;
BEGIN (* dGamma *)
(*
VAR x:LONGREAL;
x:=1.0;
REPEAT
y := Dcsevl(x,Gamcs,ngam);
y := Dcsevl(x,Algmcs,nalgm);
TIO.WrLngReal(x,12,5); TIO.WrChar(";"); TIO.WrLngReal(y,30,20); TIO.WrLn;
x:=x - 0.01;
UNTIL (x < -1.0);
*)
IF (ABS(X) < 0.03125) AND (ABS(X) > 0.0) THEN (* MRi, 29.05.16 *)
RETURN GamSmall(X);
END;
y := ABS(X);
IF (y > 10.0) THEN
(* gamma(x) for ABS(x) .gt. 10.0. recall y = ABS(x). *)
IF (X > xmax) THEN
ErrOut('SpezFunkt2.dGamma: x so big gamma overflows');
RETURN INF;
END;
IF (X < xmin) THEN
ErrOut('SpezFunkt2.dGamma: x so small gamma underflows');
RETURN 0.0;
END;
result := exp((y-0.5)*ln(y) - y + sq2pil + D9lnMC(y));
IF (X > 0.0) THEN RETURN result;END;
IF (ABS((X - DINT(X-0.5)) / X) < dxrel ) THEN
ErrOut('SpezFunkt2.dGamma answer lt half precision');
ErrOut(' x too near negative integer');
END;
sinpiy := sin(pi*y);
IF (sinpiy = 0.0) THEN
ErrOut('SpezFunkt2.dGamma x is a negative integer');
END;
RETURN -pi/(y*sinpiy*result);
ELSE
(* compute gamma(x) for -xbnd .le. x .le. xbnd. reduce interval *)
(* and find gamma(1+y) for 0.0 .le. y .lt. 1.0 first of all. *)
n := round(X);
IF (X < 0.0) THEN DEC(n); END;
y := X - VAL(LONGREAL,n);
DEC(n);
result := 0.9375 + Dcsevl(2.0*y-1.0,Gamcs,ngam);
IF (n = 0) THEN
RETURN result;
ELSIF (n <= 0) THEN
(* compute gamma(x) for x .lt. 1.0 *)
n:=-n;
IF (X = 0.0) THEN ErrOut('SpezFunkt2.dGamma: x is 0'); END;
IF (X < 0.0) AND ((X+VAL(LONGREAL,(n-2))) = 0.0) THEN
FatalError('dgamma: x is a negative integer');
END;
IF (X < -0.5) AND (ABS((X-DINT(X-0.5)) / X) < dxrel) THEN
ErrOut('SpezFunkt2.dGamma answer lt half precision because');
ErrOut(' x too near negative integer');
END;
IF (y < xsml) THEN
ErrOut('SpezFunkt2.dGamma x is so close to 0.0 that');
ErrOut(' the result overflows');
RETURN INF;
END;
FOR i:=1 TO n DO
result :=result / (X + VAL(LONGREAL,(i-1)));
END;
RETURN result;
END; (* IF n *)
END; (* IF x < 10 *)
(* gamma(x) for x .ge. 2.0 and x .le. 10.0 *)
FOR i:=1 TO n DO
result:=result*(y+VAL(LONGREAL,i));
END;
RETURN result;
END dGamma;
PROCEDURE InitDS(VAR Koeff : ARRAY OF LONGREAL;
NKoeff : CARDINAL;
eta : LONGREAL) : CARDINAL;
(*--------------------------------------------------------------*)
(* Global für Initialisierung-Code *)
(* *)
(* Determine the number of terms needed in an orthogonal *)
(* polynomial series so that it meets a specified accuracy. *)
(* *)
(* Initialize the orthogonal series, represented by the *)
(* array Koeff, so that InitDS is the number of terms needed to *)
(* insure the error is no larger than eta. Ordinarily, eta will *)
(* be chosen to be one-tenth machine precision. *)
(* *)
(* KoeffVek : Array of coefficients in an orthogonal series. *)
(* NKoeff : Number of coefficients in OS. *)
(* eta : Scalar containing requested accuracy of series. *)
(*--------------------------------------------------------------*)
VAR err : LONGREAL;
i,ii : CARDINAL;
BEGIN
IF (NKoeff < 1) THEN
FatalError('InitDS: Number of coefficients is less than 1');
END;
err:=0.0; ii:=0;
REPEAT (* Im Original ist eta & err REAL*4 *)
INC(ii);
i := NKoeff - ii;
err:=err + ABS(Koeff[i]);
(*
TIO.WrStr("InitDS: i,ii = "); TIO.WrCard(i,3); TIO.WrCard(ii,3);
TIO.WrLn;
*)
UNTIL (err > eta) OR (ii = NKoeff);
IF (i = 0) THEN
ErrOut('InitDS: Chebyshev series too short for specified accuracy');
END;
(*
TIO.WrLn;
TIO.WrStr("InitDS: NKoeff,i = "); TIO.WrCard(NKoeff,3); TIO.WrCard(i,3);
TIO.WrLn;
*)
RETURN i+1;
END InitDS;
PROCEDURE D9GamL(VAR Xmin,Xmax : LONGREAL);
(*----------------------------------------------------------------*)
(* Calculate the minimum and maximum legal bounds for x in *)
(* dGamma(x). Xmin and Xmax are not the only bounds, but they are *)
(* the only non trivial ones to calculate. *)
(* *)
(* Output arguments *)
(* *)
(* xmin minimum legal value of x in gamma(x). Any smaller *)
(* value of x might result in underflow. *)
(* xmax maximum legal value of x in gamma(x). Any larger *)
(* value of x might cause overflow. *)
(* *)
(* Only called once in the Module initialisation code *)
(*----------------------------------------------------------------*)
VAR i : INTEGER;
alnbig,alnsml,xln,xold : LONGREAL;
found : BOOLEAN;
BEGIN
alnsml := ln(D1Mach(1));
Xmin:=-alnsml;
i:=0;
REPEAT
INC(i);
xold := Xmin; xln := ln(Xmin);
Xmin := Xmin - Xmin*((Xmin+0.5)*xln - Xmin - 0.2258 + alnsml) /
(Xmin*xln + 0.5);
found := (ABS(Xmin - xold) < 0.005);
UNTIL found OR (i > 10);
IF NOT found THEN FatalError('d9gaml: unable to find xmin'); END;
(* LOOP
* FOR i:=1 TO 10 DO
* xold := Xmin; xln := ln(Xmin);
* Xmin := Xmin - Xmin*((Xmin+0.5)*xln - Xmin - 0.2258 + alnsml) /
* (Xmin*xln + 0.5);
* IF (ABS(Xmin - xold) < 0.005) THEN EXIT; END;
* END;
* FatalError('d9gaml: unable to find xmin'); HALT;
* END; *)
Xmin := -Xmin + 0.01;
alnbig := ln(D1Mach(2));
Xmax := alnbig;
i:=0;
REPEAT
INC(i);
xold := Xmax; xln := ln(Xmax);
Xmax := Xmax - Xmax*((Xmax - 0.5)*xln - Xmax + 0.9189 - alnbig) /
(Xmax*xln - 0.5);
found := (ABS(Xmax - xold) < 0.005);
UNTIL found OR (i > 10);
IF NOT found THEN FatalError('d9gaml: unable to find xmax'); END;
(* LOOP
* FOR i:=1 TO 10 DO
* xold := Xmax; xln := ln(Xmax);
* Xmax := Xmax - Xmax*((Xmax - 0.5)*xln - Xmax + 0.9189 - alnbig) /
* (Xmax*xln - 0.5);
* IF (ABS(Xmax - xold) < 0.005) THEN EXIT; END;
* END;
* FatalError('d9gaml: unable to find xmax'); HALT;
* END; *)
Xmax := Xmax - 0.01;
Xmin := DMAX1(Xmin,-Xmax + 1.0);
END D9GamL;
BEGIN
Gamcs:=gamcs; (* For efficiency reasons *)
Algmcs:=algmcs;
(* Initialisierungcode für dGaD1mma *)
ngam := InitDS(Gamcs,42,0.1*D1Mach(3));
D9GamL(xmin,xmax);
(* xsml = 2.247436222560E-308 *)
xsml := exp(DMAX1(ln(D1Mach(1)),-ln(D1Mach(2))) + 0.01);
dxrel := sqrt(D1Mach(4));
(* Initialisierungcode für D9lnMC *)
nalgm := InitDS(Algmcs,15,0.1*D1Mach(3));
xbig := 1.0/sqrt(D1Mach(3));
D9lnMCxmax := exp(DMIN(ln(D1Mach(2)/12.0),-ln(12.0*D1Mach(1))));
END LibdGamma;
PROCEDURE Errf(x : LONGREAL) : LONGREAL;
(*----------------------------------------------------------------*)
(* Ben"otigt die PROCEDURE GammaU *)
(* Reihenentwicklung : *)
(* Errf x = {2 \over {\sqrt \pi}} \cdot *)
(* \sum_{i=0}^\infty {{-1}^i x^{2i+1} \over {(2i+1) i!}} *)
(*----------------------------------------------------------------*)
CONST genau = 1.0E-15;
Konst = 1.128379167095512570; (* 2.0/sqrt(pi) *)
VAR i : CARDINAL;
Sum,y : LONGREAL;
x2,x2i1,fact : LONGREAL;
BEGIN
IF (x < 0.0) THEN RETURN - GammaU(0.5,x*x); END;
IF (x = 0.0) THEN
RETURN 0.0;
ELSIF (x > 6.0) THEN
RETURN 1.0;
ELSIF (x < 1.0) THEN (* Reihenentwicklung *)
x2:=x*x; Sum:=x;
i:=1; fact:=1.0; x2i1:=-x*x2;
LOOP
y := x2i1 / (LFLOAT(2*i+1)*fact);
Sum:=Sum+y;
IF (ABS(y) < genau*ABS(Sum)) THEN EXIT END;
x2i1:=-x2i1*x2;
INC(i);
fact:=fact*LFLOAT(i);
END;
RETURN Konst*Sum;
ELSE (* x \in [1,6] *)
RETURN GammaU(0.5,x*x);
END;
END Errf;
MODULE ErrorFunction; (* Kapselt Daten für ErrorF,CompErrorF *)
IMPORT LowLong;
IMPORT sqrtpi,exp;
<* IF (__XDS__) THEN *>
IMPORT ErrorF,CompErrorF; (* EXPORT *)
<* ELSE *>
EXPORT ErrorF,CompErrorF;
<* END *>
TYPE Vek18 = ARRAY [1..18] OF LONGREAL;
Vek17 = ARRAY [1..17] OF LONGREAL;
CONST xup = 6.25;
A = Vek18{
1.9449071068178803E+00, 4.20186582324414E-02,
-1.86866103976769E-02, 5.1281061839107E-03,
-1.0683107461726E-03, 1.744737872522E-04,
-2.15642065714E-05, 1.7282657974E-06,
-2.00479241E-08, -1.64782105E-08,
2.0008475E-09, 2.57716E-11,
-3.06343E-11, 1.9158E-12,
3.703E-13, -5.43E-14,
-4.0E-15, 1.2E-15
};
B = Vek17{
1.4831105640848036E+00, -3.010710733865950E-01,
6.89948306898316E-02, -1.39162712647222E-02,
2.4207995224335E-03, -3.658639685849E-04,
4.86209844323E-05, -5.7492565580E-06,
6.113243578E-07, -5.89910153E-08,
5.2070091E-09, -4.232976E-010,
3.18811E-11, -2.2361E-012,
1.467E-13, -9.0E-15,
5.0E-16
};
PROCEDURE ErrorF(x : LONGREAL) : LONGREAL;
VAR t,s,s1,s2,x2 : LONGREAL;
i : CARDINAL;
BEGIN
t:=ABS(x); (* t = XV *)
IF (t <= 2.0) THEN
x2 := x*x - 2.0;
s1 := B[17];
s2 := 0.0;
s := x2*B[17] + B[16]; (* s = BJP2 *)
FOR i:=15 TO 1 BY -1 DO
s2 := s1;
s1 := s;
s := x2*s1 - s2 + B[i];
END;
RETURN (s - s2)*x / 2.0;
ELSIF (t < xup) THEN
x2 := 2.0 - 20.0 / (t + 3.0);
s1 := A[18];
s2 := 0.0;
s := x2*A[18] + A[17]; (* s = BJ *)
FOR i:=16 TO 1 BY -1 DO
s2 := s1;
s1 := s;
s := x2*s1 - s2 + A[i];
END;
x2 := ((s - s2) / (2.0*t))*exp(-x*x) / sqrtpi;
RETURN (1.0 - x2)*LowLong.sign(x);
ELSE
RETURN LowLong.sign(x);
END;
END ErrorF;
TYPE Vek23 = ARRAY [1..23] OF LONGREAL;
CONST xlow = -6.25;
xhigh = 27.28;
C = Vek23{
1.455897212750385E-01, -2.734219314954260E-01,
2.260080669166197E-01, -1.635718955239687E-01,
1.026043120322792E-01, -5.480232669380236E-02,
2.414322397093253E-02, -8.220621168415435E-03,
1.802962431316418E-03, -2.553523453642242E-05,
-1.524627476123466E-04, 4.789838226695987E-05,
3.584014089915968E-06, -6.182369348098529E-06,
7.478317101785790E-07, 6.575825478226343E-07,
-1.822565715362025E-07, -7.043998994397452E-08,
3.026547320064576E-08, 7.532536116142436E-09,
-4.066088879757269E-09, -5.718639670776992E-10,
3.328130055126039E-10
};
PROCEDURE CompErrorF(x : LONGREAL) : LONGREAL;
(*----------------------------------------------------------------*)
(* Maximale Abweichung mit F77 Routine ist *)
(* *)
(* im Intervall (0 ,2 ) *)
(* im Intervall (2 ,5 ) *)
(* im Intervall (5 ,10) *)
(* im Intervall (10,15) *)
(* im Intervall (15,25) *)
(*----------------------------------------------------------------*)
CONST xlow = -6.25;
xhigh = 27.28;
VAR t,s,polx : LONGREAL;
i : CARDINAL;
BEGIN
IF (x <= xlow) THEN
RETURN 2.0;
ELSIF (x >= xhigh) THEN
RETURN 0.0;
ELSE
t := 1.0 - 7.5 / (ABS(x) + 3.75);
polx := 0.0;
FOR i:=23 TO 1 BY -1 DO polx:=polx*t + C[i]; END;
s := exp(-x*x)*polx;
IF (x < 0.0) THEN
RETURN 2.0 - s;
ELSE
RETURN s;
END;
END;
END CompErrorF;
END ErrorFunction;
PROCEDURE CLnGamma(Z : LONGCOMPLEX) : LONGCOMPLEX; (* CERN *)
TYPE V10 = ARRAY [1..10] OF LONGREAL;
CONST C1 = 9.18938533204672742E-01;
C2 = 1.14472988584940017E+00;
Ci = V10{
8.33333333333333333E-02, -2.77777777777777778E-03,
7.93650793650793651E-04, -5.95238095238095238E-04,
8.41750841750841751E-04, -1.91752691752691753E-03,
6.41025641025641026E-03, -2.95506535947712418E-02,
1.79644372368830573E-01, -1.39243221690590112E+00
};
VAR a,x,y,t,ui,ur,ya : LONGREAL;
i : INTEGER;
u,h,p,r : LONGCOMPLEX;
BEGIN
x := RE(Z);
y := IM(Z);
IF (x <= 0.0) AND (y = 0.0) AND (ABS(x) = LFLOAT(TRUNC(ABS(x)))) THEN
h := CMPLX(INF,INF);
ErrOut("Argument equals non-positive integer (CLnGamma)");
ELSE
ya := ABS(y);
u := CMPLX(x,ya);
IF (x < 0.0) THEN
u := one - u;
END;
h := zero;
ur := RE(u);
IF (ur < 7.0) THEN
ui := IM(u);
a := arctan2(ui,ur);
h := u;
FOR i:=1 TO 6-VAL(INTEGER,ur) DO
ur := ur + 1.0;
u := CMPLX(ur,ui);
h := h*u;
a := a + arctan2(ui,ur);
END;
h := CMPLX(0.5*ln(sqr(RE(h)) + sqr(IM(h))),a);
u := u + one;
END; (* IF *)
r := one / (u*u);
p := scalarMult(Ci[10],r);
FOR i:=9 TO 2 BY - 1 DO
p := r*CMPLX(Ci[i] + RE(p), IM(p));
END;
h := CMPLX(C1,0.0) + CMPLX(RE(u)-0.5,IM(u))*LongComplexMath.ln(u) - u +
CMPLX(Ci[1] + RE(p),IM(p)) / u - h;
IF (x < 0.0) THEN
ur := LFLOAT(VAL(INTEGER,x) - 1);
ui := pi*(x - ur);
x := pi*ya;
t := exp(-x-x);
a := sin(ui);
t := x + 0.5*ln(t*(a*a) + sqr(0.5*(1.0 - t)));
a := arctan2(cos(ui)*tanh(x),a) - ur*pi;
h := CMPLX(C2 - t,a) - h;
END; (* IF *)
IF (y < 0.0) THEN
h := conj(h);
END;
END; (* IF *)
RETURN h;
END CLnGamma;
PROCEDURE CGamma(Z : LONGCOMPLEX) : LONGCOMPLEX; (* MRi *)
VAR lnZ : LONGCOMPLEX;
BEGIN
IF (Z = zero) THEN
RETURN CMPLX(INF, INF);
ELSIF (RE(Z) >= MaxGam) THEN
RETURN CMPLX(INF, INF);
ELSE
lnZ := CLnGamma(Z);
IF (lnZ = CMPLX(INF,INF)) THEN
RETURN lnZ;
ELSE
RETURN LongComplexMath.exp(lnZ);
END;
END;
END CGamma;
PROCEDURE CdGamma(z : LONGCOMPLEX) : LONGCOMPLEX;
(*----------------------------------------------------------------*)
(* This code is a Modula-2 translation of the fftlog subroutine *)
(* form github.com/coccoinomane/fftlog-f90/blob/master/cdgamma.f *)
(* published under the MIT licence. It is based on work of *)
(* Takuya Ooura, Research Institute for Mathematical Sciences, *)
(* Kyoto University with modification by Andrew Hamilton *)
(*----------------------------------------------------------------*)
PROCEDURE ANINT(x : LONGREAL) : LONGREAL;
(* Fortran ANINT function *)
BEGIN
IF (x < 0.0) THEN
RETURN LFLOAT(INT(x-0.5));
ELSE
RETURN LFLOAT(INT(x+0.5));
END;
END ANINT;
CONST IfLn = FALSE; (* calculate gamma(z) *)
Huge = LFLOAT(MAX(REAL));
OneP = 1.0 + 4.0*MachEps;
CONST PV = 7.31790632447016203E+00; PU = 3.48064577727581257E+00;
PR = 3.27673720261526849E-02;
P1 = 1.05400280458730808E+01; P2 = 4.73821439163096063E+01;
P3 = 9.11395751189899762E+01; P4 = 6.62756400966213521E+01;
P5 = 1.32280130755055088E+01; P6 = 2.93729529320536228E-01;
Q1 = 9.99999999999975753E-01; Q2 = 2.00000000000603851E+00;
Q3 = 2.99999999944915534E+00; Q4 = 4.00000003016801681E+00;
Q5 = 4.99999857982434025E+00; Q6 = 6.00009857740312429E+00;
VAR xr,xi,yr,yi,ur,ui,vr,vi,wr,wi,t : LONGREAL;
BEGIN
xr := RE(z); xi := IM(z);
IF (ABS(xr) < 2.0*MachEps) AND (ABS(xi) < 2.0*MachEps) THEN
ErrOut("CdGamma) is (+infinity,+infinity) (CdGamma)");
yi := INF; yr := INF;
ELSIF (xr <= 0.0) AND (ABS(xi) <= 2.0*MachEps) AND IsNearInteger(xr) THEN
(* x = 0, -1, -2 *)
IF IfLn THEN (* lngamma *)
(* real part is +infinity *)
ErrOut("RE(CdGamma) is + infinity (CdGamma)");
yr := INF;
yi := pi*LFLOAT(INT(xr));
ELSE (* gamma *)
IF ODD(INT(xr*OneP)) THEN
yr := -INF;
ErrOut("RE(CdGamma) is - infinity (CdGamma)");
ELSE (* + infinity at even negative integers *)
yr := INF;
ErrOut("RE(CdGamma) is + infinity (CdGamma)");
END;
yi := 0.0;
END; (* IfLn *)
ELSE
(* re(x) < 1/2 : use reflection formula *)
IF (xr < 0.5) THEN
wr := 1.0 - xr;
wi := -xi;
ELSE
wr := xr;
wi := xi;
END;
(* large |x|: keep only leading term of rational function *)
t := wr*wr + wi*wi;
IF (t > Huge) THEN (* Assume we calculete with LONGREAL *)
(* rational function v *)
vr := wr/t + PR;
vi := -wi/t;
(* ln(overall factor) u = ln(x + pv) - 1 *)
yr := wr + PV;
ur := yr; IF (ur < 0.0) THEN ur := -ur; END;
ui := wi; IF (ui < 0.0) THEN ui := -ui; END;
IF (ur >= ui) THEN
t := wi / yr;
ur := ln(ur) + 0.5*ln(1.0 + t*t) - 1.0;
ELSE
t := yr/wi;
ur := ln(ui) + 0.5*ln(1.0 + t*t) - 1.0;
END;
ui := arctan2(wi,yr);
ELSE
(* not large |x| *)
ur := wr + Q6;
vr := ur*(wr + Q5) - wi*wi;
vi := wi*(wr + Q5) + ur*wi;
yr := P6 + (P5*ur + P4*vr);
yi := P5*wi + P4*vi;
ur := vr*(wr + Q4) - vi*wi;
ui := vi*(wr + Q4) + vr*wi;
vr := ur*(wr + Q3) - ui*wi;
vi := ui*(wr + Q3) + ur*wi;
yr := yr + (P3*ur + P2*vr);
yi := yi + (P3*ui + P2*vi);
ur := vr*(wr + Q2) - vi*wi;
ui := vi*(wr + Q2) + vr*wi;
vr := ur*(wr + Q1) - ui*wi;
vi := ui*(wr + Q1) + ur*wi;
yr := yr + (P1*ur + vr);
yi := yi + (P1*ui + vi);
ur := vr*wr - vi*wi;
ui := vi*wr + vr*wi;
t := ur*ur + ui*ui;
vr := (yr*ur + yi*ui)/t + PR;
vi := (yi*ur - yr*ui)/t;
yr := wr + PV;
ur := 0.5*ln(yr*yr + wi*wi) - 1.0;
ui := arctan2(wi,yr);
END; (* IF t *)
yr := ur*(wr - 0.5) - ui*wi - PU;
yi := ui*(wr - 0.5) + ur*wi;
yr := yr + 0.5*ln(vr*vr + vi*vi);
yi := yi + arctan2(vi,vr);
t := 1.0;
IF (xr < 0.5) THEN
wi := ANINT(xr);
wr := xr - wi;
IF (wi > xr) THEN wi := wi - 1.0; END;
IF (xi = 0.0) THEN (* case of real x *)
wr := ln(sin(pi*ABS(wr)));
IF IfLn THEN
wi := -pi*wi;
ELSE
IF (wi # 2.0*LFLOAT(INT(0.5*wi))) THEN t := -1.0; END;
wi:=0.0;
END;
ELSIF (ABS(xi) < 1.0) THEN
(* imaginary part of x is < 1 in absolute value *)
IF IfLn THEN
ui := -pi*wi;
IF (xi < 0.0) THEN ui := -ui; END;
ELSE
IF (wi # 2.0*LFLOAT(INT(0.5*wi))) THEN
t := -1.0;
END;
ui := 0.0;
END; (* IfLn *)
wr := pi*wr;
wi := pi*xi;
vr := sin(wr)*cosh(wi);
vi := cos(wr)*sinh(wi);
IF (wr < 0.0) THEN
vr := -vr;
vi := -vi;
END;
wr := 0.5*ln(vr*vr + vi*vi);
wi := ui + arctan2(vi,vr);
ELSE
(* imaginary part of x is >= 1 in absolute value *)
IF IfLn THEN
ui := pi*(0.5 - xr);
ELSE
IF (wi # 2.0*LFLOAT(INT(0.5*wi))) THEN
t := -1.0;
END; (* IF *)
IF (wr >= 0.0) THEN
ui := pi*(+0.5 - wr);
ELSE
ui := pi*(-0.5 - wr);
END;
END;
wi := exp(-2.0*pi*ABS(xi));
wr := 2.0*pi*wr;
vr := 0.5*(1.0 - cos(wr)*wi);
vi := -0.5*sin(wr)*wi;
ur := pi*xi;
(* w = ln[sin(pi x)] *)
IF (xi > 0.0) THEN
wr := ur + 0.5*ln(vr*vr + vi*vi);
wi := ui + arctan2(vi,vr);
ELSIF (xi < 0.0) THEN
wr := -ur + 0.5*ln(vr*vr + vi*vi);
wi := -ui - arctan2(vi,vr);
END;
END; (* IF xi *)
(* y = ln[gamma(x)] *)
yr := ln(pi) - yr - wr;
yi := -yi - wi;
END; (* IF *)
IF NOT IfLn THEN (* gamma *)
ur := exp(yr);
IF (xi = 0.0) THEN
yr := t*ur;
yi := 0.0;
ELSE
yr := t*ur*cos(yi);
yi := ur*sin(yi);
END;
END;
END; (* IF NOT z = 0, -1, -2, ... *)
RETURN CMPLX(yr,yi);
END CdGamma;
PROCEDURE IBeta(A,B : LONGREAL;
X : LONGREAL) : LONGREAL;
PROCEDURE PSeries(A,B : LONGREAL;
X : LONGREAL) : LONGREAL;
(*---------------------------------------------------------------*)
(* Power series for incomplete beta integral. *)
(* Use when B*X is small *)
(*---------------------------------------------------------------*)
VAR S,T,U,V,T1,Z,Ai : LONGREAL;
N : INTEGER;
BEGIN
Ai := 1.0 / A;
U := (1.0 - B) * X;
V := U / (A + 1.0);
T1 := V;
T := U;
N := 2;
S := 0.0;
Z := MachEps*Ai;
WHILE (ABS(V) > Z) DO
U := (LFLOAT(N) - B) * X / LFLOAT(N);
T := T * U;
V := T / (A + LFLOAT(N));
S := S + V;
INC(N);
END;
S := S + T1;
S := S + Ai;
U := A*ln(X);
IF (A + B < MaxGam) AND (ABS(U) < MaxLog) THEN
T := dGamma(A + B) / (dGamma(A)*dGamma(B));
S:=S * T*power(X,A);
ELSE
T := dLnGamma(A + B) - dLnGamma(A) - dLnGamma(B) + U + ln(S);
IF T < MinLog THEN
S := 0.0
ELSE
S := exp(T);
END;
END;
RETURN S;
END PSeries;
PROCEDURE CFrac1(A,B : LONGREAL;
X : LONGREAL) : LONGREAL;
(* Continued fraction expansion #1 for incomplete beta integral *)
CONST Big = 1.0 / MachEps;
Thresh = 3.0 * MachEps;
VAR Xk,Pk,Pkm1,Pkm2,Qk,Qkm1,Qkm2 : LONGREAL;
K1,K2,K3,K4,K5,K6,K7,K8 : LONGREAL;
R,T,Ans : LONGREAL;
N : INTEGER;
BEGIN
K1 := A;
K2 := A + B;
K3 := A;
K4 := A + 1.0;
K5 := 1.0;
K6 := B - 1.0;
K7 := K4;
K8 := A + 2.0;
Pkm2 := 0.0;
Qkm2 := 1.0;
Pkm1 := 1.0;
Qkm1 := 1.0;
Ans := 1.0;
R := 1.0;
N := 0;
LOOP
Xk := - (X * K1 * K2) / (K3 * K4);
Pk := Pkm1 + Pkm2 * Xk;
Qk := Qkm1 + Qkm2 * Xk;
Pkm2 := Pkm1;
Pkm1 := Pk;
Qkm2 := Qkm1;
Qkm1 := Qk;
Xk := (X * K5 * K6) / (K7 * K8);
Pk := Pkm1 + Pkm2 * Xk;
Qk := Qkm1 + Qkm2 * Xk;
Pkm2 := Pkm1;
Pkm1 := Pk;
Qkm2 := Qkm1;
Qkm1 := Qk;
IF (Qk <> 0.0) THEN R := Pk / Qk; END;
IF (R <> 0.0) THEN
T := ABS((Ans - R) / R);
Ans := R;
ELSE T := 1.0;END;
IF (T < Thresh) THEN EXIT; END;
K1 := K1 + 1.0;
K2 := K2 + 1.0;
K3 := K3 + 2.0;
K4 := K4 + 2.0;
K5 := K5 + 1.0;
K6 := K6 - 1.0;
K7 := K7 + 2.0;
K8 := K8 + 2.0;
IF (ABS(Qk) + ABS(Pk) > Big) THEN
Pkm2:=Pkm2 * MachEps;
Pkm1:=Pkm1 * MachEps;
Qkm2:=Qkm2 * MachEps;
Qkm1:=Qkm1 * MachEps;
END;
IF (ABS(Qk) < MachEps) OR (ABS(Pk) < MachEps) THEN
Pkm2:=Pkm2 * Big;
Pkm1:=Pkm1 * Big;
Qkm2:=Qkm2 * Big;
Qkm1:=Qkm1 * Big;
END;
INC(N);
IF (N > 400) THEN
Fehler:=TRUE;
Fehlerflag:="FPLoss (CFrac1)";
END;
END; (* LOOP *)
RETURN Ans;
END CFrac1;
PROCEDURE CFrac2(A,B : LONGREAL;
X : LONGREAL) : LONGREAL;
(* Continued fraction expansion #2 for incomplete beta integral *)
CONST Big = 1.0 / MachEps;
Thresh = 3.0 * MachEps;
VAR Xk,Pk,Pkm1,Pkm2,Qk,Qkm1,Qkm2 : LONGREAL;
K1,K2,K3,K4,K5,K6,K7,K8 : LONGREAL;
R,T,Z,Ans : LONGREAL;
N : INTEGER;
BEGIN
K1 := A;
K2 := B - 1.0;
K3 := A;
K4 := A + 1.0;
K5 := 1.0;
K6 := A + B;
K7 := A + 1.0;
K8 := A + 2.0;
Pkm2 := 0.0;
Qkm2 := 1.0;
Pkm1 := 1.0;
Qkm1 := 1.0;
Z := X / (1.0 - X);
Ans := 1.0;
R := 1.0;
N := 0;
LOOP
Xk := - (Z * K1 * K2) / (K3 * K4);
Pk := Pkm1 + Pkm2 * Xk;
Qk := Qkm1 + Qkm2 * Xk;
Pkm2 := Pkm1;
Pkm1 := Pk;
Qkm2 := Qkm1;
Qkm1 := Qk;
Xk := (Z * K5 * K6) / (K7 * K8);
Pk := Pkm1 + Pkm2 * Xk;
Qk := Qkm1 + Qkm2 * Xk;
Pkm2 := Pkm1;
Pkm1 := Pk;
Qkm2 := Qkm1;
Qkm1 := Qk;
IF (Qk <> 0.0) THEN R := Pk / Qk; END;
IF (R <> 0.0) THEN
T := ABS((Ans - R) / R);
Ans := R;
ELSE
T := 1.0;
END;
IF (T < Thresh) THEN EXIT; END;
K1 := K1 + 1.0;
K2 := K2 - 1.0;
K3 := K3 + 2.0;
K4 := K4 + 2.0;
K5 := K5 + 1.0;
K6 := K6 + 1.0;
K7 := K7 + 2.0;
K8 := K8 + 2.0;
IF (ABS(Qk) + ABS(Pk) > Big) THEN
Pkm2:=Pkm2 * MachEps;
Pkm1:=Pkm1 * MachEps;
Qkm2:=Qkm2 * MachEps;
Qkm1:=Qkm1 * MachEps;
END;
IF (ABS(Qk) < MachEps) OR (ABS(Pk) < MachEps) THEN
Pkm2:=Pkm2 * Big;
Pkm1:=Pkm1 * Big;
Qkm2:=Qkm2 * Big;
Qkm1:=Qkm1 * Big;
END;
INC(N);
N := N + 1;
IF (N > 400) THEN
Fehler:=TRUE;
Fehlerflag:="FPLoss (CFrac1)";
END;
END; (* LOOP *)
RETURN Ans;
END CFrac2;
VAR A1,B1,X1,T,W,Xc,Y : LONGREAL;
flag : BOOLEAN;
BEGIN
IF (A <= 0.0) OR (B <= 0.0) OR (X < 0.0) THEN
RETURN NAN;
END;
IF (X > 1.0) THEN
RETURN 1.0;
END;
IF (X = 0.0) OR (X = 1.0) THEN
RETURN X;
END;
flag := FALSE;
IF (B * X <= 1.0) AND (X <= 0.95) THEN
T := PSeries(A, B, X);
ELSE
W := 1.0 - X;
IF (X > (A / (A + B))) THEN
(* Reverse a and b if x is greater than the mean. *)
flag := TRUE;
A1 := B;
B1 := A;
Xc := X;
X1 := W;
ELSE
A1 := A;
B1 := B;
Xc := W;
X1 := X;
END;
IF flag AND (B1 * X1 <= 1.0) AND (X1 <= 0.95) THEN
T := PSeries(A1, B1, X1);
ELSE
(* Choose expansion for optimal convergence *)
Y := X1 * (A1 + B1 - 2.0) - (A1 - 1.0);
IF Y < 0.0 THEN
W := CFrac1(A1, B1, X1)
ELSE
W := CFrac2(A1, B1, X1) / Xc;
END;
(* Multiply w by the factor *)
(* a b _ _ _ *)
(* x (1-x) | (a+b) / ( a | (a) | (b) ) *)
Y := A1*ln(X1);
T := B1*ln(Xc);
IF (A1 + B1 < MaxGam) AND (ABS(Y) < MaxLog) AND (ABS(T) < MaxLog) THEN
T:= power(Xc,B1);
T:=T * power(X1,A1);
T:=T / A1;
T:=T * W;
T:=T * dGamma(A1 + B1) / (dGamma(A1)*dGamma(B1));
ELSE
(* Resort to logarithms *)
Y:=Y + T + dLnGamma(A1 + B1) - dLnGamma(A1)
- dLnGamma(B1) + ln(W/A1);
IF Y < MinLog THEN
T := 0.0
ELSE
T := exp(Y);
END;
END;
END; (* IF ELSE *)
END; (* IF ELSE *)
IF flag THEN
IF (T <= MachEps) THEN
T := 1.0 - MachEps
ELSE
T := 1.0 - T;
END;
END;
RETURN T;
END IBeta;
PROCEDURE SetGamitFast(opt : BOOLEAN);
BEGIN
fast:=opt;
END SetGamitFast;
PROCEDURE dGamit(A,X : LONGREAL): LONGREAL; (* SLATEC *)
PROCEDURE dLnGamSmall( X : LONGREAL;
VAR dLnGam : LONGREAL;
VAR SignGam : LONGREAL);
(*--------------------------------------------------------------*)
(* dLnGamSmall(x,dLnGam,SignGam) calculates the natural *)
(* logarithm of the absolute value of the gamma function for *)
(* real argumentx and stores the result in argument dLnGam. *)
(*--------------------------------------------------------------*)
VAR iFehl : INTEGER;
BEGIN (* +++ *)
IF fast THEN
dLnGam := aLnGamma(ABS(X),iFehl);
ELSE
dLnGam := dLnGamma(ABS(X));
END;
SignGam := 1.0;
IF (X <= 0.0) AND (TRUNC((X / 2.0) + 0.1) = 0) THEN
SignGam := - 1.0;
END;
END dLnGamSmall;
PROCEDURE d9Gamit( A,X : LONGREAL; (* SLATEC *)
VAR aLnAp1,SignGam : LONGREAL) : LONGREAL;
(*--------------------------------------------------------------*)
(* Compute Tricomis incomplete gamma function for small x. *)
(*--------------------------------------------------------------*)
CONST eps = 0.5*1.11E-16;
VAR RetVal,s,sgng2,t,te : LONGREAL;
ae,aeps,algs,alg2,fk : LONGREAL;
k,m,ma : INTEGER;
BEGIN
IF (X <= 0.0) THEN
ErrOut('x should be > 0 (d9gmit)');
(* Fehler:=TRUE; *)
RETURN MAX(LONGREAL);
END;
ma := VAL(INTEGER,(A + 0.5));
IF (A < 0.0) THEN
ma := VAL(INTEGER,(A - 0.5));
END;
aeps := A - VAL(LONGREAL,ma);
ae := A;
IF (A < -0.5) THEN
ae := aeps;
END;
te := ae;
s := 1.0;
k:=0;
REPEAT
k:=k + 1;
fk := VAL(LONGREAL,k);
te := -X*te / fk;
t := te / (ae + fk);
s := s + t;
UNTIL (ABS(t) < eps*ABS(s)) OR (k > 200);
IF (k > 200) THEN
ErrOut('No convergence in 200 terms of taylor-s series');
END;
IF (A >= -0.5) THEN
algs := -aLnAp1 + ln(s);
RetVal := exp(algs);
ELSE
algs := - dLnGamma(1.0 + aeps) + ln(s);
s := 1.0;
m := -ma - 1;
IF (m # 0) THEN
t := 1.0;
k:=0;
REPEAT
k:=k + 1;
t := X*t / (aeps - VAL(LONGREAL,(m+1-k)));
s := s + t;
UNTIL (ABS(t) < eps*ABS(s)) OR (k = m);
END;
RetVal := 0.0;
algs := -VAL(LONGREAL,ma)*ln(X) + algs;
IF (s = 0.0) OR (aeps = 0.0) THEN
RetVal := exp(algs);
ELSE
sgng2 := SignGam*DSIGN(1.0,s);
alg2 := -X - aLnAp1 + ln(ABS(s));
IF (alg2 > bot) THEN
RetVal := sgng2*exp(alg2);
END;
IF (algs > bot) THEN
RetVal:=RetVal + exp(algs);
END;
END;
END;
RETURN RetVal;
END d9Gamit;
PROCEDURE dGamR(X : LONGREAL): LONGREAL; (* SLATEC *)
(*--------------------------------------------------------------*)
(* dGamR(x) calculates the reciprocal of the complete gamma *)
(* function for real argument x. *)
(*--------------------------------------------------------------*)
VAR alngx,sgngx : LONGREAL;
BEGIN
IF (X <= 0.0) AND (VAL(LONGREAL,(VAL(INTEGER,X))) = X ) THEN
RETURN 0.0;
ELSE
IF (ABS(X) > 10.0) THEN
dLnGamSmall(X,alngx,sgngx);
RETURN sgngx*exp(-alngx);
ELSE
RETURN 1.0 / dGamma(X);
END;
END;
END dGamR;
VAR aeps,ainta,aLnAp1 : LONGREAL;
sga,SignGam,RetVal : LONGREAL;
iFehl : INTEGER;
BEGIN
sga := 1.0;
IF (A # 0.0) THEN
sga := DSIGN(1.0,A);
END;
ainta := VAL(LONGREAL,VAL(INTEGER,A + 0.5*sga));
aeps := A - ainta;
IF (X <= 0.0) THEN
RetVal := 0.0;
IF (ainta > 0.0) OR (aeps # 0.0) THEN
RetVal := dGamR(A+1.0);
END;
ELSIF (X <= 1.0) THEN (* Special case for small X *)
IF (A >= -0.5) OR (aeps # 0.0) THEN
dLnGamSmall(A+1.0,aLnAp1,SignGam);
END;
RetVal := d9Gamit(A,X,aLnAp1,SignGam);
ELSE
IF (A = 0.0) THEN
(* Bitte pruefen, das kann nicht funktionieren !!! *)
RetVal := GammaU(0.0,X) / sqrt(X);
ELSE
IF fast THEN
RetVal := GammaIn(X,A,iFehl) / power(X,A);
ELSE
RetVal := GammaU(A ,X) / power(X,A);
END;
END;
END;
RETURN RetVal;
END dGamit;
BEGIN
Sqrt2Pi := sqrt(2.0*pi);
bot := ln(1.0/MAX(LONGREAL));
fast := TRUE;
END SpezFunkt2.