GNU Octave  3.8.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
gamlim.f
Go to the documentation of this file.
1 *DECK GAMLIM
2  SUBROUTINE gamlim (XMIN, XMAX)
3 C***BEGIN PROLOGUE GAMLIM
4 C***PURPOSE Compute the minimum and maximum bounds for the argument in
5 C the Gamma function.
6 C***LIBRARY SLATEC (FNLIB)
7 C***CATEGORY C7A, R2
8 C***TYPE SINGLE PRECISION (GAMLIM-S, DGAMLM-D)
9 C***KEYWORDS COMPLETE GAMMA FUNCTION, FNLIB, LIMITS, SPECIAL FUNCTIONS
10 C***AUTHOR Fullerton, W., (LANL)
11 C***DESCRIPTION
12 C
13 C Calculate the minimum and maximum legal bounds for X in GAMMA(X).
14 C XMIN and XMAX are not the only bounds, but they are the only non-
15 C trivial ones to calculate.
16 C
17 C Output Arguments --
18 C XMIN minimum legal value of X in GAMMA(X). Any smaller value of
19 C X might result in underflow.
20 C XMAX maximum legal value of X in GAMMA(X). Any larger value will
21 C cause overflow.
22 C
23 C***REFERENCES (NONE)
24 C***ROUTINES CALLED R1MACH, XERMSG
25 C***REVISION HISTORY (YYMMDD)
26 C 770401 DATE WRITTEN
27 C 890531 Changed all specific intrinsics to generic. (WRB)
28 C 890531 REVISION DATE from Version 3.2
29 C 891214 Prologue converted to Version 4.0 format. (BAB)
30 C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
31 C***END PROLOGUE GAMLIM
32 C***FIRST EXECUTABLE STATEMENT GAMLIM
33  alnsml = log(r1mach(1))
34  xmin = -alnsml
35  DO 10 i=1,10
36  xold = xmin
37  xln = log(xmin)
38  xmin = xmin - xmin*((xmin+0.5)*xln - xmin - 0.2258 + alnsml)
39  1 / (xmin*xln + 0.5)
40  IF (abs(xmin-xold).LT.0.005) go to 20
41  10 CONTINUE
42  CALL xermsg('SLATEC', 'GAMLIM', 'UNABLE TO FIND XMIN', 1, 2)
43 C
44  20 xmin = -xmin + 0.01
45 C
46  alnbig = log(r1mach(2))
47  xmax = alnbig
48  DO 30 i=1,10
49  xold = xmax
50  xln = log(xmax)
51  xmax = xmax - xmax*((xmax-0.5)*xln - xmax + 0.9189 - alnbig)
52  1 / (xmax*xln - 0.5)
53  IF (abs(xmax-xold).LT.0.005) go to 40
54  30 CONTINUE
55  CALL xermsg('SLATEC', 'GAMLIM', 'UNABLE TO FIND XMAX', 2, 2)
56 C
57  40 xmax = xmax - 0.01
58  xmin = max(xmin, -xmax+1.)
59 C
60  RETURN
61  END