C GAMMA SOURCE KK2000 13/11/08 21:15:41 7860 SUBROUTINE GAMMA(X,GA) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) -INC CCREEL C C ================================================== C Programme de : http://jin.ece.uiuc.edu/routines/routines.html C Purpose: Compute the gamma function G(x) C Input : x --- Argument of G(x) C (x is not equal to 0,-1,-2,...) C Output: GA --- G(x) C ================================================== C DIMENSION G(26) DATA G/ 1D0 , .5772156649015329D0, & -.6558780715202538D0 ,-.420026350340952D-1, & .1665386113822915D0 ,-.421977345555443D-1, & -.96219715278770D-2 , .72189432466630D-2, & -.11651675918591D-2 ,-.2152416741149D-3, & .1280502823882D-3 ,-.201348547807D-4, & -.12504934821D-5 , .11330272320D-5, & -.2056338417D-6 , .61160950D-8, & .50020075D-8 ,-.11812746D-8, & .1043427D-9 , .77823D-11, & -.36968D-11 , .51D-12, & -.206D-13 ,-.54D-14, & .14D-14 , .1D-15 / C XINT=FLOAT(INT(X)) IF (X.EQ.XINT) THEN IF (X.GT.0D0) THEN GA=1D0 M1=INT(X-1.D0) DO 10 K=2,M1 10 GA=GA*FLOAT(K) ELSE GA=1D0+300D0 ENDIF ELSE R=1D0 IF (ABS(X).GT.1.D0) THEN Z=ABS(X) M=INT(Z) DO 15 K=1,M 15 R=R*(Z-FLOAT(K)) Z=Z-FLOAT(M) ELSE Z=X ENDIF GR=G(26) DO 20 K=25,1,-1 20 GR=GR*Z+G(K) GA=1D0/(GR*Z) IF (ABS(X).GT.1.D0) THEN GA=GA*R IF (X.LT.0.D0) GA=-XPI/(X*GA*SIN(XPI*X)) ENDIF ENDIF RETURN END