GNU Octave  4.0.0 A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
dgamit.f
Go to the documentation of this file.
1 *DECK DGAMIT
2  DOUBLE PRECISION FUNCTION dgamit (A, X)
3 C***BEGIN PROLOGUE DGAMIT
4 C***PURPOSE Calculate Tricomi's form of the incomplete Gamma function.
5 C***LIBRARY SLATEC (FNLIB)
6 C***CATEGORY C7E
7 C***TYPE DOUBLE PRECISION (GAMIT-S, DGAMIT-D)
8 C***KEYWORDS COMPLEMENTARY INCOMPLETE GAMMA FUNCTION, FNLIB,
9 C SPECIAL FUNCTIONS, TRICOMI
10 C***AUTHOR Fullerton, W., (LANL)
11 C***DESCRIPTION
12 C
13 C Evaluate Tricomi's incomplete Gamma function defined by
14 C
15 C DGAMIT = X**(-A)/GAMMA(A) * integral from 0 to X of EXP(-T) *
16 C T**(A-1.)
17 C
18 C for A .GT. 0.0 and by analytic continuation for A .LE. 0.0.
19 C GAMMA(X) is the complete gamma function of X.
20 C
21 C DGAMIT is evaluated for arbitrary real values of A and for non-
22 C negative values of X (even though DGAMIT is defined for X .LT.
23 C 0.0), except that for X = 0 and A .LE. 0.0, DGAMIT is infinite,
24 C which is a fatal error.
25 C
26 C The function and both arguments are DOUBLE PRECISION.
27 C
28 C A slight deterioration of 2 or 3 digits accuracy will occur when
29 C DGAMIT is very large or very small in absolute value, because log-
30 C arithmic variables are used. Also, if the parameter A is very
31 C close to a negative integer (but not a negative integer), there is
32 C a loss of accuracy, which is reported if the result is less than
33 C half machine precision.
34 C
35 C***REFERENCES W. Gautschi, A computational procedure for incomplete
36 C gamma functions, ACM Transactions on Mathematical
37 C Software 5, 4 (December 1979), pp. 466-481.
38 C W. Gautschi, Incomplete gamma functions, Algorithm 542,
39 C ACM Transactions on Mathematical Software 5, 4
40 C (December 1979), pp. 482-489.
41 C***ROUTINES CALLED D1MACH, D9GMIT, D9LGIC, D9LGIT, DGAMR, DLGAMS,
42 C DLNGAM, XERCLR, XERMSG
43 C***REVISION HISTORY (YYMMDD)
44 C 770701 DATE WRITTEN
45 C 890531 Changed all specific intrinsics to generic. (WRB)
46 C 890531 REVISION DATE from Version 3.2
47 C 891214 Prologue converted to Version 4.0 format. (BAB)
48 C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
49 C 920528 DESCRIPTION and REFERENCES sections revised. (WRB)
50 C***END PROLOGUE DGAMIT
51  DOUBLE PRECISION A, X, AEPS, AINTA, ALGAP1, ALNEPS, ALNG, ALX,
52  1 bot, h, sga, sgngam, sqeps, t, d1mach, dgamr, d9gmit, d9lgit,
53  2 dlngam, d9lgic
54  LOGICAL FIRST
55  SAVE alneps, sqeps, bot, first
56  DATA first /.true./
57 C***FIRST EXECUTABLE STATEMENT DGAMIT
58  IF (first) THEN
59  alneps = -log(d1mach(3))
60  sqeps = sqrt(d1mach(4))
61  bot = log(d1mach(1))
62  ENDIF
63  first = .false.
64 C
65  IF (x .LT. 0.d0) CALL xermsg('SLATEC', 'DGAMIT', 'X IS NEGATIVE'
66  + , 2, 2)
67 C
68  IF (x.NE.0.d0) alx = log(x)
69  sga = 1.0d0
70  IF (a.NE.0.d0) sga = sign(1.0d0, a)
71  ainta = aint(a + 0.5d0*sga)
72  aeps = a - ainta
73 C
74  IF (x.GT.0.d0) go to 20
75  dgamit = 0.0d0
76  IF (ainta.GT.0.d0 .OR. aeps.NE.0.d0) dgamit = dgamr(a+1.0d0)
77  RETURN
78 C
79  20 IF (x.GT.1.d0) go to 30
80  IF (a.GE.(-0.5d0) .OR. aeps.NE.0.d0) CALL dlgams(a+1.0d0, algap1,
81  1 sgngam)
82  dgamit = d9gmit(a, x, algap1, sgngam, alx)
83  RETURN
84 C
85  30 IF (a.LT.x) go to 40
86  t = d9lgit(a, x, dlngam(a+1.0d0))
87  IF (t.LT.bot) CALL xerclr
88  dgamit = exp(t)
89  RETURN
90 C
91  40 alng = d9lgic(a, x, alx)
92 C
93 C EVALUATE DGAMIT IN TERMS OF LOG (DGAMIC (A, X))
94 C
95  h = 1.0d0
96  IF (aeps.EQ.0.d0 .AND. ainta.LE.0.d0) go to 50
97 C
98  CALL dlgams(a+1.0d0, algap1, sgngam)
99  t = log(abs(a)) + alng - algap1
100  IF (t.GT.alneps) go to 60
101 C
102  IF (t.GT.(-alneps)) h = 1.0d0 - sga * sgngam * exp(t)
103  IF (abs(h).GT.sqeps) go to 50
104 C
105  CALL xerclr
106  CALL xermsg('SLATEC', 'DGAMIT', 'RESULT LT HALF PRECISION', 1,
107  + 1)
108 C
109  50 t = -a*alx + log(abs(h))
110  IF (t.LT.bot) CALL xerclr
111  dgamit = sign(exp(t), h)
112  RETURN
113 C
114  60 t = t - a*alx
115  IF (t.LT.bot) CALL xerclr
116  dgamit = -sga * sgngam * exp(t)
117  RETURN
118 C
119  END
subroutine dlgams(X, DLGAM, SGNGAM)
Definition: dlgams.f:2
double precision function dgamr(X)
Definition: dgamr.f:2
double precision function d9lgic(A, X, ALX)
Definition: d9lgic.f:2
double precision function d9gmit(A, X, ALGAP1, SGNGAM, ALX)
Definition: d9gmit.f:2
int exp(void) const
Definition: DET.h:64
octave_value log(void) const
Definition: ov.h:1190
subroutine xermsg(LIBRAR, SUBROU, MESSG, NERR, LEVEL)
Definition: xermsg.f:2
double precision function d9lgit(A, X, ALGAP1)
Definition: d9lgit.f:2
subroutine xerclr
Definition: xerclr.f:2
T abs(T x)
Definition: pr-output.cc:3062
octave_value sqrt(void) const
Definition: ov.h:1200