GNU Octave  6.2.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
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
#define eps(C)
double precision function d1mach(i)
Definition: d1mach.f:23
double precision function d9lgit(A, X, ALGAP1)
Definition: d9lgit.f:3
static T abs(T x)
Definition: pr-output.cc:1678
subroutine xermsg(LIBRAR, SUBROU, MESSG, NERR, LEVEL)
Definition: xermsg.f:3