d9gaml.f
1.42 KB
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
subroutine d9gaml (xmin, xmax)
c june 1977 edition. w. fullerton, c3, los alamos scientific lab.
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 dble prec minimum legal value of x in gamma(x). any smaller
c value of x might result in underflow.
c xmax dble prec maximum legal value of x in gamma(x). any larger
c value of x might cause overflow.
c
double precision xmin, xmax, alnbig, alnsml, xln, xold, d1mach,
1 dlog
external d1mach
c
alnsml = dlog(d1mach(1))
xmin = -alnsml
do 10 i=1,10
xold = xmin
xln = dlog(xmin)
xmin = xmin - xmin*((xmin+0.5d0)*xln - xmin - 0.2258d0 + alnsml)
1 / (xmin*xln+0.5d0)
if (dabs(xmin-xold).lt.0.005d0) go to 20
10 continue
call seteru (27hd9gaml unable to find xmin, 27, 1, 2)
c
20 xmin = -xmin + 0.01d0
c
alnbig = dlog (d1mach(2))
xmax = alnbig
do 30 i=1,10
xold = xmax
xln = dlog(xmax)
xmax = xmax - xmax*((xmax-0.5d0)*xln - xmax + 0.9189d0 - alnbig)
1 / (xmax*xln-0.5d0)
if (dabs(xmax-xold).lt.0.005d0) go to 40
30 continue
call seteru (27hd9gaml unable to find xmax, 27, 2, 2)
c
40 xmax = xmax - 0.01d0
xmin = dmax1 (xmin, -xmax+1.d0)
c
return
end