r9chu.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
67
68
69
70
function r9chu (a, b, z)
c august 1977 edition. w. fullerton, c3, los alamos scientific lab.
c
c evaluate for large z z**a * u(a,b,z) where u is the logarithmic
c confluent hypergeometric function. a rational approximation due to y.
c l. luke is used. when u is not in the asymptotic region, i.e., when a
c or b is large compared with z, considerable significance loss occurs.
c a warning is provided when the computed result is less than half
c precision.
c
dimension aa(4), bb(4)
external r1mach
data eps, sqeps / 2*0.0 /
c
if (eps.ne.0.0) go to 10
eps = 4.0*r1mach(4)
sqeps = sqrt (r1mach(4))
c
10 bp = 1.0 + a - b
ab = a*bp
ct2 = 2.0*(z-ab)
sab = a + bp
c
bb(1) = 1.0
aa(1) = 1.0
c
ct3 = sab + 1.0 + ab
bb(2) = 1.0 + 2.0*z/ct3
aa(2) = 1.0 + ct2/ct3
c
anbn = ct3 + sab + 3.0
ct1 = 1.0 + 2.0*z/anbn
bb(3) = 1.0 + 6.0*ct1*z/ct3
aa(3) = 1.0 + 6.0*ab/anbn + 3.0*ct1*ct2/ct3
c
do 30 i=4,300
x2i1 = 2*i - 3
ct1 = x2i1/(x2i1-2.0)
anbn = anbn + x2i1 + sab
ct2 = (x2i1 - 1.0) / anbn
c2 = x2i1*ct2 - 1.0
d1z = x2i1*2.0*z/anbn
c
ct3 = sab*ct2
g1 = d1z + ct1*(c2+ct3)
g2 = d1z - c2
g3 = ct1*(1.0 - ct3 - 2.0*ct2)
c
bb(4) = g1*bb(3) + g2*bb(2) + g3*bb(1)
aa(4) = g1*aa(3) + g2*aa(2) + g3*aa(1)
c
r9chu = aa(4)/bb(4)
if (abs(r9chu-aa(1)/bb(1)).lt.eps*abs(r9chu)) go to 40
c
c if overflows or underflows prove to be a problem, the statements
c below could be altered to incorporate a dynamically adjusted scale
c factor.
c
do 20 j=1,3
bb(j) = bb(j+1)
aa(j) = aa(j+1)
20 continue
30 continue
call seteru (35hr9chu no convergence in 300 terms, 35, 1, 2)
c
40 if (r9chu.lt.sqeps .or. r9chu.gt.1.0/sqeps) call seteru (
1 32hr9chu answer lt half precision, 32, 2, 1)
c
return
end