      SUBROUTINE GAMLIM(XMIN,XMAX)
C***BEGIN PROLOGUE  GAMLIM
C***DATE WRITTEN   770401   (YYMMDD)
C***REVISION DATE  820801   (YYMMDD)
C***CATEGORY NO.  C7A,R2
C***KEYWORDS  GAMMA FUNCTION,LIMITS,SPECIAL FUNCTION
C***AUTHOR  FULLERTON, W., (LANL)
C***PURPOSE  Computes the minimum and maximum bounds for X in GAMMA(X).
C***DESCRIPTION
C
C Calculate the minimum and maximum legal bounds for X in GAMMA(X).
C XMIN and XMAX are not the only bounds, but they are the only non-
C trivial ones to calculate.
C
C             Output Arguments --
C XMIN   minimum legal value of X in GAMMA(X).  Any smaller value of
C        X might result in underflow.
C XMAX   maximum legal value of X in GAMMA(X).  Any larger value will
C        cause overflow.
C***REFERENCES  (NONE)
C***ROUTINES CALLED  R1MACH,XERROR
C***END PROLOGUE  GAMLIM
C***FIRST EXECUTABLE STATEMENT  GAMLIM
      ALNSML = ALOG(R1MACH(1))
      XMIN = -ALNSML
      DO 10 I=1,10
        XOLD = XMIN
        XLN = ALOG(XMIN)
        XMIN = XMIN - XMIN*((XMIN+0.5)*XLN - XMIN - 0.2258 + ALNSML)
     1    / (XMIN*XLN + 0.5)
        IF (ABS(XMIN-XOLD).LT.0.005) GO TO 20
10    CONTINUE
      CALL XERROR ( 'GAMLIM  UNABLE TO FIND XMIN', 27, 1, 2)
C
20    XMIN = -XMIN + 0.01
C
      ALNBIG = ALOG(R1MACH(2))
      XMAX = ALNBIG
      DO 30 I=1,10
        XOLD = XMAX
        XLN = ALOG(XMAX)
        XMAX = XMAX - XMAX*((XMAX-0.5)*XLN - XMAX + 0.9189 - ALNBIG)
     1    / (XMAX*XLN - 0.5)
        IF (ABS(XMAX-XOLD).LT.0.005) GO TO 40
30    CONTINUE
      CALL XERROR ( 'GAMLIM  UNABLE TO FIND XMAX', 27, 2, 2)
C
40    XMAX = XMAX - 0.01
      XMIN = AMAX1 (XMIN, -XMAX+1.)
C
      RETURN
      END
      FUNCTION INITS(OS,NOS,ETA)
C***BEGIN PROLOGUE  INITS
C***DATE WRITTEN   770401   (YYMMDD)
C***REVISION DATE  820801   (YYMMDD)
C***CATEGORY NO.  C3A2
C***KEYWORDS  INITIALIZE,ORTHOGONAL SERIES,SPECIAL FUNCTION
C***AUTHOR  FULLERTON, W., (LANL)
C***PURPOSE  Initializes an orthogonal series so that it defines the
C            number of terms to carry in the series to meet a specified
C            error.
C***DESCRIPTION
C
C Initialize the orthogonal series so that INITS is the number of terms
C needed to insure the error is no larger than ETA.  Ordinarily, ETA
C will be chosen to be one-tenth machine precision.
C
C             Input Arguments --
C OS     array of NOS coefficients in an orthogonal series.
C NOS    number of coefficients in OS.
C ETA    requested accuracy of series.
C***REFERENCES  (NONE)
C***ROUTINES CALLED  XERROR
C***END PROLOGUE  INITS
      DIMENSION OS(NOS)
C***FIRST EXECUTABLE STATEMENT  INITS
      IF (NOS.LT.1) CALL XERROR ( 'INITS   NUMBER OF COEFFICIENTS LT 1',
     1 35, 2, 2)
C
      ERR = 0.
      DO 10 II=1,NOS
        I = NOS + 1 - II
        ERR = ERR + ABS(OS(I))
        IF (ERR.GT.ETA) GO TO 20
10    CONTINUE
C
20    IF (I.EQ.NOS) CALL XERROR ( 'INITS   ETA MAY BE TOO SMALL', 28,
     1  1, 2)
      INITS = I
C
      RETURN
      END
      FUNCTION R9LGMC(X)
C***BEGIN PROLOGUE  R9LGMC
C***DATE WRITTEN   770801   (YYMMDD)
C***REVISION DATE  820801   (YYMMDD)
C***CATEGORY NO.  C7E
C***KEYWORDS  CORRECTION FACTOR,LOG GAMMA,SPECIAL FUNCTION
C***AUTHOR  FULLERTON, W., (LANL)
C***PURPOSE  Computes the log Gamma correction factor so that
C            ALOG(GAMMA(X)) = ALOG(SQRT(2*PI)) + (X-.5)*ALOG(X) - X
C            + R9LGMC(X)
C***DESCRIPTION
C
C Compute the log gamma correction factor for X .GE. 10.0 so that
C  ALOG (GAMMA(X)) = ALOG(SQRT(2*PI)) + (X-.5)*ALOG(X) - X + R9LGMC(X)
C
C Series for ALGM       on the interval  0.          to  1.00000D-02
C                                        with weighted error   3.40E-16
C                                         log weighted error  15.47
C                               significant figures required  14.39
C                                    decimal places required  15.86
C***REFERENCES  (NONE)
C***ROUTINES CALLED  CSEVL,INITS,R1MACH,XERROR
C***END PROLOGUE  R9LGMC
      DIMENSION ALGMCS(6)
      DATA ALGMCS( 1) /    .1666389480 45186E0 /
      DATA ALGMCS( 2) /   -.0000138494 817606E0 /
      DATA ALGMCS( 3) /    .0000000098 108256E0 /
      DATA ALGMCS( 4) /   -.0000000000 180912E0 /
      DATA ALGMCS( 5) /    .0000000000 000622E0 /
      DATA ALGMCS( 6) /   -.0000000000 000003E0 /
      DATA NALGM, XBIG, XMAX / 0, 2*0.0 /
C***FIRST EXECUTABLE STATEMENT  R9LGMC
      IF (NALGM.NE.0) GO TO 10
      NALGM = INITS (ALGMCS, 6, R1MACH(3))
      XBIG = 1.0/SQRT(R1MACH(3))
      XMAX = EXP (AMIN1(ALOG(R1MACH(2)/12.0), -ALOG(12.0*R1MACH(1))) )
C
10    IF (X.LT.10.0) CALL XERROR ( 'R9LGMC  X MUST BE GE 10', 23, 1, 2)
      IF (X.GE.XMAX) GO TO 20
C
      R9LGMC = 1.0/(12.0*X)
      IF (X.LT.XBIG) R9LGMC = CSEVL (2.0*(10./X)**2-1., ALGMCS, NALGM)/X
      RETURN
C
20    R9LGMC = 0.0
      CALL XERROR ( 'R9LGMC  X SO BIG R9LGMC UNDERFLOWS', 34, 2, 1)
      RETURN
C
      END
      FUNCTION CSEVL(X,CS,N)
C***BEGIN PROLOGUE  CSEVL
C***DATE WRITTEN   770401   (YYMMDD)
C***REVISION DATE  820801   (YYMMDD)
C***CATEGORY NO.  C3A2
C***KEYWORDS  CHEBYSHEV,FNLIB,SPECIAL FUNCTION
C***AUTHOR  FULLERTON, W., (LANL)
C***PURPOSE  Evaluate the N-term Chebyshev series CS at X.
C***DESCRIPTION
C
C Evaluate the N-term Chebyshev series CS at X.  Adapted from
C R. Broucke, Algorithm 446, C.A.C.M., 16, 254 (1973). Also see Fox
C and Parker, Chebyshev Polynomials in Numerical Analysis, Oxford Press,
C page 56.
C
C       Input Arguments --
C X    value at which the series is to be evaluated.
C CS   array of N terms of a Chebyshev series.  In eval-
C      uating CS, only half the first coefficient is summed.
C N    number of terms in array CS.
C***REFERENCES  (NONE)
C***ROUTINES CALLED  XERROR
C***END PROLOGUE  CSEVL
C
       DIMENSION CS(1)
C***FIRST EXECUTABLE STATEMENT  CSEVL
       IF(N.LT.1) CALL XERROR( 'CSEVL   NUMBER OF TERMS LE 0', 28, 2,2)
       IF(N.GT.1000) CALL XERROR ( 'CSEVL   NUMBER OF TERMS GT 1000',
     1   31,3,2)
       IF (X.LT. -1.0 .OR. X.GT. 1.0) CALL XERROR( 'CSEVL   X OUTSIDE (-
     11,+1)', 25, 1, 1)
C
       B1=0.
       B0=0.
       TWOX=2.*X
       DO 10 I=1,N
       B2=B1
       B1=B0
       NI=N+1-I
       B0=TWOX*B1-B2+CS(NI)
10     CONTINUE
C
       CSEVL = 0.5 * (B0-B2)
C
       RETURN
      END
      REAL FUNCTION R1MACH(I)
C***BEGIN PROLOGUE  R1MACH
C***DATE WRITTEN   790101   (YYMMDD)
C***REVISION DATE  831014   (YYMMDD)
C***CATEGORY NO.  R1
C***KEYWORDS  MACHINE CONSTANTS
C***AUTHOR  FOX, P. A., (BELL LABS)
C           HALL, A. D., (BELL LABS)
C           SCHRYER, N. L., (BELL LABS)
C***PURPOSE  Returns single precision machine dependent constants
C***DESCRIPTION
C
C     R1MACH can be used to obtain machine-dependent parameters
C     for the local machine environment.  It is a function
C     subroutine with one (input) argument, and can be called
C     as follows, for example
C
C          A = R1MACH(I)
C
C     where I=1,...,5.  The (output) value of A above is
C     determined by the (input) value of I.  The results for
C     various values of I are discussed below.
C
C  Single-Precision Machine Constants
C  R1MACH(1) = B**(EMIN-1), the smallest positive magnitude.
C  R1MACH(2) = B**EMAX*(1 - B**(-T)), the largest magnitude.
C  R1MACH(3) = B**(-T), the smallest relative spacing.
C  R1MACH(4) = B**(1-T), the largest relative spacing.
C  R1MACH(5) = LOG10(B)
C***REFERENCES  FOX, P.A., HALL, A.D., SCHRYER, N.L, *FRAMEWORK FOR
C                 A PORTABLE LIBRARY*, ACM TRANSACTIONS ON MATHE-
C                 MATICAL SOFTWARE, VOL. 4, NO. 2, JUNE 1978,
C                 PP. 177-188.
C***ROUTINES CALLED  XERROR
C***END PROLOGUE  R1MACH
C
      INTEGER SMALL(2)
      INTEGER LARGE(2)
      INTEGER RIGHT(2)
      INTEGER DIVER(2)
      INTEGER LOG10(2)
C
      REAL RMACH(5)
C
      EQUIVALENCE (RMACH(1),SMALL(1))
      EQUIVALENCE (RMACH(2),LARGE(1))
      EQUIVALENCE (RMACH(3),RIGHT(1))
      EQUIVALENCE (RMACH(4),DIVER(1))
      EQUIVALENCE (RMACH(5),LOG10(1))
C
C     MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM.
C
C     DATA RMACH(1) / Z400800000 /
C     DATA RMACH(2) / Z5FFFFFFFF /
C     DATA RMACH(3) / Z4E9800000 /
C     DATA RMACH(4) / Z4EA800000 /
C     DATA RMACH(5) / Z500E730E8 /
C
C     MACHINE CONSTANTS FOR THE BURROUGHS 5700/6700/7700 SYSTEMS.
C
C     DATA RMACH(1) / O1771000000000000 /
C     DATA RMACH(2) / O0777777777777777 /
C     DATA RMACH(3) / O1311000000000000 /
C     DATA RMACH(4) / O1301000000000000 /
C     DATA RMACH(5) / O1157163034761675 /
C
C     MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES.
C
C     DATA RMACH(1) / 00564000000000000000B /
C     DATA RMACH(2) / 37767777777777777776B /
C     DATA RMACH(3) / 16414000000000000000B /
C     DATA RMACH(4) / 16424000000000000000B /
C     DATA RMACH(5) / 17164642023241175720B /
C
C     MACHINE CONSTANTS FOR THE CRAY 1
C
C     DATA RMACH(1) / 200034000000000000000B /
C     DATA RMACH(2) / 577767777777777777776B /
C     DATA RMACH(3) / 377224000000000000000B /
C     DATA RMACH(4) / 377234000000000000000B /
C     DATA RMACH(5) / 377774642023241175720B /
C
C     MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200
C
C     NOTE - IT MAY BE APPROPRIATE TO INCLUDE THE FOLLOWING CARD -
C     STATIC RMACH(5)
C
C     DATA SMALL/20K,0/,LARGE/77777K,177777K/
C     DATA RIGHT/35420K,0/,DIVER/36020K,0/
C     DATA LOG10/40423K,42023K/
C
C     MACHINE CONSTANTS FOR THE HARRIS 220
C
C     DATA SMALL(1),SMALL(2) / '20000000, '00000201 /
C     DATA LARGE(1),LARGE(2) / '37777777, '00000177 /
C     DATA RIGHT(1),RIGHT(2) / '20000000, '00000352 /
C     DATA DIVER(1),DIVER(2) / '20000000, '00000353 /
C     DATA LOG10(1),LOG10(2) / '23210115, '00000377 /
C
C     MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 SERIES.
C
C     DATA RMACH(1) / O402400000000 /
C     DATA RMACH(2) / O376777777777 /
C     DATA RMACH(3) / O714400000000 /
C     DATA RMACH(4) / O716400000000 /
C     DATA RMACH(5) / O776464202324 /
C
C     MACHINE CONSTANTS FOR THE HP 2100
C
C     3 WORD DOUBLE PRECISION WITH FTN4
C
C     DATA SMALL(1), SMALL(2) / 40000B,       1 /
C     DATA LARGE(1), LARGE(2) / 77777B, 177776B /
C     DATA RIGHT(1), RIGHT(2) / 40000B,    325B /
C     DATA DIVER(1), DIVER(2) / 40000B,    327B /
C     DATA LOG10(1), LOG10(2) / 46420B,  46777B /
C
C     MACHINE CONSTANTS FOR THE HP 2100
C     4 WORD DOUBLE PRECISION WITH FTN4
C
C     DATA SMALL(1), SMALL(2) / 40000B,       1 /
C     DATA LARGE91), LARGE(2) / 77777B, 177776B /
C     DATA RIGHT(1), RIGHT(2) / 40000B,    325B /
C     DATA DIVER(1), DIVER(2) / 40000B,    327B /
C     DATA LOG10(1), LOG10(2) / 46420B,  46777B /
C
C     MACHINE CONSTANTS FOR THE IBM 360/370 SERIES,
C     THE XEROX SIGMA 5/7/9, THE SEL SYSTEMS 85/86  AND
C     THE PERKIN ELMER (INTERDATA) 7/32.
C
C     DATA RMACH(1) / Z00100000 /
C     DATA RMACH(2) / Z7FFFFFFF /
C     DATA RMACH(3) / Z3B100000 /
C     DATA RMACH(4) / Z3C100000 /
C     DATA RMACH(5) / Z41134413 /
C
C     MACHINE CONSTANTS FOR THE PDP-10 (KA OR KI PROCESSOR).
C
C     DATA RMACH(1) / "000400000000 /
C     DATA RMACH(2) / "377777777777 /
C     DATA RMACH(3) / "146400000000 /
C     DATA RMACH(4) / "147400000000 /
C     DATA RMACH(5) / "177464202324 /
C
C     MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING
C     32-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL).
C
C     DATA SMALL(1) /    8388608 /
C     DATA LARGE(1) / 2147483647 /
C     DATA RIGHT(1) /  880803840 /
C     DATA DIVER(1) /  889192448 /
C     DATA LOG10(1) / 1067065499 /
C
C     DATA RMACH(1) / O00040000000 /
C     DATA RMACH(2) / O17777777777 /
C     DATA RMACH(3) / O06440000000 /
C     DATA RMACH(4) / O06500000000 /
C     DATA RMACH(5) / O07746420233 /
C
C     MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING
C     16-BIT INTEGERS  (EXPRESSED IN INTEGER AND OCTAL).
C
C     DATA SMALL(1),SMALL(2) /   128,     0 /
C     DATA LARGE(1),LARGE(2) / 32767,    -1 /
C     DATA RIGHT(1),RIGHT(2) / 13440,     0 /
C     DATA DIVER(1),DIVER(2) / 13568,     0 /
C     DATA LOG10(1),LOG10(2) / 16282,  8347 /
C
C     DATA SMALL(1),SMALL(2) / O000200, O000000 /
C     DATA LARGE(1),LARGE(2) / O077777, O177777 /
C     DATA RIGHT(1),RIGHT(2) / O032200, O000000 /
C     DATA DIVER(1),DIVER(2) / O032400, O000000 /
C     DATA LOG10(1),LOG10(2) / O037632, O020233 /
C
C     MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES.
C
C     DATA RMACH(1) / O000400000000 /
C     DATA RMACH(2) / O377777777777 /
C     DATA RMACH(3) / O146400000000 /
C     DATA RMACH(4) / O147400000000 /
C     DATA RMACH(5) / O177464202324 /
C
C     MACHINE CONSTANTS FOR THE VAX 11/780
C    (EXPRESSED IN INTEGER AND HEXADECIMAL)
C  ***THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSTEMS***
C  *** THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS***
C
      DATA SMALL(1) /       128 /
      DATA LARGE(1) /    -32769 /
      DATA RIGHT(1) /     13440 /
      DATA DIVER(1) /     13568 /
      DATA LOG10(1) / 547045274 /
C
C     DATA SMALL(1) / Z00000080 /
C     DATA LARGE(1) / ZFFFF7FFF /
C     DATA RIGHT(1) / Z00003480 /
C     DATA DIVER(1) / Z00003500 /
C     DATA LOG10(1) / Z209B3F9A /
C
C     MACHINE CONSTANTS FOR THE Z80 MICROPROCESSOR
C
C     DATA SMALL(1),SMALL(2) /     0,    256/
C     DATA LARGE(1),LARGE(2) /    -1,   -129/
C     DATA RIGHT(1),RIGHT(2) /     0,  26880/
C     DATA DIVER(1),DIVER(2) /     0,  27136/
C     DATA LOG10(1),LOG10(2) /  8347,  32538/
C
C
C***FIRST EXECUTABLE STATEMENT  R1MACH
      IF (I .LT. 1  .OR.  I .GT. 5)
     1   CALL XERROR ( 'R1MACH -- I OUT OF BOUNDS',25,1,2)
C
      R1MACH = RMACH(I)
      RETURN
C
      END
      SUBROUTINE XERROR(MESSG,NMESSG,NERR,LEVEL)
C***BEGIN PROLOGUE  XERROR
C***DATE WRITTEN   790801   (YYMMDD)
C***REVISION DATE  820801   (YYMMDD)
C***CATEGORY NO.  R3C
C***KEYWORDS  ERROR,XERROR PACKAGE
C***AUTHOR  JONES, R. E., (SNLA)
C***PURPOSE  Processes an error (diagnostic) message
C***DESCRIPTION
C
C     Abstract
C        XERROR processes a diagnostic message, in a manner
C        determined by the value of LEVEL and the current value
C        of the library error control flag, KONTRL.
C        (See subroutine XSETF for details.)
C
C     Description of Parameters
C      --Input--
C        MESSG - the Hollerith message to be processed, containing
C                no more than 72 characters.
C        NMESSG- the actual number of characters in MESSG.
C        NERR  - the error number associated with this message.
C                NERR must not be zero.
C        LEVEL - error category.
C                =2 means this is an unconditionally fatal error.
C                =1 means this is a recoverable error.  (I.e., it is
C                   non-fatal if XSETF has been appropriately called.)
C                =0 means this is a warning message only.
C                =-1 means this is a warning message which is to be
C                   printed at most once, regardless of how many
C                   times this call is executed.
C
C     Examples
C        CALL XERROR(23HSMOOTH -- NUM WAS ZERO.,23,1,2)
C        CALL XERROR(43HINTEG  -- LESS THAN FULL ACCURACY ACHIEVED.,
C                    43,2,1)
C        CALL XERROR(65HROOTER -- ACTUAL ZERO OF F FOUND BEFORE INTERVAL
C    1 FULLY COLLAPSED.,65,3,0)
C        CALL XERROR(39HEXP    -- UNDERFLOWS BEING SET TO ZERO.,39,1,-1)
C
C     Written by Ron Jones, with SLATEC Common Math Library Subcommittee
C     Latest revision ---  7 Feb 1979
C***REFERENCES  JONES R.E., *SLATEC COMMON MATHEMATICAL LIBRARY ERROR
C                 HANDLING PACKAGE*, SAND78-1189, SANDIA LABORATORIES,
C                 1978.
C***ROUTINES CALLED  XERRWV
C***END PROLOGUE  XERROR
C
C     DIMENSION MESSG(NMESSG)
C***FIRST EXECUTABLE STATEMENT  XERROR
C     CALL XERRWV(MESSG,NMESSG,NERR,LEVEL,0,0,0,0,0.,0.)
C     RETURN
C     END
C *
C *  MODIFIED S.A.HOWARD, XERRWV WAS NOT OBTAINED FROM LANL THEREFORE
C *  THIS SUBROUTINE WAS KLUDGED.  THIS VERSION WILL PRINT THE MESSAGE
C *  AND ACT BASED ON THE FLAG *LEVEL*.  IF LEVEL=2 EXECUTION IS HALTED.
C *
      CHARACTER MESSG(NMESSG)
      WRITE(6,1) LEVEL,NERR,(MESSG(I),I=1,NMESSG)
1     FORMAT(' -GAMMA FUNCTION ROUTINES: ERROR LEVEL ',I2,'  NUMBER ',
     1 I2,/,1X,72A1)
      IF(LEVEL.EQ.2) STOP'--GAMMA FUNCTION ROUTINES: ERROR LEVEL 2'
      RETURN
      END
      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
      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 LANL DEPENDENT CODE REMOVED 81.02.04
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

