relibc_openlibm/slatec/gamma.f

*DECK GAMMA
      FUNCTION GAMMA (X)
C***BEGIN PROLOGUE  GAMMA
C***PURPOSE  Compute the complete Gamma function.
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C7A
C***TYPE      SINGLE PRECISION (GAMMA-S, DGAMMA-D, CGAMMA-C)
C***KEYWORDS  COMPLETE GAMMA FUNCTION, FNLIB, SPECIAL FUNCTIONS
C***AUTHOR  Fullerton, W., (LANL)
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
C***REFERENCES  (NONE)
C***ROUTINES CALLED  CSEVL, GAMLIM, INITS, R1MACH, R9LGMC, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   770601  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890531  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C***END PROLOGUE  GAMMA
      DIMENSION GCS(23)
      LOGICAL FIRST
      SAVE GCS, PI, SQ2PIL, NGCS, XMIN, XMAX, DXREL, FIRST
      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 LOG (SQRT (2.*PI) )
      DATA SQ2PIL /0.91893 85332 04672 74E0/
      DATA FIRST /.TRUE./
C
C LANL DEPENDENT CODE REMOVED 81.02.04
C
C***FIRST EXECUTABLE STATEMENT  GAMMA
      IF (FIRST) THEN
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
      ENDIF
      FIRST = .FALSE.
C
      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 - 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 XERMSG ('SLATEC', 'GAMMA', 'X IS 0', 4, 2)
      IF (X .LT. 0. .AND. X+N-2 .EQ. 0.) CALL XERMSG ('SLATEC', 'GAMMA'
     1, 'X IS A NEGATIVE INTEGER', 4, 2)
      IF (X .LT. (-0.5) .AND. ABS((X-AINT(X-0.5))/X) .LT. DXREL) CALL
     1XERMSG ( 'SLATEC', 'GAMMA',
     2'ANSWER LT HALF PRECISION BECAUSE X TOO NEAR NEGATIVE INTEGER'
     3, 1, 1)
C
      DO 20 I=1,N
        GAMMA = GAMMA / (X+I-1)
 20   CONTINUE
      RETURN
C
C GAMMA(X) FOR X .GE. 2.
C
 30   DO 40 I=1,N
        GAMMA = (Y+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 XERMSG ('SLATEC', 'GAMMA',
     +   'X SO BIG GAMMA OVERFLOWS', 3, 2)
C
      GAMMA = 0.
      IF (X .LT. XMIN) CALL XERMSG ('SLATEC', 'GAMMA',
     +   'X SO SMALL GAMMA UNDERFLOWS', 2, 1)
      IF (X.LT.XMIN) RETURN
C
      GAMMA = EXP((Y-0.5)*LOG(Y) - Y + SQ2PIL + R9LGMC(Y) )
      IF (X.GT.0.) RETURN
C
      IF (ABS((X-AINT(X-0.5))/X) .LT. DXREL) CALL XERMSG ('SLATEC',
     +   'GAMMA',
     +   'ANSWER LT HALF PRECISION, X TOO NEAR NEGATIVE INTEGER', 1, 1)
C
      SINPIY = SIN (PI*Y)
      IF (SINPIY .EQ. 0.) CALL XERMSG ('SLATEC', 'GAMMA',
     +   'X IS A NEGATIVE INTEGER', 4, 2)
C
      GAMMA = -PI / (Y*SINPIY*GAMMA)
C
      RETURN
      END