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
d9lgit.f
Go to the documentation of this file.
1 *DECK D9LGIT
2  DOUBLE PRECISION FUNCTION d9lgit (A, X, ALGAP1)
3 C***BEGIN PROLOGUE D9LGIT
4 C***SUBSIDIARY
5 C***PURPOSE Compute the logarithm of Tricomi's incomplete Gamma
6 C function with Perron's continued fraction for large X and
7 C A .GE. X.
8 C***LIBRARY SLATEC (FNLIB)
9 C***CATEGORY C7E
10 C***TYPE DOUBLE PRECISION (R9LGIT-S, D9LGIT-D)
11 C***KEYWORDS FNLIB, INCOMPLETE GAMMA FUNCTION, LOGARITHM,
12 C PERRON'S CONTINUED FRACTION, SPECIAL FUNCTIONS, TRICOMI
13 C***AUTHOR Fullerton, W., (LANL)
14 C***DESCRIPTION
15 C
16 C Compute the log of Tricomi's incomplete gamma function with Perron's
17 C continued fraction for large X and for A .GE. X.
18 C
19 C***REFERENCES (NONE)
20 C***ROUTINES CALLED D1MACH, XERMSG
21 C***REVISION HISTORY (YYMMDD)
22 C 770701 DATE WRITTEN
23 C 890531 Changed all specific intrinsics to generic. (WRB)
24 C 890531 REVISION DATE from Version 3.2
25 C 891214 Prologue converted to Version 4.0 format. (BAB)
26 C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
27 C 900720 Routine changed from user-callable to subsidiary. (WRB)
28 C***END PROLOGUE D9LGIT
29  DOUBLE PRECISION a, x, algap1, ax, a1x, eps, fk, hstar, p, r, s,
30  1 sqeps, t, d1mach
31  LOGICAL first
32  SAVE eps, sqeps, first
33  DATA first /.true./
34 C***FIRST EXECUTABLE STATEMENT D9LGIT
35  IF (first) THEN
36  eps = 0.5d0*d1mach(3)
37  sqeps = sqrt(d1mach(4))
38  ENDIF
39  first = .false.
40 C
41  IF (x .LE. 0.d0 .OR. a .LT. x) CALL xermsg('SLATEC', 'D9LGIT',
42  + 'X SHOULD BE GT 0.0 AND LE A', 2, 2)
43 C
44  ax = a + x
45  a1x = ax + 1.0d0
46  r = 0.d0
47  p = 1.d0
48  s = p
49  DO 20 k=1,200
50  fk = k
51  t = (a+fk)*x*(1.d0+r)
52  r = t/((ax+fk)*(a1x+fk)-t)
53  p = r*p
54  s = s + p
55  IF (abs(p).LT.eps*s) go to 30
56  20 CONTINUE
57  CALL xermsg('SLATEC', 'D9LGIT',
58  + 'NO CONVERGENCE IN 200 TERMS OF CONTINUED FRACTION', 3, 2)
59 C
60  30 hstar = 1.0d0 - x*s/a1x
61  IF (hstar .LT. sqeps) CALL xermsg('SLATEC', 'D9LGIT',
62  + 'RESULT LESS THAN HALF PRECISION', 1, 1)
63 C
64  d9lgit = -x - algap1 - log(hstar)
65  RETURN
66 C
67  END