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
r9lgit.f
Go to the documentation of this file.
1 *DECK R9LGIT
2  FUNCTION r9lgit (A, X, ALGAP1)
3 C***BEGIN PROLOGUE R9LGIT
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 SINGLE 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 R1MACH, 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 R9LGIT
29  SAVE eps, sqeps
30  DATA eps, sqeps / 2*0.0 /
31 C***FIRST EXECUTABLE STATEMENT R9LGIT
32  IF (eps.EQ.0.0) eps = 0.5*r1mach(3)
33  IF (sqeps.EQ.0.0) sqeps = sqrt(r1mach(4))
34 C
35  IF (x .LE. 0.0 .OR. a .LT. x) CALL xermsg('SLATEC', 'R9LGIT',
36  + 'X SHOULD BE GT 0.0 AND LE A', 2, 2)
37 C
38  ax = a + x
39  a1x = ax + 1.0
40  r = 0.0
41  p = 1.0
42  s = p
43  DO 20 k=1,200
44  fk = k
45  t = (a+fk)*x*(1.0+r)
46  r = t/((ax+fk)*(a1x+fk)-t)
47  p = r*p
48  s = s + p
49  IF (abs(p).LT.eps*s) go to 30
50  20 CONTINUE
51  CALL xermsg('SLATEC', 'R9LGIT',
52  + 'NO CONVERGENCE IN 200 TERMS OF CONTINUED FRACTION', 3, 2)
53 C
54  30 hstar = 1.0 - x*s/a1x
55  IF (hstar .LT. sqeps) CALL xermsg('SLATEC', 'R9LGIT',
56  + 'RESULT LESS THAN HALF PRECISION', 1, 1)
57 C
58  r9lgit = -x - algap1 - log(hstar)
59 C
60  RETURN
61  END