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
dgamlm.f
Go to the documentation of this file.
1 *DECK DGAMLM
2  SUBROUTINE dgamlm (XMIN, XMAX)
3 C***BEGIN PROLOGUE DGAMLM
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 DOUBLE 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 double precision minimum legal value of X in gamma(X). Any
19 C smaller value of X might result in underflow.
20 C XMAX double precision maximum legal value of X in gamma(X). Any
21 C larger value of X might cause overflow.
22 C
23 C***REFERENCES (NONE)
24 C***ROUTINES CALLED D1MACH, XERMSG
25 C***REVISION HISTORY (YYMMDD)
26 C 770601 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 DGAMLM
32  DOUBLE PRECISION xmin, xmax, alnbig, alnsml, xln, xold, d1mach
33 C***FIRST EXECUTABLE STATEMENT DGAMLM
34  alnsml = log(d1mach(1))
35  xmin = -alnsml
36  DO 10 i=1,10
37  xold = xmin
38  xln = log(xmin)
39  xmin = xmin - xmin*((xmin+0.5d0)*xln - xmin - 0.2258d0 + alnsml)
40  1 / (xmin*xln+0.5d0)
41  IF (abs(xmin-xold).LT.0.005d0) go to 20
42  10 CONTINUE
43  CALL xermsg('SLATEC', 'DGAMLM', 'UNABLE TO FIND XMIN', 1, 2)
44 C
45  20 xmin = -xmin + 0.01d0
46 C
47  alnbig = log(d1mach(2))
48  xmax = alnbig
49  DO 30 i=1,10
50  xold = xmax
51  xln = log(xmax)
52  xmax = xmax - xmax*((xmax-0.5d0)*xln - xmax + 0.9189d0 - alnbig)
53  1 / (xmax*xln-0.5d0)
54  IF (abs(xmax-xold).LT.0.005d0) go to 40
55  30 CONTINUE
56  CALL xermsg('SLATEC', 'DGAMLM', 'UNABLE TO FIND XMAX', 2, 2)
57 C
58  40 xmax = xmax - 0.01d0
59  xmin = max(xmin, -xmax+1.d0)
60 C
61  RETURN
62  END