c8lgmc.f
1.92 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
53
54
55
56
57
58
59
60
61
62
63
64
65
66
complex function c8lgmc (zin)
c may 1978 edition. w. fullerton, c3, los alamos scientific lab.
c
c compute the log gamma correction term for almost all z so that
c clog(cgamma(z)) = 0.5*alog(2.*pi) + (z-0.5)*clog(z) - z + c8lgmc(z)
c this routine is valid even when real(z) is negative or when cabs(z)
c is small.
c when real(z) is negative, c8lgmc merely returns a correction which
c may be wrong by a multiple of 2*pi*i.
c
complex zin, z, corr, c9lgmc, c0lgmc, cexp, clnrel
external c0lgmc, c9lgmc, clnrel,
1 r1mach
data pi / 3.1415926535 8979324e0 /
data bound, sqeps, eps / 3*0.0 /
c
if (bound.ne.0.0) go to 10
nterm = -0.30*alog(r1mach(3))
bound = 0.1170*float(nterm)*
1 (0.1*r1mach(3))**(-1.0/(2.0*float(nterm)-1.0))
eps = 10.0*r1mach(4)
sqeps = sqrt (r1mach(4))
c
10 z = zin
x = real(z)
y = aimag(z)
if (x.lt.0.0) z = -z
c
corr = (0.0, 0.0)
if (x.gt.bound .or. abs(y).ge.bound) go to 30
cabsz = cabs(z)
if (cabsz.ge.bound) go to 30
c
if (x.eq.0.0 .and. y.eq.0.0) call seteru (
1 57hc8lgmc log gamma correction term is infinite for z = 0.0,
2 57, 3, 2)
c
n = sqrt (bound**2-y**2) - abs(x) + 1.0
do 20 i=1,n
corr = corr + c0lgmc(z)
z = z + 1.0
20 continue
c
30 c8lgmc = c9lgmc(z) + corr
if (x.gt.0.0) return
c
call entsrc (irold, 1)
test = 1.0
if (x.lt.(-0.5)) test = cabs ((z-aint(x-0.5))/z)
if (test.lt.eps) call seteru (
1 31hc8lgmc z is a negative integer, 31, 2, 2)
c
z = zin
if (y.lt.0.0) z = conjg(z)
corr = -cmplx(0.0,pi) + clnrel(-cexp(cmplx(0.,2.*pi)*z))
if (y.lt.0.0) corr = conjg(corr)
c
c8lgmc = corr - c8lgmc
c
call erroff
call entsrc (ir, irold)
if (test.lt.sqeps) call seteru ( 63hc8lgmc answer lt half precisi
1on, z too near a negative integer, 63, 1, 1)
c
return
end