| 1 | SUBROUTINE GAMMAFUNC(X,GA) |
|---|
| 2 | C |
|---|
| 3 | C ================================================== |
|---|
| 4 | C Purpose: Compute the gamma function â(x) |
|---|
| 5 | C Input : x --- Argument of â(x) |
|---|
| 6 | C ( x is not equal to 0,-1,-2,úúú ) |
|---|
| 7 | C Output: GA --- â(x) |
|---|
| 8 | C ================================================== |
|---|
| 9 | C |
|---|
| 10 | IMPLICIT DOUBLE PRECISION (A-H,O-Z) |
|---|
| 11 | DIMENSION G(26) |
|---|
| 12 | PI=3.141592653589793D0 |
|---|
| 13 | IF (X.EQ.INT(X)) THEN |
|---|
| 14 | IF (X.GT.0.0D0) THEN |
|---|
| 15 | GA=1.0D0 |
|---|
| 16 | M1=X-1 |
|---|
| 17 | DO 10 K=2,M1 |
|---|
| 18 | 10 GA=GA*K |
|---|
| 19 | ELSE |
|---|
| 20 | GA=1.0D+300 |
|---|
| 21 | ENDIF |
|---|
| 22 | ELSE |
|---|
| 23 | IF (DABS(X).GT.1.0D0) THEN |
|---|
| 24 | Z=DABS(X) |
|---|
| 25 | M=INT(Z) |
|---|
| 26 | R=1.0D0 |
|---|
| 27 | DO 15 K=1,M |
|---|
| 28 | 15 R=R*(Z-K) |
|---|
| 29 | Z=Z-M |
|---|
| 30 | ELSE |
|---|
| 31 | Z=X |
|---|
| 32 | ENDIF |
|---|
| 33 | DATA G/1.0D0,0.5772156649015329D0, |
|---|
| 34 | & -0.6558780715202538D0, -0.420026350340952D-1, |
|---|
| 35 | & 0.1665386113822915D0,-.421977345555443D-1, |
|---|
| 36 | & -.96219715278770D-2, .72189432466630D-2, |
|---|
| 37 | & -.11651675918591D-2, -.2152416741149D-3, |
|---|
| 38 | & .1280502823882D-3, -.201348547807D-4, |
|---|
| 39 | & -.12504934821D-5, .11330272320D-5, |
|---|
| 40 | & -.2056338417D-6, .61160950D-8, |
|---|
| 41 | & .50020075D-8, -.11812746D-8, |
|---|
| 42 | & .1043427D-9, .77823D-11, |
|---|
| 43 | & -.36968D-11, .51D-12, |
|---|
| 44 | & -.206D-13, -.54D-14, .14D-14, .1D-15/ |
|---|
| 45 | GR=G(26) |
|---|
| 46 | DO 20 K=25,1,-1 |
|---|
| 47 | 20 GR=GR*Z+G(K) |
|---|
| 48 | GA=1.0D0/(GR*Z) |
|---|
| 49 | IF (DABS(X).GT.1.0D0) THEN |
|---|
| 50 | GA=GA*R |
|---|
| 51 | IF (X.LT.0.0D0) GA=-PI/(X*GA*DSIN(PI*X)) |
|---|
| 52 | ENDIF |
|---|
| 53 | ENDIF |
|---|
| 54 | RETURN |
|---|
| 55 | END |
|---|
