d9lgmc.f
2.36 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
47
48
49
50
51
52
double precision function d9lgmc (x)
c august 1977 edition. w. fullerton, c3, los alamos scientific lab.
c
c compute the log gamma correction factor for x .ge. 10. so that
c dlog (dgamma(x)) = dlog(dsqrt(2*pi)) + (x-.5)*dlog(x) - x + d9lgmc(x)
c
double precision x, algmcs(15), xbig, xmax, dcsevl, d1mach,
1 dexp, dlog, dsqrt
external d1mach, dcsevl, initds
c
c series for algm on the interval 0. to 1.00000e-02
c with weighted error 1.28e-31
c log weighted error 30.89
c significant figures required 29.81
c decimal places required 31.48
c
data algmcs( 1) / +.1666389480 4518632472 0572965082 2 d+0 /
data algmcs( 2) / -.1384948176 0675638407 3298605913 5 d-4 /
data algmcs( 3) / +.9810825646 9247294261 5717154748 7 d-8 /
data algmcs( 4) / -.1809129475 5724941942 6330626671 9 d-10 /
data algmcs( 5) / +.6221098041 8926052271 2601554341 6 d-13 /
data algmcs( 6) / -.3399615005 4177219443 0333059966 6 d-15 /
data algmcs( 7) / +.2683181998 4826987489 5753884666 6 d-17 /
data algmcs( 8) / -.2868042435 3346432841 4462239999 9 d-19 /
data algmcs( 9) / +.3962837061 0464348036 7930666666 6 d-21 /
data algmcs( 10) / -.6831888753 9857668701 1199999999 9 d-23 /
data algmcs( 11) / +.1429227355 9424981475 7333333333 3 d-24 /
data algmcs( 12) / -.3547598158 1010705471 9999999999 9 d-26 /
data algmcs( 13) / +.1025680058 0104709120 0000000000 0 d-27 /
data algmcs( 14) / -.3401102254 3167487999 9999999999 9 d-29 /
data algmcs( 15) / +.1276642195 6300629333 3333333333 3 d-30 /
c
data nalgm, xbig, xmax / 0, 2*0.d0 /
c
if (nalgm.ne.0) go to 10
nalgm = initds (algmcs, 15, sngl(d1mach(3)) )
xbig = 1.0d0/dsqrt(d1mach(3))
xmax = dexp (dmin1(dlog(d1mach(2)/12.d0), -dlog(12.d0*d1mach(1))))
c
10 if (x.lt.10.d0) call seteru (23hd9lgmc x must be ge 10, 23, 1, 2)
if (x.ge.xmax) go to 20
c
d9lgmc = 1.d0/(12.d0*x)
if (x.lt.xbig) d9lgmc = dcsevl (2.0d0*(10.d0/x)**2-1.d0, algmcs,
1 nalgm) / x
return
c
20 d9lgmc = 0.d0
call seteru (34hd9lgmc x so big d9lgmc underflows, 34, 2, 0)
return
c
end