      FUNCTION GAMMA(A)
C     ******************************************************************
C     EVALUATION OF THE GAMMA FUNCTION FOR REAL ARGUMENTS
C     ******************************************************************
      REAL P(7), Q(7)
      DOUBLE PRECISION D, G, Z, LNX, GLOG
C     ---------------------
C     WRITTEN BY ALFRED H. MORRIS
C          NAVAL SURFACE WEAPONS CENTER
C          DAHLGREN, VIRGINIA
C     ---------------------
C     AMAX IS THE LARGEST POSITIVE VALUE OF A FOR WHICH GAMMA(A) IS
C     COMPUTED. THIS IS A MACHINE DEPENDENT CONSTANT. IF B IS THE
C     BASE OF THE FLOATING POINT ARITHMETIC BEING USED AND EMAX THE
C     SMALLEST INTEGER FOR WHICH X .LT. B**EMAX FOR ALL POSITIVE
C     REAL X (SEE THE FUNCTION I1MACH), THEN THE FOLLOWING SETTINGS
C     ARE APPROPRIATE ...
C
C               B       EMAX        AMAX
C
C               2        127        34.8
C               8         76        53.4
C              16         63        57.5
C               2        255        58.0
C               2       1070       176.5
C               2       8191       966.9
C
C     IN PARTICULAR, THE FOLLOWING SETTING IS APPROPRIATE FOR THE
C     CDC 6000-7000 SERIES COMPUTERS ...
C
      DATA AMAX/176.5/
C
C     GAMMA(A) IS ASSIGNED THE VALUE 0 WHEN THE GAMMA FUNCTION
C     CANNOT BE COMPUTED. IT IS ASSUMED THAT 1/T DOES NOT OVERFLOW
C     FOR ANY NONZERO T.
C     ---------------------
C     D = 0.5*(LN(2*PI) - 1)
C     ---------------------
      DATA PI/3.1415926535898/
      DATA D/.41893853320467274178D0/
C     ---------------------
      DATA P(1)/ .539637273585445E-03/,  P(2)/ .261939260042690E-02/,
     1     P(3)/ .204493667594920E-01/,  P(4)/ .730981088720487E-01/,
     2     P(5)/ .279648642639792E+00/,  P(6)/ .553413866010467E+00/,
     3     P(7)/ 1.0/
      DATA Q(1)/-.832979206704073E-03/,  Q(2)/ .470059485860584E-02/,
     1     Q(3)/ .225211131035340E-01/,  Q(4)/-.170458969313360E+00/,
     2     Q(5)/-.567902761974940E-01/,  Q(6)/ .113062953091122E+01/,
     3     Q(7)/ 1.0/
C     ---------------------
      DATA R1/.820756370353826E-03/, R2/-.595156336428591E-03/,
     1     R3/.793650663183693E-03/, R4/-.277777777770481E-02/,
     2     R5/.833333333333333E-01/
C     ---------------------
      GAMMA = 0.0
      X = A
      IF (ABS(A) .GE. 15.0) GO TO 40
C     ------------------------------------------------------------------
C          EVALUATION OF GAMMA(A) FOR ABS(A).LT.15
C     ------------------------------------------------------------------
      T = 1.0
      M = INT(A) - 1
C
C     LET T BE THE PRODUCT OF A-J WHEN A.GE.2
C
      IF (M) 20,12,10
   10 DO 11 J = 1,M
        X = X - 1.0
   11   T = X*T
   12 X = X - 1.0
      GO TO 30
C
C     LET T BE THE PRODUCT OF A+J WHEN A.LT.1
C
   20 T = A
      IF (A) 21,60,30
   21 M = -M
      IF (M .EQ. 1) GO TO 23
      DO 22 J = 1,M
        X = X + 1.0
   22   T = X*T
      IF (T) 30,60,30
   23 X = (A + 0.5) + 0.5
      T = X*T
      IF (T .EQ. 0.0) RETURN
C
C     COMPUTE GAMMA(1 + X) FOR 0.LE.X.LT.1
C
   30 TOP = P(1)
      BOT = Q(1)
      DO 31 I = 2,7
        TOP = P(I) + X*TOP
   31   BOT = Q(I) + X*BOT
      GAMMA = TOP/BOT
C
C     TERMINATION
C
      IF (A .LT. 1.0) GO TO 35
      GAMMA = GAMMA*T
      RETURN
   35 GAMMA = GAMMA/T
      RETURN
C     ------------------------------------------------------------------
C          EVALUATION OF GAMMA(A) FOR ABS(A).GE.15
C     ------------------------------------------------------------------
   40 IF (A .GT. 0.0) GO TO 50
      IF (A .LT. 1.0 - AMAX) RETURN
      X = -A
      N = X
      T = X - N
      IF (T .GT. 0.9) T = 1.0 - T
      S = SIN(PI*T)/PI
      IF (MOD(N,2) .EQ. 0) S = -S
      IF (S .EQ. 0.0) RETURN
C
C     COMPUTE THE MODIFIED ASYMPTOTIC SUM
C
   50 IF (A .GT. AMAX) RETURN
      T = 1.0/(X*X)
      G = ((((R1*T + R2)*T + R3)*T + R4)*T + R5)/X
C
C     FINAL ASSEMBLY
C
      LNX = GLOG(X)
      Z = X
      G = (D + G) + (Z - 0.5D0)*(LNX - 1.D0)
      GAMMA = G
      T = G - DBLE(GAMMA)
      GAMMA = EXP(GAMMA)*(1.0 + T)
      IF (A .LT. 0.0) GAMMA = 1.0/(S*X*GAMMA)
   60 RETURN
      END
