]>
Commit | Line | Data |
---|---|---|
fe4da5cc | 1 | * |
2 | * $Id$ | |
3 | * | |
4 | * $Log$ | |
5 | * Revision 1.1.1.1 1996/04/01 15:02:05 mclareni | |
6 | * Mathlib gen | |
7 | * | |
8 | * | |
9 | #include "gen/pilot.h" | |
10 | #if !defined(CERNLIB_DOUBLE) | |
11 | FUNCTION RGAPNC(A,X) | |
12 | #endif | |
13 | #if defined(CERNLIB_DOUBLE) | |
14 | FUNCTION DGAPNC(A,X) | |
15 | #include "gen/imp64.inc" | |
16 | #endif | |
17 | C Calculates the incomplete gamma function P(A,X) as defined by | |
18 | C GSTAR(A,X) in Ref. 1. Based on | |
19 | C 1. W. Gautschi, ALGORITHM 542 Incomplete Gamma Functions, | |
20 | C ACM Trans. Math. Software 5 (1979) 482-489 | |
21 | C 2. W. Gautschi, A computational procedure for incomplete gamma | |
22 | C functions, ACM Trans. Math. Software 5 (1979) 466-481 | |
23 | ||
24 | CHARACTER NAME*(*) | |
25 | CHARACTER*80 ERRTXT | |
26 | #if !defined(CERNLIB_DOUBLE) | |
27 | PARAMETER (NAME = 'RGAPNC') | |
28 | #endif | |
29 | #if defined(CERNLIB_DOUBLE) | |
30 | PARAMETER (NAME = 'RGAPNC/DGAPNC') | |
31 | #endif | |
32 | PARAMETER (EPS = 5D-14) | |
33 | PARAMETER (ALH = -0.69314 71805 59945 31D0) | |
34 | PARAMETER (Z1 = 1, HALF = Z1/2, QUAR = Z1/4) | |
35 | PARAMETER (C1 = 3*Z1/2, KMAX = 600, EPS1 = EPS/100) | |
36 | ||
37 | DIMENSION C(25) | |
38 | ||
39 | DATA C( 1) / 0.57721 56649 01532 86D0/ | |
40 | DATA C( 2) /-0.65587 80715 20253 88D0/ | |
41 | DATA C( 3) /-0.04200 26350 34095 24D0/ | |
42 | DATA C( 4) / 0.16653 86113 82291 49D0/ | |
43 | DATA C( 5) /-0.04219 77345 55544 34D0/ | |
44 | DATA C( 6) /-0.00962 19715 27876 97D0/ | |
45 | DATA C( 7) / 0.00721 89432 46663 10D0/ | |
46 | DATA C( 8) /-0.00116 51675 91859 07D0/ | |
47 | DATA C( 9) /-0.00021 52416 74114 95D0/ | |
48 | DATA C(10) / 0.00012 80502 82388 12D0/ | |
49 | DATA C(11) /-0.00002 01348 54780 79D0/ | |
50 | DATA C(12) /-0.00000 12504 93482 14D0/ | |
51 | DATA C(13) / 0.00000 11330 27231 98D0/ | |
52 | DATA C(14) /-0.00000 02056 33841 70D0/ | |
53 | DATA C(15) / 0.00000 00061 16095 10D0/ | |
54 | DATA C(16) / 0.00000 00050 02007 64D0/ | |
55 | DATA C(17) /-0.00000 00011 81274 57D0/ | |
56 | DATA C(18) / 0.00000 00001 04342 67D0/ | |
57 | DATA C(19) / 0.00000 00000 07782 26D0/ | |
58 | DATA C(20) /-0.00000 00000 03696 81D0/ | |
59 | DATA C(21) / 0.00000 00000 00510 04D0/ | |
60 | DATA C(22) /-0.00000 00000 00020 58D0/ | |
61 | DATA C(23) /-0.00000 00000 00005 35D0/ | |
62 | DATA C(24) / 0.00000 00000 00001 23D0/ | |
63 | DATA C(25) /-0.00000 00000 00000 12D0/ | |
64 | ||
65 | #if !defined(CERNLIB_DOUBLE) | |
66 | GLGAMA(V)=ALGAMA(V) | |
67 | #endif | |
68 | #if defined(CERNLIB_DOUBLE) | |
69 | GLGAMA(V)=DLGAMA(V) | |
70 | #endif | |
71 | ||
72 | HST=0 | |
73 | IF(X .LT. 0) THEN | |
74 | WRITE(ERRTXT,101) X | |
75 | CALL MTLPRT(NAME,'C334.1',ERRTXT) | |
76 | GO TO 99 | |
77 | ELSEIF(X .EQ. 0) THEN | |
78 | IF(A .LT. 0) THEN | |
79 | WRITE(ERRTXT,102) A | |
80 | CALL MTLPRT(NAME,'C334.2',ERRTXT) | |
81 | ELSEIF(A .EQ. 0) THEN | |
82 | HST=1 | |
83 | ENDIF | |
84 | GO TO 99 | |
85 | ELSE | |
86 | ALX=LOG(X) | |
87 | ENDIF | |
88 | IF(X .LT. QUAR) THEN | |
89 | ALFA=ALH/ALX | |
90 | ELSE | |
91 | ALFA=X+QUAR | |
92 | ENDIF | |
93 | MA=HALF-A | |
94 | AEPS=A+MA | |
95 | ||
96 | IF(MA .GT. 0) THEN | |
97 | IF(AEPS .NE. 0) THEN | |
98 | SG=(-1)**(MA-1)*SIGN(Z1,A)*SIGN(Z1,AEPS) | |
99 | ALGP1=GLGAMA(1+AEPS)-LOG(ABS(AEPS)) | |
100 | IF(MA .NE. 1) ALGP1=ALGP1+GLGAMA(1-AEPS)-GLGAMA(MA-AEPS) | |
101 | ELSE | |
102 | SG=0 | |
103 | ALGP1=0 | |
104 | ENDIF | |
105 | ELSE | |
106 | SG=SIGN(Z1,A) | |
107 | ALGP1=GLGAMA(1+A) | |
108 | ENDIF | |
109 | IF(A .GT. ALFA) THEN | |
110 | TERM=1 | |
111 | SUM=1 | |
112 | DO 1 K = 1,KMAX | |
113 | TERM=X*TERM/(A+K) | |
114 | SUM=SUM+TERM | |
115 | IF(ABS(TERM) .LE. EPS*SUM) GO TO 2 | |
116 | 1 CONTINUE | |
117 | GO TO 98 | |
118 | 2 HST=EXP(A*ALX-X+LOG(SUM)-ALGP1) | |
119 | ELSEIF(X .GT. C1) THEN | |
120 | P=0 | |
121 | S=1-A | |
122 | Q=(X+S)*(X-1-A) | |
123 | R=4*(X+S) | |
124 | TERM=1 | |
125 | SUM=1 | |
126 | RHO=0 | |
127 | DO 3 K = 2,KMAX | |
128 | P=P+S | |
129 | Q=Q+R | |
130 | R=R+8 | |
131 | S=S+2 | |
132 | T=P*(1+RHO) | |
133 | RHO=T/(Q-T) | |
134 | TERM=RHO*TERM | |
135 | SUM=SUM+TERM | |
136 | IF(ABS(TERM) .LE. EPS*SUM) GO TO 4 | |
137 | 3 CONTINUE | |
138 | GO TO 98 | |
139 | 4 IF(A .LT. 0) THEN | |
140 | HST=1 | |
141 | IF(MA .LT. 0 .OR. AEPS .NE. 0) | |
142 | 1 HST=1-SG*EXP(LOG(ABS(A)*SUM/(X+1-A))-X+A*ALX-ALGP1) | |
143 | ELSEIF(A .EQ. 0) THEN | |
144 | HST=1 | |
145 | ELSE | |
146 | HST=1-EXP(A*ALX-X+LOG(A*SUM/(X+1-A))-ALGP1) | |
147 | ENDIF | |
148 | ELSE | |
149 | AE=A | |
150 | IF(A .LT. HALF) THEN | |
151 | IF(A .LT. -HALF) AE=AEPS | |
152 | SUM=C(25) | |
153 | DO 12 K = 24,1,-1 | |
154 | 12 SUM=AE*SUM+C(K) | |
155 | GA=-SUM/(1+AE*SUM) | |
156 | Y=AE*ALX | |
157 | IF(ABS(Y) .GE. 1) THEN | |
158 | U=GA-(EXP(Y)-1)/AE | |
159 | ELSE | |
160 | SUM=1 | |
161 | TERM=1 | |
162 | DO 7 K = 2,KMAX | |
163 | TERM=Y*TERM/K | |
164 | SUM=SUM+TERM | |
165 | IF(ABS(TERM) .LE. EPS1*SUM) GO TO 8 | |
166 | 7 CONTINUE | |
167 | GO TO 98 | |
168 | 8 U=GA-SUM*ALX | |
169 | ENDIF | |
170 | ELSE | |
171 | U=EXP(GLGAMA(A))-X**A/A | |
172 | ENDIF | |
173 | P=AE*X | |
174 | Q=AE+1 | |
175 | R=AE+3 | |
176 | TERM=1 | |
177 | SUM=1 | |
178 | DO 9 K = 2,KMAX | |
179 | P=P+X | |
180 | Q=Q+R | |
181 | R=R+2 | |
182 | TERM=-P*TERM/Q | |
183 | SUM=SUM+TERM | |
184 | IF(ABS(TERM) .LE. EPS1*SUM) GO TO 10 | |
185 | 9 CONTINUE | |
186 | GO TO 98 | |
187 | 10 H=U+X**(AE+1)*SUM/(AE+1) | |
188 | IF(A .LT. -HALF) THEN | |
189 | H=H*EXP(X-AE*ALX) | |
190 | DO 13 J = 1,MA | |
191 | 13 H=(1-X*H)/(J-AE) | |
192 | HST=1 | |
193 | IF(MA .LT. 0 .OR. AEPS .NE. 0) | |
194 | 1 HST=1-SG*EXP(LOG(ABS(A)*H)-X+A*ALX-ALGP1) | |
195 | ELSEIF(A .EQ. 0) THEN | |
196 | HST=1 | |
197 | ELSE | |
198 | HST=1-A*H*EXP(-ALGP1) | |
199 | ENDIF | |
200 | ENDIF | |
201 | #if defined(CERNLIB_DOUBLE) | |
202 | 99 DGAPNC=HST | |
203 | #endif | |
204 | #if !defined(CERNLIB_DOUBLE) | |
205 | 99 RGAPNC=HST | |
206 | #endif | |
207 | RETURN | |
208 | ||
209 | 98 WRITE(ERRTXT,103) A,X | |
210 | CALL MTLPRT(NAME,'C334.3',ERRTXT) | |
211 | GO TO 99 | |
212 | 101 FORMAT('ILLEGAL ARGUMENT X = ',1P,D15.8,' < 0') | |
213 | 102 FORMAT('ILLEGAL ARGUMENTS A = ',1P,D15.8,' < 0, X = 0') | |
214 | 103 FORMAT('PROBLEMS WITH CONVERGENCE, A = ',1P,D15.8,' X = ',D15.8) | |
215 | END |