SUBROUTINE GAMCDF(X,GAMMA,CDF)
C
C     PURPOSE--THIS SUBROUTINE COMPUTES THE CUMULATIVE DISTRIBUTION
C              FUNCTION VALUE FOR THE GAMMA
C              DISTRIBUTION WITH SINGLE PRECISION 
C              TAIL LENGTH PARAMETER = GAMMA.
C              THE GAMMA DISTRIBUTION USED
C              HEREIN HAS MEAN = GAMMA
C              AND STANDARD DEVIATION = SQRT(GAMMA).
C              THIS DISTRIBUTION IS DEFINED FOR ALL POSITIVE X,
C              AND HAS THE PROBABILITY DENSITY FUNCTION
C              F(X) = (1/CONSTANT) * (X**(GAMMA-1)) * EXP(-X)
C              WHERE THE CONSTANT = THE GAMMA FUNCTION EVALUATED
C              AT THE VALUE GAMMA.
C     INPUT  ARGUMENTS--X      = THE SINGLE PRECISION VALUE 
C                                AT WHICH THE CUMULATIVE DISTRIBUTION 
C                                FUNCTION IS TO BE EVALUATED.
C                                X SHOULD BE POSITIVE.
C                     --GAMMA  = THE SINGLE PRECISION VALUE 
C                                OF THE TAIL LENGTH PARAMETER.
C                                GAMMA SHOULD BE POSITIVE.
C     OUTPUT ARGUMENTS--CDF    = THE SINGLE PRECISION CUMULATIVE
C                                DISTRIBUTION FUNCTION VALUE.
C     OUTPUT--THE SINGLE PRECISION CUMULATIVE DISTRIBUTION
C             FUNCTION VALUE CDF FOR THE GAMMA DISTRIBUTION 
C             WITH TAIL LENGTH PARAMETER VALUE = GAMMA.
C     PRINTING--NONE UNLESS AN INPUT ARGUMENT ERROR CONDITION EXISTS. 
C     RESTRICTIONS--GAMMA SHOULD BE POSITIVE.
C                 --X SHOULD BE POSITIVE.
C     OTHER DATAPAC   SUBROUTINES NEEDED--NONE.
C     FORTRAN LIBRARY SUBROUTINES NEEDED--DEXP, DLOG.
C     MODE OF INTERNAL OPERATIONS--DOUBLE PRECISION.
C     LANGUAGE--ANSI FORTRAN. 
C     ACCURACY--(ON THE UNIVAC 1108, EXEC 8 SYSTEM AT NBS)
C               COMPARED TO THE KNOWN GAMMA = 1 (EXPONENTIAL)
C               RESULTS, AGREEMENT WAS HAD OUT TO 7 SIGNIFICANT
C               DIGITS FOR ALL TESTED X.
C               THE TESTED X VALUES COVERED THE ENTIRE
C               RANGE OF THE DISTRIBUTION--FROM THE 0.00001 
C               PERCENT POINT UP TO THE 99.99999 PERCENT POINT
C               OF THE DISTRIBUTION.
C     REFERENCES--WILK, GNANADESIKAN, AND HUYETT, 'PROBABILITY
C                 PLOTS FOR THE GAMMA DISTRIBUTION',
C                 TECHNOMETRICS, 1962, PAGES 1-15,
C                 ESPECIALLY PAGES 3-5. 
C               --NATIONAL BUREAU OF STANDARDS APPLIED MATHEMATICS
C                 SERIES 55, 1964, PAGE 257, FORMULA 6.1.41.
C               --JOHNSON AND KOTZ, CONTINUOUS UNIVARIATE
C                 DISTRIBUTIONS--1, 1970, PAGES 166-206.
C               --HASTINGS AND PEACOCK, STATISTICAL
C                 DISTRIBUTIONS--A HANDBOOK FOR
C                 STUDENTS AND PRACTITIONERS, 1975,
C                 PAGES 68-73.
C     WRITTEN BY--JAMES J. FILLIBEN
C                 STATISTICAL ENGINEERING LABORATORY (205.03)
C                 NATIONAL BUREAU OF STANDARDS
C                 WASHINGTON, D. C. 20234
C                 PHONE:  301-921-2315
C     ORIGINAL VERSION--NOVEMBER  1975. 
C
C---------------------------------------------------------------------
C
      DOUBLE PRECISION DX,DGAMMA,AI,TERM,SUM,CUT1,CUT2,CUTOFF,T
      DOUBLE PRECISION Z,Z2,Z3,Z4,Z5,DEN,A,B,C,D,G
      DOUBLE PRECISION DEXP,DLOG
      DIMENSION D(10)
      DATA C/ .918938533204672741D0/
      DATA D(1),D(2),D(3),D(4),D(5)
     1                 /+.833333333333333333D-1,-.277777777777777778D-2,
     1+.793650793650793651D-3,-.595238095238095238D-3,+.8417508417508417
     151D-3/
      DATA D(6),D(7),D(8),D(9),D(10)
     1     /-.191752691752691753D-2,+.641025641025641025D-2,-.2955065359
     147712418D-1,+.179644372368830573D0,-.139243221690590111D1/
C
C
      IPR=6
C
C     CHECK THE INPUT ARGUMENTS FOR ERRORS
C
      IF(X.LE.0.0)GOTO50
      IF(GAMMA.LE.0.0)GOTO55
      GOTO90
   50 WRITE(IPR,4)
      WRITE(IPR,46)X
      CDF=0.0
      RETURN
   55 WRITE(IPR,15) 
      WRITE(IPR,46)GAMMA
      CDF=0.0
      RETURN
   90 CONTINUE
    4 FORMAT(1H ,100H***** NON-FATAL DIAGNOSTIC--THE FIRST  INPUT ARGUME
     1NT TO THE GAMCDF SUBROUTINE IS NON-POSITIVE *****)
   15 FORMAT(1H , 91H***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE
     1 GAMCDF SUBROUTINE IS NON-POSITIVE *****)
   46 FORMAT(1H , 35H***** THE VALUE OF THE ARGUMENT IS ,E15.8,6H *****)
C
C-----START POINT-----------------------------------------------------
C
      DX=X
      DGAMMA=GAMMA
      MAXIT=10000
C
C     COMPUTE THE GAMMA FUNCTION USING THE ALGORITHM IN THE 
C     NBS APPLIED MATHEMATICS SERIES REFERENCE.
C
      Z=DGAMMA
      DEN=1.0D0
  300 IF(Z.GE.10.0D0)GOTO400
      DEN=DEN*Z
      Z=Z+1
      GOTO300
  400 Z2=Z*Z
      Z3=Z*Z2
      Z4=Z2*Z2
      Z5=Z2*Z3
      A=(Z-0.5D0)*DLOG(Z)-Z+C 
      B=D(1)/Z+D(2)/Z3+D(3)/Z5+D(4)/(Z2*Z5)+D(5)/(Z4*Z5)+
     1D(6)/(Z*Z5*Z5)+D(7)/(Z3*Z5*Z5)+D(8)/(Z5*Z5*Z5)+D(9)/(Z2*Z5*Z5*Z5)
      G=DEXP(A+B)/DEN
C
C     COMPUTE T-SUB-Q AS DEFINED ON PAGE 4 OF THE WILK, GNANADESIKAN, 
C     AND HUYETT REFERENCE
C
      SUM=1.0D0/DGAMMA
      TERM=1.0D0/DGAMMA
      CUT1=DX-DGAMMA
      CUT2=DX*10000000000.0D0 
      DO200I=1,MAXIT
      AI=I
      TERM=DX*TERM/(DGAMMA+AI)
      SUM=SUM+TERM
      CUTOFF=CUT1+(CUT2*TERM/SUM)
      IF(AI.GT.CUTOFF)GOTO250 
  200 CONTINUE
      WRITE(IPR,205)MAXIT
      WRITE(IPR,206)X
      WRITE(IPR,207)GAMMA
      WRITE(IPR,208)
      CDF=1.0
      RETURN
C
  250 T=SUM
      CDF=(DX**DGAMMA)*(DEXP(-DX))*T/G
C
  205 FORMAT(1H ,48H*****ERROR IN INTERNAL OPERATIONS IN THE GAMCDF , 
     1 45HSUBROUTINE--THE NUMBER OF ITERATIONS EXCEEDS ,I7) 
  206 FORMAT(1H ,33H     THE INPUT VALUE OF X     IS ,E15.8)
  207 FORMAT(1H ,33H     THE INPUT VALUE OF GAMMA IS ,E15.8)
  208 FORMAT(1H ,48H     THE OUTPUT VALUE OF CDF HAS BEEN SET TO 1.0) 
C
      RETURN
      END