gamma.f

Go to the documentation of this file.

*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