GNU Octave  4.0.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
d9lgic.f
Go to the documentation of this file.
1 *DECK D9LGIC
2  DOUBLE PRECISION FUNCTION d9lgic (A, X, ALX)
3 C***BEGIN PROLOGUE D9LGIC
4 C***SUBSIDIARY
5 C***PURPOSE Compute the log complementary incomplete Gamma function
6 C for large X and for A .LE. X.
7 C***LIBRARY SLATEC (FNLIB)
8 C***CATEGORY C7E
9 C***TYPE DOUBLE PRECISION (R9LGIC-S, D9LGIC-D)
10 C***KEYWORDS COMPLEMENTARY INCOMPLETE GAMMA FUNCTION, FNLIB, LARGE X,
11 C LOGARITHM, SPECIAL FUNCTIONS
12 C***AUTHOR Fullerton, W., (LANL)
13 C***DESCRIPTION
14 C
15 C Compute the log complementary incomplete gamma function for large X
16 C and for A .LE. X.
17 C
18 C***REFERENCES (NONE)
19 C***ROUTINES CALLED D1MACH, XERMSG
20 C***REVISION HISTORY (YYMMDD)
21 C 770701 DATE WRITTEN
22 C 890531 Changed all specific intrinsics to generic. (WRB)
23 C 890531 REVISION DATE from Version 3.2
24 C 891214 Prologue converted to Version 4.0 format. (BAB)
25 C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
26 C 900720 Routine changed from user-callable to subsidiary. (WRB)
27 C***END PROLOGUE D9LGIC
28  DOUBLE PRECISION A, X, ALX, EPS, FK, P, R, S, T, XMA, XPA, D1MACH
29  SAVE eps
30  DATA eps / 0.d0 /
31 C***FIRST EXECUTABLE STATEMENT D9LGIC
32  IF (eps.EQ.0.d0) eps = 0.5d0*d1mach(3)
33 C
34  xpa = x + 1.0d0 - a
35  xma = x - 1.d0 - a
36 C
37  r = 0.d0
38  p = 1.d0
39  s = p
40  DO 10 k=1,300
41  fk = k
42  t = fk*(a-fk)*(1.d0+r)
43  r = -t/((xma+2.d0*fk)*(xpa+2.d0*fk)+t)
44  p = r*p
45  s = s + p
46  IF (abs(p).LT.eps*s) go to 20
47  10 CONTINUE
48  CALL xermsg('SLATEC', 'D9LGIC',
49  + 'NO CONVERGENCE IN 300 TERMS OF CONTINUED FRACTION', 1, 2)
50 C
51  20 d9lgic = a*alx - x + log(s/xpa)
52 C
53  RETURN
54  END
double precision function d9lgic(A, X, ALX)
Definition: d9lgic.f:2
octave_value log(void) const
Definition: ov.h:1190
subroutine xermsg(LIBRAR, SUBROU, MESSG, NERR, LEVEL)
Definition: xermsg.f:2
T abs(T x)
Definition: pr-output.cc:3062