      REAL FUNCTION GAMMA(X)
C***BEGIN PROLOGUE  GAMMA
C***DATE WRITTEN   770601   (YYMMDD)
C***REVISION DATE  820801   (YYMMDD)
C***CATEGORY NO.  C7A
C***KEYWORDS  GAMMA FUNCTION,SPECIAL FUNCTION
C***AUTHOR  FULLERTON, W., (LANL)
C***PURPOSE  Computes the Gamma function.
C***DESCRIPTION
C
C GAMMA computes the gamma function at X, where X is not 0, -1, -2, ....
C GAMMA and X are single precision.
C***REFERENCES  (NONE)
C***ROUTINES CALLED  CSEVL,GAMLIM,INITS,R1MACH,R9LGMC,XERROR
C***END PROLOGUE  GAMMA
      INTEGER NGCS,INITS,N,I
      REAL GCS,PI,SQ2PIL,XMIN,XMAX,DXREL,R1MACH,X,Y,CSEVL
      REAL R9LGMC,SINPIY
      DIMENSION GCS(23)
      DATA GCS   ( 1) / .0085711955 90989331E0/
      DATA GCS   ( 2) / .0044153813 24841007E0/
      DATA GCS   ( 3) / .0568504368 1599363E0/
      DATA GCS   ( 4) /-.0042198353 96418561E0/
      DATA GCS   ( 5) / .0013268081 81212460E0/
      DATA GCS   ( 6) /-.0001893024 529798880E0/
      DATA GCS   ( 7) / .0000360692 532744124E0/
      DATA GCS   ( 8) /-.0000060567 619044608E0/
      DATA GCS   ( 9) / .0000010558 295463022E0/
      DATA GCS   (10) /-.0000001811 967365542E0/
      DATA GCS   (11) / .0000000311 772496471E0/
      DATA GCS   (12) /-.0000000053 542196390E0/
      DATA GCS   (13) / .0000000009 193275519E0/
      DATA GCS   (14) /-.0000000001 577941280E0/
      DATA GCS   (15) / .0000000000 270798062E0/
      DATA GCS   (16) /-.0000000000 046468186E0/
      DATA GCS   (17) / .0000000000 007973350E0/
      DATA GCS   (18) /-.0000000000 001368078E0/
      DATA GCS   (19) / .0000000000 000234731E0/
      DATA GCS   (20) /-.0000000000 000040274E0/
      DATA GCS   (21) / .0000000000 000006910E0/
      DATA GCS   (22) /-.0000000000 000001185E0/
      DATA GCS   (23) / .0000000000 000000203E0/
      DATA PI /3.14159 26535 89793 24E0/
C SQ2PIL IS ALOG (SQRT (2.*PI) )
      DATA SQ2PIL /0.91893 85332 04672 74E0/
      DATA NGCS, XMIN, XMAX, DXREL /0, 3*0.0 /
C
C
C***FIRST EXECUTABLE STATEMENT  GAMMA
      IF (NGCS.NE.0) GO TO 10 
C
C ---------------------------------------------------------------------
C INITIALIZE.  FIND LEGAL BOUNDS FOR X, AND DETERMINE THE NUMBER OF
C TERMS IN THE SERIES REQUIRED TO ATTAIN AN ACCURACY TEN TIMES BETTER 
C THAN MACHINE PRECISION.
C
      NGCS = INITS (GCS, 23, 0.1*R1MACH(3))
C
      CALL GAMLIM (XMIN, XMAX)
      DXREL = SQRT (R1MACH(4))
C
C ---------------------------------------------------------------------
C FINISH INITIALIZATION.  START EVALUATING GAMMA(X).
C
 10   Y = ABS(X)
      IF (Y.GT.10.0) GO TO 50 
C
C COMPUTE GAMMA(X) FOR ABS(X) .LE. 10.0.  REDUCE INTERVAL AND
C FIND GAMMA(1+Y) FOR 0. .LE. Y .LT. 1. FIRST OF ALL.
C
      N = X
      IF (X.LT.0.) N = N - 1
      Y = X - FLOAT(N)
      N = N - 1
      GAMMA = 0.9375 + CSEVL(2.*Y-1., GCS, NGCS)
      IF (N.EQ.0) RETURN
C
      IF (N.GT.0) GO TO 30
C
C COMPUTE GAMMA(X) FOR X .LT. 1.
C
      N = -N
      IF (X.EQ.0.) CALL XERROR ( 'GAMMA   X IS 0', 14, 4, 2)
      IF (X.LT.0. .AND. X+FLOAT(N-2).EQ.0.) CALL XERROR (  'GAMMA   X IS
     1 A NEGATIVE INTEGER', 31, 4, 2)
      IF (X.LT.(-0.5) .AND. ABS((X-AINT(X-0.5))/X).LT.DXREL) CALL
     1  XERROR ( 'GAMMA   ANSWER LT HALF PRECISION BECAUSE X TOO NEAR NE
     2GATIVE INTEGER', 68, 1, 1)
C
      DO 20 I=1,N
        GAMMA = GAMMA / (X+FLOAT(I-1))
 20   CONTINUE
      RETURN
C
C GAMMA(X) FOR X .GE. 2.
C
 30   DO 40 I=1,N
        GAMMA = (Y+FLOAT(I))*GAMMA
 40   CONTINUE
      RETURN
C
C COMPUTE GAMMA(X) FOR ABS(X) .GT. 10.0.  RECALL Y = ABS(X).
C
 50   IF (X.GT.XMAX) CALL XERROR ( 'GAMMA   X SO BIG GAMMA OVERFLOWS',
     1  32, 3, 2)
C
      GAMMA = 0.
      IF (X.LT.XMIN) CALL XERROR ( 'GAMMA   X SO SMALL GAMMA UNDERFLOWS'
     1  , 35, 2, 1) 
      IF (X.LT.XMIN) RETURN
C
      GAMMA = EXP((Y-0.5)*ALOG(Y) - Y + SQ2PIL + R9LGMC(Y) )
      IF (X.GT.0.) RETURN
C
      IF (ABS((X-AINT(X-0.5))/X).LT.DXREL) CALL XERROR ( 'GAMMA   ANSWER
     1 LT HALF PRECISION, X TOO NEAR NEGATIVE INTEGER'  , 61, 1, 1)
C
      SINPIY = SIN (PI*Y)
      IF (SINPIY.EQ.0.) CALL XERROR ( 'GAMMA   X IS A NEGATIVE INTEGER',
     1 31, 4, 2)
C
      GAMMA = -PI / (Y*SINPIY*GAMMA)
C
      RETURN
      END 
