GNU Octave  4.4.1
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
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
OCTAVE_EXPORT octave_value_list or N dimensional array whose elements are all equal to the base of natural logarithms The constant ex $e satisfies the equation log(e)
static T abs(T x)
Definition: pr-output.cc:1696
#define eps(C)
real function r1mach(i)
Definition: r1mach.f:20