clngam.f
2.55 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
67
68
69
70
71
72
73
74
75
76
77
78
79
complex function clngam (zin)
c august 1980 edition. w. fullerton c3, los alamos scientific lab.
c eventually clngam should make use of c8lgmc for all z except for
c z in the vicinity of 1 and 2.
complex zin, z, corr, cexp, clog, clnrel, c9lgmc
external c9lgmc, carg, clnrel,
1 r1mach
data pi / 3.1415926535 8979324e0 /
data sq2pil / 0.9189385332 0467274e0 /
c
data bound, dxrel, rmax / 3*0.0 /
c
if (bound.ne.0.) go to 10
n = -0.30*alog(r1mach(3))
c bound = n*(0.1*eps)**(-1/(2*n-1))/(pi*exp(1))
bound = 0.1171*float(n)*(0.1*r1mach(3))**(-1./(2.*float(n)-1.))
dxrel = sqrt (r1mach(4))
rmax = r1mach(2)/alog(r1mach(2))
c
10 z = zin
x = real(zin)
y = aimag(zin)
c
corr = (0.0, 0.0)
cabsz = cabs(z)
if (cabsz.gt.rmax) call seteru (
1 44hclngam cabs(z) so large result may overflow, 44, 2, 2)
if (x.ge.0.0 .and. cabsz.gt.bound) go to 50
if (x.lt.0.0 .and. abs(y).gt.bound) go to 50
c
if (cabsz.lt.bound) go to 20
c
c use the reflection formula for real(z) negative, cabs(z) large, and
c abs(aimag(y)) small.
c
call entsrc (irold, 1)
if (y.gt.0.0) z = conjg (z)
corr = cexp (-cmplx(0.0,2.0*pi)*z)
if (real(corr).eq.1.0 .and. aimag(corr).eq.0.0) call seteru (
1 31hclngam z is a negative integer, 31, 3, 2)
c
clngam = sq2pil + 1.0 - cmplx(0.0,pi)*(z-0.5) - clnrel(-corr)
1 + (z-0.5)*clog(1.0-z) - z - c9lgmc(1.0-z)
if (y.gt.0.0) clngam = conjg (clngam)
c
call erroff
call entsrc (ir, irold)
if (abs(y).gt.dxrel) return
c X in sequence field on next line is to preserve trailing blank
if (0.5*cabs((1.-corr)*clngam/z).lt.dxrel) call seteru (68hclngam X
1 answer lt half precision because z too near negative integer, 68,
2 1, 1)
return
c
c use the recursion relation for cabs(z) small.
c
20 if (x.ge.(-0.5) .or. abs(y).gt.dxrel) go to 30
if (cabs((z-aint(x-0.5))/x).lt.dxrel) call seteru ( 68hclngam ans
1wer lt half precision because z too near negative integer, 68,1,1)
c
30 n = sqrt (bound**2 - y**2) - x + 1.0
argsum = 0.0
corr = (1.0, 0.0)
do 40 i=1,n
argsum = argsum + carg(z)
corr = z*corr
z = 1.0 + z
40 continue
c
if (real(corr).eq.0.0 .and. aimag(corr).eq.0.0) call seteru (
1 31hclngam z is a negative integer, 31, 3, 2)
corr = -cmplx (alog(cabs(corr)), argsum)
c
c use stirling-s approximation for large z.
c
50 clngam = sq2pil + (z-0.5)*clog(z) - z + corr + c9lgmc(z)
return
c
end