Esope source files

   1 : C GAMMA     SOURCE    KK2000    13/11/08    21:15:41     7860           
   2 :       SUBROUTINE GAMMA(X,GA)
   3 :       IMPLICIT INTEGER(I-N)
   4 :       IMPLICIT REAL*8 (A-H,O-Z)
   5 : 
   6 : -INC CCREEL
   7 : C
   8 : C     ==================================================
   9 : C     Programme de : http://jin.ece.uiuc.edu/routines/routines.html
  10 : C     Purpose: Compute the gamma function G(x)
  11 : C     Input :  x  --- Argument of G(x)
  12 : C                     (x is not equal to 0,-1,-2,...)
  13 : C     Output:  GA --- G(x)
  14 : C     ==================================================
  15 : C
  16 :       DIMENSION G(26)
  17 :       DATA G/ 1D0                 , .5772156649015329D0,
  18 :      &       -.6558780715202538D0 ,-.420026350340952D-1,
  19 :      &        .1665386113822915D0 ,-.421977345555443D-1,
  20 :      &       -.96219715278770D-2  , .72189432466630D-2,
  21 :      &       -.11651675918591D-2  ,-.2152416741149D-3,
  22 :      &        .1280502823882D-3   ,-.201348547807D-4,
  23 :      &       -.12504934821D-5     , .11330272320D-5,
  24 :      &       -.2056338417D-6      , .61160950D-8,
  25 :      &        .50020075D-8        ,-.11812746D-8,
  26 :      &        .1043427D-9         , .77823D-11,
  27 :      &       -.36968D-11          , .51D-12,
  28 :      &       -.206D-13            ,-.54D-14,
  29 :      &        .14D-14             , .1D-15 /
  30 : 
  31 : C
  32 :       XINT=FLOAT(INT(X))
  33 :       IF (X.EQ.XINT) THEN
  34 :          IF (X.GT.0D0) THEN
  35 :             GA=1D0
  36 :             M1=INT(X-1.D0)
  37 :             DO 10 K=2,M1
  38 : 10             GA=GA*FLOAT(K)
  39 :          ELSE
  40 :             GA=1D0+300D0
  41 :          ENDIF
  42 :       ELSE
  43 :          R=1D0
  44 :          IF (ABS(X).GT.1.D0) THEN
  45 :             Z=ABS(X)
  46 :             M=INT(Z)
  47 :             DO 15 K=1,M
  48 : 15             R=R*(Z-FLOAT(K))
  49 :             Z=Z-FLOAT(M)
  50 :          ELSE
  51 :             Z=X
  52 :          ENDIF
  53 :          GR=G(26)
  54 :          DO 20 K=25,1,-1
  55 : 20          GR=GR*Z+G(K)
  56 :          GA=1D0/(GR*Z)
  57 :          IF (ABS(X).GT.1.D0) THEN
  58 :             GA=GA*R
  59 :             IF (X.LT.0.D0) GA=-XPI/(X*GA*SIN(XPI*X))
  60 :          ENDIF
  61 :       ENDIF
  62 :       RETURN
  63 :       END
  64 :  
  65 :