GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
dgamit.f
Go to the documentation of this file.
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