C This file: http://ftp.aset.psu.edu/pub/ger/fortran/hdk/gamma.for
C
C===Sample (interactive) driver to test DGAMMA function.
C
      DOUBLE PRECISION X,DGAMMA,VAL,ARG
      REAL Y
      LOGICAL ERROR
      EXTERNAL DGAMMA
      DO 1 I=1,999999999
        WRITE(*,*) 'Give argument of GAMMA (or Zero) to quit:'
        READ(*,*,END=88) ARG
        IF(ARG.EQ.0) GOTO 88
        VAL=ARG
        X=DGAMMA(VAL,15,ERROR)
        IF(ERROR) THEN
        WRITE(*,*) '$$GAMMA Argument ',ARG,' is <= 0 or > 171.'
        ELSE
          WRITE(*,*) 'GAMMA(',ARG,')=',X
          WRITE(*,2) 'GAMMA(',ARG,')=',X
2         FORMAT(A,F5.3,A,F18.0)
        ENDIF
1     CONTINUE
88    STOP
      END
      DOUBLE PRECISION FUNCTION DGAMMA(Z,MSIZE,ERROR)
C=== ACM CALGO Algorithm 309, Double Precision Gamma Function.
C    Z is the argument; must be > 0 and < 57 for IBM mainframes, or
C                                       < 171 for IBM PC.
C    MSIZE is the number of decimal digits representable; 15 for both.
C    ERROR is returned LOGICAL; .TRUE. if argument was out of range, or
C                               .FALSE. otherwise (i.e. no error).
C    Translated from Algol 68, March 6, 1989, HDK.
C
      DOUBLE PRECISION X,Y,Z,PI, DELTA,EX,DECIMA,MAXZ,ARGUND,MAXEXP
      DOUBLE PRECISION DLOGAM
      INTEGER MSIZE,PARITY,ENTIER
      LOGICAL ERROR
C---IBM Mainframes:
C     DATA MAXZ/57.D0/, ARGUND/-1.D-74/, MAXEXP/174.D0/
C---IBM PC.
      DATA MAXZ/171.D0/, ARGUND/-1.D-307/, MAXEXP/702.D0/
      DATA PI/3.1415926535897932384626433832795028841971693993751D0/
C     ENTIER(Z)=FLOOR(Z)
      ENTIER(X)=DNINT(X)-DMOD(2.D0+DSIGN(1.D0,X),3.D0)
      PARITY(Z)=MOD(ENTIER(Z),2)
      DECIMA(Z)=DABS(Z-IDINT(Z))
      IF(Z.GT.MAXZ) THEN
        GOTO 99
      ELSE
        IF(Z.EQ.ENTIER(Z)) THEN
          IF(Z.LE.0) THEN
            GOTO 99
          ELSE
            Y=1.D0
            IF(Z.GT.2.D0) THEN
1             Z=Z-1
              Y=Y*Z
              IF(Z.GT.2.D0) GOTO 1
            ENDIF
          ENDIF
        ELSE
          IF(DABS(Z).LT.10**(-MSIZE)) THEN
            Y=1.D0/Z
          ELSE
            IF(Z.LT.0) THEN
              IF(Z.LT.ARGUND) THEN
                Y=0
              ELSE
                DELTA=DECIMA(Z)*PI
                IF(DELTA.LT.10**(-MSIZE*0.5D0)) THEN
                  EX=-DLOG(DECIMA(Z))
                ELSE
                  EX=DLOG(PI/DSIN(DELTA)) - DLOGAM(1.D0-Z, MSIZE)
                ENDIF
                IF(DABS(EX).GT.MAXEXP) THEN
                  Y=0
                ELSE
                  IF(PARITY(Z).EQ.1) THEN
                    Y=DEXP(EX)
                  ELSE
                    Y=-DEXP(EX)
                  ENDIF
                ENDIF
              ENDIF
            ELSE
              Y=DEXP(DLOGAM(Z,MSIZE))
            ENDIF
          ENDIF
        ENDIF
      ENDIF
      ERROR=.FALSE.
      DGAMMA=Y
      RETURN
99    ERROR=.TRUE.
      DGAMMA=0
      RETURN
      END
      DOUBLE PRECISION FUNCTION DLOGAM(T,MSIZE)
C===Log base e of GAMMA(T). See DGAMMA.
      DOUBLE PRECISION F,T,DLGM
      INTEGER TMIN,MSIZE
      TMIN=7
      IF(MSIZE.GT.18) TMIN=MSIZE-10
      IF(T.GT.TMIN) THEN
        DLOGAM=DLGM(T)
      ELSE
        F=T
1       T=T+1
        IF(T.LT.TMIN) THEN
          F=F*T
          GOTO 1
        ENDIF
        DLOGAM=DLGM(T)-DLOG(F)
      ENDIF
      RETURN
      END
      DOUBLE PRECISION FUNCTION DLGM(W)
C===Companion function for DGAMMA and DLOGAM.
      DOUBLE PRECISION C(20), W,W2,PRESUM,CONST,DEN,SUM
      INTEGER I
      C(1)=1.D0/12.D0
      C(2)=-1.D0/360.D0
      C(3)=1.D0/1260.D0
      C(4)=-1.D0/1680.D0
      C(5)=1.D0/1188.D0
      C(6)=-691.D0/360360.D0
      C(7)=1.D0/156.D0
      C(8)=-3617.D0/122400.D0
      C(9)=43867.D0/244188.D0
      C(10)=-174611.D0/125400.D0
      C(11)=77683.D0/5796.D0
      C(12)=-236364091.D0/1506960.D0
      C(13)=657931.D0/300.D0
      C(14)=-3392780147.D0/93960.D0
      C(15)=1723168255201.D0/2492028.D0
      C(16)=-7709321041217.D0/505920.D0
      C(17)=1516286977551.D0/396.D0
      C(18)=-26315271553053477373.D0/2418179400.D0
      C(19)=1542102059911661.D0/444.D0
      C(20)=-261082718496449122051.D0/21106800.D0
      CONST=.91893853320467274178032973640561763986139747363778D0
C---CONST=DLOG(SQRT(2*PI))
      DEN=W
      W2=W*W
      PRESUM=(W-.5D0)*DLOG(W)-W+CONST
      DO 1 I=1,20
        SUM=PRESUM+C(I)/DEN
        IF(SUM.EQ.PRESUM) GOTO 2
        DEN=DEN*W2
1     PRESUM=SUM
2     DLGM=SUM
      RETURN
      END