*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