GNU Octave  6.2.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
r9gmit.f
Go to the documentation of this file.
1 *DECK R9GMIT
2  FUNCTION r9gmit (A, X, ALGAP1, SGNGAM, ALX)
3 C***BEGIN PROLOGUE R9GMIT
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 SINGLE 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 ALNGAM, R1MACH, XERMSG
19 C***REVISION HISTORY (YYMMDD)
20 C 770701 DATE WRITTEN
21 C 890531 Changed all specific intrinsics to generic. (WRB)
22 C 890531 REVISION DATE from Version 3.2
23 C 891214 Prologue converted to Version 4.0 format. (BAB)
24 C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
25 C 900720 Routine changed from user-callable to subsidiary. (WRB)
26 C***END PROLOGUE R9GMIT
27  SAVE eps, bot
28  DATA eps, bot / 2*0.0 /
29 C***FIRST EXECUTABLE STATEMENT R9GMIT
30  IF (eps.EQ.0.0) eps = 0.5*r1mach(3)
31  IF (bot.EQ.0.0) bot = log(r1mach(1))
32 C
33  IF (x .LE. 0.0) CALL xermsg ('SLATEC', 'R9GMIT',
34  + 'X SHOULD BE GT 0', 1, 2)
35 C
36  ma = a + 0.5
37  IF (a.LT.0.0) ma = a - 0.5
38  aeps = a - ma
39 C
40  ae = a
41  IF (a.LT.(-0.5)) ae = aeps
42 C
43  t = 1.0
44  te = ae
45  s = t
46  DO 20 k=1,200
47  fk = k
48  te = -x*te/fk
49  t = te/(ae+fk)
50  s = s + t
51  IF (abs(t).LT.eps*abs(s)) GO TO 30
52  20 CONTINUE
53  CALL xermsg ('SLATEC', 'R9GMIT',
54  + 'NO CONVERGENCE IN 200 TERMS OF TAYLOR-S SERIES', 2, 2)
55 C
56  30 IF (a.GE.(-0.5)) algs = -algap1 + log(s)
57  IF (a.GE.(-0.5)) GO TO 60
58 C
59  algs = -alngam(1.0+aeps) + log(s)
60  s = 1.0
61  m = -ma - 1
62  IF (m.EQ.0) GO TO 50
63  t = 1.0
64  DO 40 k=1,m
65  t = x*t/(aeps-m-1+k)
66  s = s + t
67  IF (abs(t).LT.eps*abs(s)) GO TO 50
68  40 CONTINUE
69 C
70  50 r9gmit = 0.0
71  algs = -ma*log(x) + algs
72  IF (s.EQ.0.0 .OR. aeps.EQ.0.0) GO TO 60
73 C
74  sgng2 = sgngam*sign(1.0,s)
75  alg2 = -x - algap1 + log(abs(s))
76 C
77  IF (alg2.GT.bot) r9gmit = sgng2*exp(alg2)
78  IF (algs.GT.bot) r9gmit = r9gmit + exp(algs)
79  RETURN
80 C
81  60 r9gmit = exp(algs)
82  RETURN
83 C
84  END
#define eps(C)
function alngam(X)
Definition: alngam.f:3
static T abs(T x)
Definition: pr-output.cc:1678
real function r1mach(i)
Definition: r1mach.f:23
function r9gmit(A, X, ALGAP1, SGNGAM, ALX)
Definition: r9gmit.f:3
subroutine xermsg(LIBRAR, SUBROU, MESSG, NERR, LEVEL)
Definition: xermsg.f:3