GNU Octave  4.4.1
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
d9gmit.f
Go to the documentation of this file.
1 *DECK D9GMIT
2  DOUBLE PRECISION FUNCTION d9gmit (A, X, ALGAP1, SGNGAM, ALX)
3 C***BEGIN PROLOGUE D9GMIT
4 C***SUBSIDIARY
5 C***PURPOSE Compute Tricomi's incomplete Gamma function for small
6 C arguments.
7 C***LIBRARY SLATEC (FNLIB)
8 C***CATEGORY C7E
9 C***TYPE DOUBLE PRECISION (R9GMIT-S, D9GMIT-D)
10 C***KEYWORDS COMPLEMENTARY INCOMPLETE GAMMA FUNCTION, FNLIB, SMALL X,
11 C SPECIAL FUNCTIONS, TRICOMI
12 C***AUTHOR Fullerton, W., (LANL)
13 C***DESCRIPTION
14 C
15 C Compute Tricomi's incomplete gamma function for small X.
16 C
17 C***REFERENCES (NONE)
18 C***ROUTINES CALLED D1MACH, DLNGAM, XERMSG
19 C***REVISION HISTORY (YYMMDD)
20 C 770701 DATE WRITTEN
21 C 890531 Changed all specific intrinsics to generic. (WRB)
22 C 890911 Removed unnecessary intrinsics. (WRB)
23 C 890911 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 D9GMIT
28  DOUBLE PRECISION a, x, algap1, sgngam, alx, ae, aeps, algs, alg2,
29  1 bot, eps, fk, s, sgng2, t, te, d1mach, dlngam
30  LOGICAL first
31  SAVE eps, bot, first
32  DATA first /.true./
33 C***FIRST EXECUTABLE STATEMENT D9GMIT
34  IF (first) THEN
35  eps = 0.5d0*d1mach(3)
36  bot = log(d1mach(1))
37  ENDIF
38  first = .false.
39 C
40  IF (x .LE. 0.d0) CALL xermsg ('SLATEC', 'D9GMIT',
41  + 'X SHOULD BE GT 0', 1, 2)
42 C
43  ma = a + 0.5d0
44  IF (a.LT.0.d0) ma = a - 0.5d0
45  aeps = a - ma
46 C
47  ae = a
48  IF (a.LT.(-0.5d0)) ae = aeps
49 C
50  t = 1.d0
51  te = ae
52  s = t
53  DO 20 k=1,200
54  fk = k
55  te = -x*te/fk
56  t = te/(ae+fk)
57  s = s + t
58  IF (abs(t).LT.eps*abs(s)) GO TO 30
59  20 CONTINUE
60  CALL xermsg ('SLATEC', 'D9GMIT',
61  + 'NO CONVERGENCE IN 200 TERMS OF TAYLOR-S SERIES', 2, 2)
62 C
63  30 IF (a.GE.(-0.5d0)) algs = -algap1 + log(s)
64  IF (a.GE.(-0.5d0)) GO TO 60
65 C
66  algs = -dlngam(1.d0+aeps) + log(s)
67  s = 1.0d0
68  m = -ma - 1
69  IF (m.EQ.0) GO TO 50
70  t = 1.0d0
71  DO 40 k=1,m
72  t = x*t/(aeps-(m+1-k))
73  s = s + t
74  IF (abs(t).LT.eps*abs(s)) GO TO 50
75  40 CONTINUE
76 C
77  50 d9gmit = 0.0d0
78  algs = -ma*log(x) + algs
79  IF (s.EQ.0.d0 .OR. aeps.EQ.0.d0) GO TO 60
80 C
81  sgng2 = sgngam * sign(1.0d0, s)
82  alg2 = -x - algap1 + log(abs(s))
83 C
84  IF (alg2.GT.bot) d9gmit = sgng2 * exp(alg2)
85  IF (algs.GT.bot) d9gmit = d9gmit + exp(algs)
86  RETURN
87 C
88  60 d9gmit = exp(algs)
89  RETURN
90 C
91  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)
double precision function d1mach(i)
Definition: d1mach.f:20
static T abs(T x)
Definition: pr-output.cc:1696
#define eps(C)