 GNU Octave  4.2.1 A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
gamit.f
Go to the documentation of this file.
1 *DECK GAMIT
2  REAL FUNCTION gamit (A, X)
3 C***BEGIN PROLOGUE GAMIT
4 C***PURPOSE Calculate Tricomi's form of the incomplete Gamma function.
5 C***LIBRARY SLATEC (FNLIB)
6 C***CATEGORY C7E
7 C***TYPE SINGLE 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 GAMIT = 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 GAMIT is evaluated for arbitrary real values of A and for non-
22 C negative values of X (even though GAMIT is defined for X .LT.
23 C 0.0), except that for X = 0 and A .LE. 0.0, GAMIT is infinite,
24 C which is a fatal error.
25 C
26 C The function and both arguments are REAL.
27 C
28 C A slight deterioration of 2 or 3 digits accuracy will occur when
29 C GAMIT 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 ALGAMS, ALNGAM, GAMR, R1MACH, R9GMIT, R9LGIC,
42 C R9LGIT, 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 GAMIT
51  LOGICAL FIRST
52  SAVE alneps, sqeps, bot, first
53  DATA first /.true./
54 C***FIRST EXECUTABLE STATEMENT GAMIT
55  IF (first) THEN
56  alneps = -log(r1mach(3))
57  sqeps = sqrt(r1mach(4))
58  bot = log(r1mach(1))
59  ENDIF
60  first = .false.
61 C
62  IF (x .LT. 0.0) CALL xermsg('SLATEC', 'GAMIT', 'X IS NEGATIVE',
63  + 2, 2)
64 C
65  IF (x.NE.0.0) alx = log(x)
66  sga = 1.0
67  IF (a.NE.0.0) sga = sign(1.0, a)
68  ainta = aint(a+0.5*sga)
69  aeps = a - ainta
70 C
71  IF (x.GT.0.0) go to 20
72  gamit = 0.0
73  IF (ainta.GT.0.0 .OR. aeps.NE.0.0) gamit = gamr(a+1.0)
74  RETURN
75 C
76  20 IF (x.GT.1.0) go to 40
77  IF (a.GE.(-0.5) .OR. aeps.NE.0.0) CALL algams(a+1.0, algap1,
78  1 sgngam)
79  gamit = r9gmit(a, x, algap1, sgngam, alx)
80  RETURN
81 C
82  40 IF (a.LT.x) go to 50
83  t = r9lgit(a, x, alngam(a+1.0))
84  IF (t.LT.bot) CALL xerclr
85  gamit = exp(t)
86  RETURN
87 C
88  50 alng = r9lgic(a, x, alx)
89 C
90 C EVALUATE GAMIT IN TERMS OF LOG(GAMIC(A,X))
91 C
92  h = 1.0
93  IF (aeps.EQ.0.0 .AND. ainta.LE.0.0) go to 60
94  CALL algams(a+1.0, algap1, sgngam)
95  t = log(abs(a)) + alng - algap1
96  IF (t.GT.alneps) go to 70
97  IF (t.GT.(-alneps)) h = 1.0 - sga*sgngam*exp(t)
98  IF (abs(h).GT.sqeps) go to 60
99  CALL xerclr
100  CALL xermsg('SLATEC', 'GAMIT', 'RESULT LT HALF PRECISION', 1, 1)
101 C
102  60 t = -a*alx + log(abs(h))
103  IF (t.LT.bot) CALL xerclr
104  gamit = sign(exp(t), h)
105  RETURN
106 C
107  70 t = t - a*alx
108  IF (t.LT.bot) CALL xerclr
109  gamit = -sga*sgngam*exp(t)
110  RETURN
111 C
112  END
function r9lgit(A, X, ALGAP1)
Definition: r9lgit.f:2
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)
function gamr(X)
Definition: gamr.f:2
function r9gmit(A, X, ALGAP1, SGNGAM, ALX)
Definition: r9gmit.f:2
int exp(void) const
Definition: DET.h:66
subroutine algams(X, ALGAM, SGNGAM)
Definition: algams.f:2
is false
Definition: cellfun.cc:398
may be zero for pure relative error test tem the relative tolerance must be greater than or equal to