GNU Octave  4.0.0 A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
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
octave_value log(void) const
Definition: ov.h:1190
charNDArray max(char d, const charNDArray &m)
Definition: chNDArray.cc:233
subroutine xermsg(LIBRAR, SUBROU, MESSG, NERR, LEVEL)
Definition: xermsg.f:2
subroutine gamlim(XMIN, XMAX)
Definition: gamlim.f:2
T abs(T x)
Definition: pr-output.cc:3062