This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / MINICERN / mathlib / gen / g / gamdis.F
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