]>
Commit | Line | Data |
---|---|---|
fe4da5cc | 1 | * |
2 | * $Id$ | |
3 | * | |
4 | * $Log$ | |
5 | * Revision 1.1.1.1 1996/04/01 15:02:43 mclareni | |
6 | * Mathlib gen | |
7 | * | |
8 | * | |
9 | #include "gen/pilot.h" | |
10 | C This will be GAMDIS,IF=DOUBLE and GAMDIS64,IF=-DOUBLE | |
11 | FUNCTION GAMDIS(X,A) | |
12 | ||
13 | C Calculates the gamma distribution function | |
14 | C | |
15 | C G(x,a) = (1/gamma(a)) * int(0,x)[exp(-t) * t**(a-1)] dt. | |
16 | C | |
17 | C Based on | |
18 | C W. Gautschi, ALGORITHM 542 Incomplete Gamma Functions, | |
19 | C ACM Trans. Math. Software 5 (1979) 482-489. | |
20 | ||
21 | CHARACTER NAME*(*) | |
22 | CHARACTER*80 ERRTXT | |
23 | PARAMETER (NAME = 'GAMDIS') | |
24 | DIMENSION C(14) | |
25 | ||
26 | PARAMETER (EPS = 1E-5, EPS1 = 5E-7) | |
27 | PARAMETER (ALH = -0.69314 72) | |
28 | PARAMETER (Z1 = 1, HALF = Z1/2, QUAR = Z1/4) | |
29 | PARAMETER (C1 = 3*Z1/2, KMAX = 300) | |
30 | ||
31 | DATA C | |
32 | 1/ 0.5772157,-0.6558781,-0.0420026, 0.1665386,-0.0421977, | |
33 | 2 -0.0096220, 0.0072189,-0.0011652,-0.0002152, 0.0001281, | |
34 | 3 -0.0000201,-0.0000013, 0.0000011,-0.0000002/ | |
35 | ||
36 | HST=0 | |
37 | IF(X .EQ. 0) GO TO 99 | |
38 | IF(X .LT. 0 .OR. A .LE. 0) THEN | |
39 | WRITE(ERRTXT,101) X,A | |
40 | CALL MTLPRT(NAME,'G106.1',ERRTXT) | |
41 | GO TO 99 | |
42 | ELSE | |
43 | ALX=LOG(X) | |
44 | ENDIF | |
45 | IF(X .LT. QUAR) THEN | |
46 | ALFA=ALH/ALX | |
47 | ELSE | |
48 | ALFA=X+QUAR | |
49 | ENDIF | |
50 | IF(A .GT. ALFA) THEN | |
51 | TERM=1 | |
52 | SUM=1 | |
53 | DO 1 K = 1,KMAX | |
54 | TERM=X*TERM/(A+K) | |
55 | SUM=SUM+TERM | |
56 | IF(ABS(TERM) .LE. EPS*SUM) GO TO 2 | |
57 | 1 CONTINUE | |
58 | GO TO 98 | |
59 | 2 HST=SUM*EXP(A*ALX-X-ALOGAM(1+A)) | |
60 | ELSEIF(X .GT. C1) THEN | |
61 | P=0 | |
62 | S=1-A | |
63 | Q=(X+S)*(X-1-A) | |
64 | R=4*(X+S) | |
65 | TERM=1 | |
66 | SUM=1 | |
67 | RHO=0 | |
68 | DO 3 K = 2,KMAX | |
69 | P=P+S | |
70 | Q=Q+R | |
71 | R=R+8 | |
72 | S=S+2 | |
73 | T=P*(1+RHO) | |
74 | RHO=T/(Q-T) | |
75 | TERM=RHO*TERM | |
76 | SUM=SUM+TERM | |
77 | IF(ABS(TERM) .LE. EPS*SUM) GO TO 4 | |
78 | 3 CONTINUE | |
79 | GO TO 98 | |
80 | 4 HST=1-(A*SUM/(X+1-A))*EXP(A*ALX-X-ALOGAM(1+A)) | |
81 | ELSE | |
82 | IF(A .LT. HALF) THEN | |
83 | SUM=C(14) | |
84 | DO 12 K = 13,1,-1 | |
85 | 12 SUM=A*SUM+C(K) | |
86 | GA=-SUM/(1+A*SUM) | |
87 | Y=A*ALX | |
88 | IF(ABS(Y) .GE. 1) THEN | |
89 | U=GA-(EXP(Y)-1)/A | |
90 | ELSE | |
91 | SUM=1 | |
92 | TERM=1 | |
93 | DO 7 K = 2,KMAX | |
94 | TERM=Y*TERM/K | |
95 | SUM=SUM+TERM | |
96 | IF(ABS(TERM) .LE. EPS1*SUM) GO TO 8 | |
97 | 7 CONTINUE | |
98 | GO TO 98 | |
99 | 8 U=GA-SUM*ALX | |
100 | ENDIF | |
101 | ELSE | |
102 | U=GAMMA(A)-EXP(A*ALX)/A | |
103 | ENDIF | |
104 | P=A*X | |
105 | Q=A+1 | |
106 | R=A+3 | |
107 | TERM=1 | |
108 | SUM=1 | |
109 | DO 9 K = 2,KMAX | |
110 | P=P+X | |
111 | Q=Q+R | |
112 | R=R+2 | |
113 | TERM=-P*TERM/Q | |
114 | SUM=SUM+TERM | |
115 | IF(ABS(TERM) .LE. EPS1*SUM) GO TO 10 | |
116 | 9 CONTINUE | |
117 | GO TO 98 | |
118 | 10 HST=1-A*(U+SUM*EXP((1+A)*ALX)/(1+A))/GAMMA(1+A) | |
119 | ENDIF | |
120 | 99 GAMDIS=HST | |
121 | RETURN | |
122 | ||
123 | 98 WRITE(ERRTXT,102) X,A | |
124 | CALL MTLPRT(NAME,'G106.2',ERRTXT) | |
125 | GO TO 99 | |
126 | 101 FORMAT('ILLEGAL ARGUMENT(S) X = ',E15.8,' A = ',E15.8) | |
127 | 102 FORMAT('PROBLEMS WITH CONVERGENCE, X = ',E15.8,' A = ',E15.8) | |
128 | END |