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