1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
|
*DECK DGAMLM
SUBROUTINE DGAMLM (XMIN, XMAX)
C***BEGIN PROLOGUE DGAMLM
C***PURPOSE Compute the minimum and maximum bounds for the argument in
C the Gamma function.
C***LIBRARY SLATEC (FNLIB)
C***CATEGORY C7A, R2
C***TYPE DOUBLE PRECISION (GAMLIM-S, DGAMLM-D)
C***KEYWORDS COMPLETE GAMMA FUNCTION, FNLIB, LIMITS, SPECIAL FUNCTIONS
C***AUTHOR Fullerton, W., (LANL)
C***DESCRIPTION
C
C Calculate the minimum and maximum legal bounds for X in gamma(X).
C XMIN and XMAX are not the only bounds, but they are the only non-
C trivial ones to calculate.
C
C Output Arguments --
C XMIN double precision minimum legal value of X in gamma(X). Any
C smaller value of X might result in underflow.
C XMAX double precision maximum legal value of X in gamma(X). Any
C larger value of X might cause overflow.
C
C***REFERENCES (NONE)
C***ROUTINES CALLED D1MACH, XERMSG
C***REVISION HISTORY (YYMMDD)
C 770601 DATE WRITTEN
C 890531 Changed all specific intrinsics to generic. (WRB)
C 890531 REVISION DATE from Version 3.2
C 891214 Prologue converted to Version 4.0 format. (BAB)
C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
C***END PROLOGUE DGAMLM
DOUBLE PRECISION XMIN, XMAX, ALNBIG, ALNSML, XLN, XOLD, D1MACH
C***FIRST EXECUTABLE STATEMENT DGAMLM
ALNSML = LOG(D1MACH(1))
XMIN = -ALNSML
DO 10 I=1,10
XOLD = XMIN
XLN = LOG(XMIN)
XMIN = XMIN - XMIN*((XMIN+0.5D0)*XLN - XMIN - 0.2258D0 + ALNSML)
1 / (XMIN*XLN+0.5D0)
IF (ABS(XMIN-XOLD).LT.0.005D0) GO TO 20
10 CONTINUE
CALL XERMSG ('SLATEC', 'DGAMLM', 'UNABLE TO FIND XMIN', 1, 2)
C
20 XMIN = -XMIN + 0.01D0
C
ALNBIG = LOG (D1MACH(2))
XMAX = ALNBIG
DO 30 I=1,10
XOLD = XMAX
XLN = LOG(XMAX)
XMAX = XMAX - XMAX*((XMAX-0.5D0)*XLN - XMAX + 0.9189D0 - ALNBIG)
1 / (XMAX*XLN-0.5D0)
IF (ABS(XMAX-XOLD).LT.0.005D0) GO TO 40
30 CONTINUE
CALL XERMSG ('SLATEC', 'DGAMLM', 'UNABLE TO FIND XMAX', 2, 2)
C
40 XMAX = XMAX - 0.01D0
XMIN = MAX (XMIN, -XMAX+1.D0)
C
RETURN
END
|
