|
a/SpezFunkt2.mod.m2pp |
|
b/SpezFunkt2.mod.m2pp |
|
... |
|
... |
| 41 |
(* 27.08.16, MRi: Argument von dGamma ist nicht mehr "VAR" *)
|
41 |
(* 27.08.16, MRi: Argument von dGamma ist nicht mehr "VAR" *)
|
| 42 |
(* 27.08.16, MRi: Erstellen der ersten Verionsn von IGamma,JGamma,IBeta *)
|
42 |
(* 27.08.16, MRi: Erstellen der ersten Verionsn von IGamma,JGamma,IBeta *)
|
| 43 |
(* 16.02.17, MRi: Erstellen der ersten Version (aLnGamma,GamamaIn,dGamit) *)
|
43 |
(* 16.02.17, MRi: Erstellen der ersten Version (aLnGamma,GamamaIn,dGamit) *)
|
| 44 |
(* 14.02.18, MRi: Hinzufuegen von CLnGamma und CGamma *)
|
44 |
(* 14.02.18, MRi: Hinzufuegen von CLnGamma und CGamma *)
|
| 45 |
(* 14.02.18, MRi: Erstellen der ersten Version von CdGamma *)
|
45 |
(* 14.02.18, MRi: Erstellen der ersten Version von CdGamma *)
|
| 46 |
(* 23.05.18, MRi: Erweiterung von CdGamma und aufnahme in SpezFunkt2 *)
|
46 |
(* 23.05.18, MRi: Erweiterung von CdGamma und Aufnahme in SpezFunkt2 *)
|
|
|
47 |
(* 19.12.18, MRi: Ersetzen von GammaIn durch Neuuebersetzng von AS239 *)
|
|
|
48 |
(* da die alte Routine (basierend auf AS32) fuer grosse X *)
|
|
|
49 |
(* teilweise falsche Ergebnisse erzielt hatte. *)
|
| 47 |
(*------------------------------------------------------------------------*)
|
50 |
(*------------------------------------------------------------------------*)
|
| 48 |
(* Testroutinen vorhanden, IBeta gegen unabhaengigen Fortran-Code ge- *)
|
51 |
(* Testroutinen vorhanden, IBeta gegen unabhaengigen Fortran-Code ge- *)
|
| 49 |
(* testet, IGamma & JGamma nur gegen Pascal-Aequivalent *)
|
52 |
(* testet, IGamma & JGamma nur gegen Pascal-Aequivalent *)
|
| 50 |
(* Testcode (parially against Fortran versions) in TstSpzFun2 *)
|
53 |
(* Testcode (parially against Fortran versions) in TstSpzFun2 *)
|
| 51 |
(* dGamit wird in BoysInt2Lib intensiv genutzt *)
|
54 |
(* dGamit wird in BoysInt2Lib intensiv genutzt *)
|
| 52 |
(*------------------------------------------------------------------------*)
|
55 |
(*------------------------------------------------------------------------*)
|
| 53 |
(* Offene Punkte *)
|
56 |
(* Offene Punkte *)
|
| 54 |
(* *)
|
57 |
(* *)
|
| 55 |
(* - Modulstruktur durchsehen und eventuell bereinigen. *)
|
58 |
(* - Modulstruktur durchsehen und eventuell bereinigen. *)
|
| 56 |
(* - besserer Test fuer CLnGamma, CGamma *)
|
59 |
(* - besserer Test fuer CLnGamma, CGamma *)
|
| 57 |
(* - Zusammenfuehren von SpezFunkt2 und SpezFunkt2b *)
|
|
|
| 58 |
(*------------------------------------------------------------------------*)
|
60 |
(*------------------------------------------------------------------------*)
|
| 59 |
(* Implementation : Michael Riedl *)
|
61 |
(* Implementation : Michael Riedl *)
|
| 60 |
(* Licence : GNU Lesser General Public License (LGPL) *)
|
62 |
(* Licence : GNU Lesser General Public License (LGPL) *)
|
| 61 |
(*------------------------------------------------------------------------*)
|
63 |
(*------------------------------------------------------------------------*)
|
| 62 |
|
64 |
|
|
... |
|
... |
| 69 |
FROM LongMath IMPORT pi,sqrt,power,ln,exp,sin,cos,round;
|
71 |
FROM LongMath IMPORT pi,sqrt,power,ln,exp,sin,cos,round;
|
| 70 |
FROM LMathLib IMPORT INF,NAN,MachEps,MinLog,MaxLog,sqrtpi,
|
72 |
FROM LMathLib IMPORT INF,NAN,MachEps,MinLog,MaxLog,sqrtpi,
|
| 71 |
sinh,cosh,tanh,arctan2;
|
73 |
sinh,cosh,tanh,arctan2;
|
| 72 |
IMPORT LongComplexMath;
|
74 |
IMPORT LongComplexMath;
|
| 73 |
FROM LongComplexMath IMPORT zero,one,conj,scalarMult;
|
75 |
FROM LongComplexMath IMPORT zero,one,conj,scalarMult;
|
| 74 |
FROM F77func IMPORT D1Mach;
|
76 |
FROM F77func IMPORT D1Mach,DMIN;
|
| 75 |
|
77 |
|
| 76 |
VAR fast : BOOLEAN;
|
78 |
VAR fast : BOOLEAN;
|
| 77 |
Sqrt2Pi : LONGREAL; (* sqrt(2 pi) *)
|
79 |
Sqrt2Pi : LONGREAL; (* sqrt(2 pi) *)
|
| 78 |
bot : LONGREAL; (* wird in d9Gamit genutzt *)
|
80 |
bot : LONGREAL; (* wird in d9Gamit genutzt *)
|
| 79 |
|
81 |
|
|
... |
|
... |
| 87 |
ELSE
|
89 |
ELSE
|
| 88 |
RETURN - ABS(a);
|
90 |
RETURN - ABS(a);
|
| 89 |
END;
|
91 |
END;
|
| 90 |
END DSIGN;
|
92 |
END DSIGN;
|
| 91 |
|
93 |
|
|
|
94 |
PROCEDURE AlNorm(X : LONGREAL;
|
|
|
95 |
Upper : BOOLEAN) : LONGREAL;
|
|
|
96 |
|
|
|
97 |
(*----------------------------------------------------------------*)
|
|
|
98 |
(* Algorithm AS66 Applied Statistics (1973) vol22 no.3 *)
|
|
|
99 |
(* *)
|
|
|
100 |
(* Evaluates the tail area of the standardised normal curve *)
|
|
|
101 |
(* from x to infinity if upper is .true. or *)
|
|
|
102 |
(* from minus infinity to x if upper is .false. *)
|
|
|
103 |
(*----------------------------------------------------------------*)
|
|
|
104 |
|
|
|
105 |
(* machine dependent constants *)
|
|
|
106 |
|
|
|
107 |
CONST ltone = 7.0; utzero = 18.66; con = 1.28;
|
|
|
108 |
|
|
|
109 |
p = 0.398942280444; q = 0.39990348504;
|
|
|
110 |
r = 0.398942280385;
|
|
|
111 |
a1 = 5.75885480458; a2 = 2.62433121679;
|
|
|
112 |
a3 = 5.92885724438;
|
|
|
113 |
b1 = -29.8213557807; b2 = 48.6959930692;
|
|
|
114 |
c1 = -3.8052E-08; c2 = 3.98064794E-04;
|
|
|
115 |
c3 = -0.151679116635; c4 = 4.8385912808;
|
|
|
116 |
c5 = 0.742380924027; c6 = 3.99019417011;
|
|
|
117 |
d1 = 1.00000615302; d2 = 1.98615381364;
|
|
|
118 |
d3 = 5.29330324926; d4 = -15.1508972451;
|
|
|
119 |
d5 = 30.789933034;
|
|
|
120 |
|
|
|
121 |
VAR z,y,alnorm : LONGREAL;
|
|
|
122 |
up : BOOLEAN;
|
|
|
123 |
BEGIN
|
|
|
124 |
up := Upper;
|
|
|
125 |
z := X;
|
|
|
126 |
IF (z < 0.0) THEN
|
|
|
127 |
up := NOT up;
|
|
|
128 |
z := -z;
|
|
|
129 |
END; (* IF *)
|
|
|
130 |
IF (((z <= ltone) OR up) AND (z <= utzero)) THEN
|
|
|
131 |
y := 0.5*z*z;
|
|
|
132 |
IF (z > con) THEN
|
|
|
133 |
alnorm := r*exp(-y) / (z + c1 + d1 / (z + c2 + d2 /
|
|
|
134 |
(z + c3 + d3 / (z + c4 + d4 /
|
|
|
135 |
(z + c5 + d5 / (z + c6))))));
|
|
|
136 |
ELSE
|
|
|
137 |
alnorm := 0.5 - z*(p-q*y / (y + a1 + b1 /
|
|
|
138 |
(y + a2 + b2 / (y + a3))));
|
|
|
139 |
END;
|
|
|
140 |
ELSE
|
|
|
141 |
alnorm := 0.0;
|
|
|
142 |
END; (* IF *)
|
|
|
143 |
IF (NOT up) THEN
|
|
|
144 |
alnorm := 1.0 - alnorm;
|
|
|
145 |
END;
|
|
|
146 |
RETURN alnorm;
|
|
|
147 |
END AlNorm;
|
|
|
148 |
|
| 92 |
PROCEDURE GammaIn( X : LONGREAL; (* SLATEC *)
|
149 |
PROCEDURE GammaIn( X : LONGREAL;
|
| 93 |
P : LONGREAL;
|
150 |
P : LONGREAL;
|
| 94 |
VAR iFehl : INTEGER) : LONGREAL;
|
151 |
VAR Ifault : INTEGER) : LONGREAL;
|
| 95 |
|
152 |
|
|
|
153 |
(*----------------------------------------------------------------*)
|
|
|
154 |
(* Auxiliary functions required: ALOGAM := logarithm of the gamma *)
|
|
|
155 |
(* function, and AlNorm := algorithm AS66 *)
|
|
|
156 |
(*----------------------------------------------------------------*)
|
|
|
157 |
|
| 96 |
CONST eps = 1.0E-08;
|
158 |
CONST Tol = 1.0E-14;
|
| 97 |
Big = 1.0E+37;
|
159 |
Xbig = 1.0E+08; Oflow = MAX(REAL);
|
| 98 |
Small = 1.0E-37;
|
160 |
Elimit = -88.0; Plimit = 1000.0;
|
| 99 |
|
161 |
|
| 100 |
VAR a,an,arg,b,dif,factor : LONGREAL;
|
162 |
VAR pn1,pn2,pn3,pn4,pn5,pn6 : LONGREAL;
|
| 101 |
g,gin,rn,term,GInc : LONGREAL;
|
163 |
a,b,c,an,rn,arg,gamma : LONGREAL;
|
| 102 |
i : INTEGER;
|
|
|
| 103 |
pn : ARRAY [1..6] OF LONGREAL;
|
|
|
| 104 |
BEGIN
|
164 |
BEGIN
|
| 105 |
(* check the input. *)
|
165 |
(* Check that we have valid values for X and P *)
|
| 106 |
IF (P <= 0.0) THEN
|
166 |
IF ((P <= 0.0) OR (X < 0.0)) THEN
|
| 107 |
iFehl := 1;
|
167 |
Ifault := 1;
|
| 108 |
RETURN 0.0;
|
168 |
RETURN 0.0;
|
|
|
169 |
END; (* IF *)
|
|
|
170 |
Ifault := 0;
|
| 109 |
ELSIF (X < 0.0) THEN
|
171 |
IF (X = 0.0) THEN
|
| 110 |
iFehl:=2;
|
|
|
| 111 |
RETURN 0.0;
|
172 |
RETURN 0.0;
|
| 112 |
ELSIF (X = 0.0) THEN
|
173 |
END; (* IF *)
|
| 113 |
iFehl:=0;
|
174 |
(* Use a normal approximation if P > Plimit *)
|
|
|
175 |
IF (P > Plimit) THEN
|
|
|
176 |
pn1 := 3.0*sqrt(P)*((X/P)**(1.0/3.0)+1.0/(9.0*P)-1.0);
|
|
|
177 |
RETURN AlNorm(pn1,FALSE);
|
|
|
178 |
END; (* IF *)
|
|
|
179 |
(* If X is extremely large compared to P then set gamma = 1 *)
|
|
|
180 |
IF (X > Xbig) THEN
|
| 114 |
RETURN 0.0;
|
181 |
RETURN 1.0;
|
| 115 |
END;
|
182 |
END; (* IF *)
|
| 116 |
|
183 |
IF (X <= 1.0) OR (X < P) THEN
|
| 117 |
g := aLnGamma(P,iFehl);
|
184 |
(* Use Pearson's series expansion. *)
|
| 118 |
|
185 |
(* Note that P is not large enough to force overflow in ALOGAM. *)
|
| 119 |
IF (iFehl # 0) THEN
|
186 |
(* No need to test IFAULT on exit since P > 0. *)
|
| 120 |
iFehl:=4;
|
187 |
arg := P*ln(X) - X - aLnGamma(P+1.0,Ifault);
|
| 121 |
RETURN 0.0;
|
|
|
| 122 |
END;
|
|
|
| 123 |
|
|
|
| 124 |
arg := P*ln(X) - X - g;
|
|
|
| 125 |
|
|
|
| 126 |
IF (arg < ln(Small)) THEN
|
|
|
| 127 |
iFehl:=3;
|
|
|
| 128 |
RETURN 0.0;
|
|
|
| 129 |
END;
|
|
|
| 130 |
|
|
|
| 131 |
iFehl:=0;
|
|
|
| 132 |
factor := exp(arg);
|
|
|
| 133 |
|
|
|
| 134 |
IF (X <= 1.0) OR (X < P) THEN (* calculation by series expansion. *)
|
|
|
| 135 |
|
|
|
| 136 |
gin := 1.0;
|
188 |
c := 1.0;
|
| 137 |
term := 1.0;
|
189 |
gamma := 1.0;
|
| 138 |
rn := P;
|
190 |
a := P;
|
| 139 |
LOOP
|
191 |
LOOP
|
| 140 |
rn := rn + 1.0;
|
192 |
a := a + 1.0;
|
| 141 |
term := term*X / rn;
|
193 |
c := c*X/a;
|
| 142 |
gin := gin + term;
|
194 |
gamma := gamma + c;
|
| 143 |
IF (term <= eps) THEN
|
195 |
IF (c <= Tol) THEN
|
|
|
196 |
arg := arg + ln(gamma);
|
|
|
197 |
gamma := 0.0;
|
|
|
198 |
IF (arg >= Elimit) THEN
|
|
|
199 |
gamma := exp(arg);
|
|
|
200 |
END; (* IF *)
|
| 144 |
EXIT;
|
201 |
EXIT;
|
| 145 |
END;
|
202 |
END; (* IF *)
|
| 146 |
END; (* LOOP *)
|
203 |
END; (* LOOP *)
|
| 147 |
GInc := gin*factor / P;
|
204 |
ELSE
|
| 148 |
|
205 |
(* Use a continued fraction expansion *)
|
| 149 |
ELSE (* calculation by continued fraction. *)
|
206 |
arg := P*ln(X) - X - aLnGamma(P,Ifault);
|
| 150 |
|
|
|
| 151 |
a := 1.0 - P;
|
207 |
a := 1.0 - P;
|
| 152 |
b := a + X + 1.0;
|
208 |
b := a + X + 1.0;
|
| 153 |
term := 0.0;
|
209 |
c := 0.0;
|
| 154 |
|
|
|
| 155 |
pn[1] := 1.0;
|
210 |
pn1 := 1.0;
|
| 156 |
pn[2] := X;
|
211 |
pn2 := X;
|
| 157 |
pn[3] := X + 1.0;
|
212 |
pn3 := X + 1.0;
|
| 158 |
pn[4] := X*b;
|
213 |
pn4 := X*b;
|
| 159 |
|
|
|
| 160 |
gin := pn[3] / pn[4];
|
214 |
gamma := pn3/pn4;
|
| 161 |
|
|
|
| 162 |
LOOP
|
215 |
LOOP
|
| 163 |
a:=a + 1.0;
|
216 |
a := a + 1.0;
|
| 164 |
b:=b + 2.0;
|
217 |
b := b + 2.0;
|
| 165 |
term:=term + 1.0;
|
218 |
c := c + 1.0;
|
| 166 |
an := a*term;
|
219 |
an := a*c;
|
| 167 |
FOR i:=1 TO 2 DO
|
|
|
| 168 |
pn[i+4] := b*pn[i+2] - an*pn[i];
|
220 |
pn5 := b*pn3 - an*pn1;
|
| 169 |
END;
|
221 |
pn6 := b*pn4 - an*pn2;
|
| 170 |
IF (pn[6] # 0.0) THEN
|
222 |
IF (ABS(pn6) > 0.0) THEN
|
| 171 |
rn := pn[5] / pn[6];
|
223 |
rn := pn5/pn6;
|
| 172 |
dif := ABS(gin - rn);
|
224 |
IF (ABS(gamma-rn) <= DMIN(Tol,Tol*rn)) THEN
|
| 173 |
IF (dif <= eps) THEN (* absolute error tolerance satisfied *)
|
225 |
arg := arg + ln(gamma);
|
| 174 |
IF (dif <= eps*rn) THEN (* relative error tolerance satisfied *)
|
226 |
gamma := 1.0;
|
| 175 |
GInc := 1.0 - factor*gin;
|
227 |
IF (arg >= Elimit) THEN
|
|
|
228 |
gamma := 1.0 - exp(arg);
|
|
|
229 |
END; (* IF *)
|
| 176 |
EXIT;
|
230 |
EXIT
|
| 177 |
END;
|
|
|
| 178 |
END;
|
231 |
ELSE
|
| 179 |
gin := rn;
|
232 |
gamma := rn;
|
|
|
233 |
END; (* IF *)
|
| 180 |
END;
|
234 |
END; (* IF *)
|
| 181 |
FOR i:=1 TO 4 DO
|
|
|
| 182 |
pn[i] := pn[i+2];
|
235 |
pn1 := pn3;
|
| 183 |
END;
|
236 |
pn2 := pn4;
|
|
|
237 |
pn3 := pn5;
|
|
|
238 |
pn4 := pn6;
|
| 184 |
IF (ABS(pn[5]) >= Big) THEN
|
239 |
IF (ABS(pn5) >= Oflow) THEN
|
| 185 |
FOR i:=1 TO 4 DO
|
240 |
(* Re-scale terms in continued fraction if terms are large *)
|
| 186 |
pn[i] := pn[i] / Big;
|
241 |
pn1 := pn1 / Oflow;
|
| 187 |
END;
|
242 |
pn2 := pn2 / Oflow;
|
|
|
243 |
pn3 := pn3 / Oflow;
|
|
|
244 |
pn4 := pn4 / Oflow;
|
| 188 |
END;
|
245 |
END; (* IF *)
|
| 189 |
END; (* LOOP *)
|
246 |
END; (* LOOP *)
|
| 190 |
|
|
|
| 191 |
END; (* IF *)
|
247 |
END; (* IF *)
|
| 192 |
RETURN GInc;
|
248 |
RETURN gamma;
|
| 193 |
END GammaIn;
|
249 |
END GammaIn;
|
| 194 |
|
250 |
|
| 195 |
PROCEDURE GammaU(a,x : LONGREAL) : LONGREAL;
|
251 |
PROCEDURE GammaU(a,x : LONGREAL) : LONGREAL;
|
| 196 |
|
252 |
|
| 197 |
(*----------------------------------------------------------------*)
|
253 |
(*----------------------------------------------------------------*)
|