|
3111
|
1 *DECK DGAMIT |
|
|
2 DOUBLE PRECISION FUNCTION DGAMIT (A, X) |
|
|
3 C***BEGIN PROLOGUE DGAMIT |
|
|
4 C***PURPOSE Calculate Tricomi's form of the incomplete Gamma function. |
|
|
5 C***LIBRARY SLATEC (FNLIB) |
|
|
6 C***CATEGORY C7E |
|
|
7 C***TYPE DOUBLE PRECISION (GAMIT-S, DGAMIT-D) |
|
|
8 C***KEYWORDS COMPLEMENTARY INCOMPLETE GAMMA FUNCTION, FNLIB, |
|
|
9 C SPECIAL FUNCTIONS, TRICOMI |
|
|
10 C***AUTHOR Fullerton, W., (LANL) |
|
|
11 C***DESCRIPTION |
|
|
12 C |
|
|
13 C Evaluate Tricomi's incomplete Gamma function defined by |
|
|
14 C |
|
|
15 C DGAMIT = X**(-A)/GAMMA(A) * integral from 0 to X of EXP(-T) * |
|
|
16 C T**(A-1.) |
|
|
17 C |
|
|
18 C for A .GT. 0.0 and by analytic continuation for A .LE. 0.0. |
|
|
19 C GAMMA(X) is the complete gamma function of X. |
|
|
20 C |
|
|
21 C DGAMIT is evaluated for arbitrary real values of A and for non- |
|
|
22 C negative values of X (even though DGAMIT is defined for X .LT. |
|
|
23 C 0.0), except that for X = 0 and A .LE. 0.0, DGAMIT is infinite, |
|
|
24 C which is a fatal error. |
|
|
25 C |
|
|
26 C The function and both arguments are DOUBLE PRECISION. |
|
|
27 C |
|
|
28 C A slight deterioration of 2 or 3 digits accuracy will occur when |
|
|
29 C DGAMIT is very large or very small in absolute value, because log- |
|
|
30 C arithmic variables are used. Also, if the parameter A is very |
|
|
31 C close to a negative integer (but not a negative integer), there is |
|
|
32 C a loss of accuracy, which is reported if the result is less than |
|
|
33 C half machine precision. |
|
|
34 C |
|
|
35 C***REFERENCES W. Gautschi, A computational procedure for incomplete |
|
|
36 C gamma functions, ACM Transactions on Mathematical |
|
|
37 C Software 5, 4 (December 1979), pp. 466-481. |
|
|
38 C W. Gautschi, Incomplete gamma functions, Algorithm 542, |
|
|
39 C ACM Transactions on Mathematical Software 5, 4 |
|
|
40 C (December 1979), pp. 482-489. |
|
|
41 C***ROUTINES CALLED D1MACH, D9GMIT, D9LGIC, D9LGIT, DGAMR, DLGAMS, |
|
|
42 C DLNGAM, XERCLR, XERMSG |
|
|
43 C***REVISION HISTORY (YYMMDD) |
|
|
44 C 770701 DATE WRITTEN |
|
|
45 C 890531 Changed all specific intrinsics to generic. (WRB) |
|
|
46 C 890531 REVISION DATE from Version 3.2 |
|
|
47 C 891214 Prologue converted to Version 4.0 format. (BAB) |
|
|
48 C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) |
|
|
49 C 920528 DESCRIPTION and REFERENCES sections revised. (WRB) |
|
|
50 C***END PROLOGUE DGAMIT |
|
|
51 DOUBLE PRECISION A, X, AEPS, AINTA, ALGAP1, ALNEPS, ALNG, ALX, |
|
|
52 1 BOT, H, SGA, SGNGAM, SQEPS, T, D1MACH, DGAMR, D9GMIT, D9LGIT, |
|
|
53 2 DLNGAM, D9LGIC |
|
|
54 LOGICAL FIRST |
|
|
55 SAVE ALNEPS, SQEPS, BOT, FIRST |
|
|
56 DATA FIRST /.TRUE./ |
|
|
57 C***FIRST EXECUTABLE STATEMENT DGAMIT |
|
|
58 IF (FIRST) THEN |
|
|
59 ALNEPS = -LOG (D1MACH(3)) |
|
|
60 SQEPS = SQRT(D1MACH(4)) |
|
|
61 BOT = LOG (D1MACH(1)) |
|
|
62 ENDIF |
|
|
63 FIRST = .FALSE. |
|
|
64 C |
|
|
65 IF (X .LT. 0.D0) CALL XERMSG ('SLATEC', 'DGAMIT', 'X IS NEGATIVE' |
|
|
66 + , 2, 2) |
|
|
67 C |
|
|
68 IF (X.NE.0.D0) ALX = LOG (X) |
|
|
69 SGA = 1.0D0 |
|
|
70 IF (A.NE.0.D0) SGA = SIGN (1.0D0, A) |
|
|
71 AINTA = AINT (A + 0.5D0*SGA) |
|
|
72 AEPS = A - AINTA |
|
|
73 C |
|
|
74 IF (X.GT.0.D0) GO TO 20 |
|
|
75 DGAMIT = 0.0D0 |
|
|
76 IF (AINTA.GT.0.D0 .OR. AEPS.NE.0.D0) DGAMIT = DGAMR(A+1.0D0) |
|
|
77 RETURN |
|
|
78 C |
|
|
79 20 IF (X.GT.1.D0) GO TO 30 |
|
|
80 IF (A.GE.(-0.5D0) .OR. AEPS.NE.0.D0) CALL DLGAMS (A+1.0D0, ALGAP1, |
|
|
81 1 SGNGAM) |
|
|
82 DGAMIT = D9GMIT (A, X, ALGAP1, SGNGAM, ALX) |
|
|
83 RETURN |
|
|
84 C |
|
|
85 30 IF (A.LT.X) GO TO 40 |
|
|
86 T = D9LGIT (A, X, DLNGAM(A+1.0D0)) |
|
|
87 IF (T.LT.BOT) CALL XERCLR |
|
|
88 DGAMIT = EXP (T) |
|
|
89 RETURN |
|
|
90 C |
|
|
91 40 ALNG = D9LGIC (A, X, ALX) |
|
|
92 C |
|
|
93 C EVALUATE DGAMIT IN TERMS OF LOG (DGAMIC (A, X)) |
|
|
94 C |
|
|
95 H = 1.0D0 |
|
|
96 IF (AEPS.EQ.0.D0 .AND. AINTA.LE.0.D0) GO TO 50 |
|
|
97 C |
|
|
98 CALL DLGAMS (A+1.0D0, ALGAP1, SGNGAM) |
|
|
99 T = LOG (ABS(A)) + ALNG - ALGAP1 |
|
|
100 IF (T.GT.ALNEPS) GO TO 60 |
|
|
101 C |
|
|
102 IF (T.GT.(-ALNEPS)) H = 1.0D0 - SGA * SGNGAM * EXP(T) |
|
|
103 IF (ABS(H).GT.SQEPS) GO TO 50 |
|
|
104 C |
|
|
105 CALL XERCLR |
|
|
106 CALL XERMSG ('SLATEC', 'DGAMIT', 'RESULT LT HALF PRECISION', 1, |
|
|
107 + 1) |
|
|
108 C |
|
|
109 50 T = -A*ALX + LOG(ABS(H)) |
|
|
110 IF (T.LT.BOT) CALL XERCLR |
|
|
111 DGAMIT = SIGN (EXP(T), H) |
|
|
112 RETURN |
|
|
113 C |
|
|
114 60 T = T - A*ALX |
|
|
115 IF (T.LT.BOT) CALL XERCLR |
|
|
116 DGAMIT = -SGA * SGNGAM * EXP(T) |
|
|
117 RETURN |
|
|
118 C |
|
|
119 END |
