dchu.f
4.28 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
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
double precision function dchu (a, b, x)
c august 1980 edition. w. fullerton, c3, los alamos scientific lab.
c this routine is not valid when 1+a-b is close to zero if x is small.
double precision a, b, x, aintb, alnx, a0, beps, b0, c0, eps,
1 factor, gamri1, gamrni, pch1ai, pch1i, pi, pochai, sum, t,
2 xeps1, xi, xi1, xn, xtoeps, d1mach, dpoch, dgamma, dgamr,
3 dpoch1, dexprl, d9chu, dint, dexp, dlog, dsin, dsqrt
external d1mach, d9chu, dexprl, dgamma, dgamr,
1 dpoch, dpoch1
data pi / 3.1415926535 8979323846 2643383279 503 d0 /
data eps / 0.0d0 /
c
if (eps.eq.0.0d0) eps = d1mach(3)
c
if (x.lt.0.0d0) call seteru (31hdchu x is negative, use cchu,
1 31, 1, 2)
c
if (x.ne.0.0d0) go to 10
if (b.ge.1.0d0) call seteru (
1 35hchu x=0, b ge 1 so chu infinite, 35, 2, 2)
dchu = dgamma(1.0d0-b)/dgamma(1.0d0+a-b)
return
c
10 if (dmax1(dabs(a),1.0d0)*dmax1(dabs(1.0d0+a-b),1.0d0).lt.
1 0.99d0*dabs(x)) go to 120
c
c the ascending series will be used, because the descending rational
c approximation (which is based on the asymptotic series) is unstable.
c
if (dabs(1.0d0+a-b).lt.dsqrt(eps)) call seteru (
1 60hdchu algorithm is bad when 1+a-b is near zero for small x,
2 60, 10, 2)
c
if (b.ge.0.0d0) aintb = dint(b+0.5d0)
if (b.lt.0.0d0) aintb = dint(b-0.5d0)
beps = b - aintb
n = aintb
c
alnx = dlog(x)
xtoeps = dexp (-beps*alnx)
c
c evaluate the finite sum. -----------------------------------------
c
if (n.ge.1) go to 40
c
c consider the case b .lt. 1.0 first.
c
sum = 1.0d0
if (n.eq.0) go to 30
c
t = 1.0d0
m = -n
do 20 i=1,m
xi1 = i - 1
t = t*(a+xi1)*x/((b+xi1)*(xi1+1.0d0))
sum = sum + t
20 continue
c
30 sum = dpoch(1.0d0+a-b, -a)*sum
go to 70
c
c now consider the case b .ge. 1.0.
c
40 sum = 0.0d0
m = n - 2
if (m.lt.0) go to 70
t = 1.0d0
sum = 1.0d0
if (m.eq.0) go to 60
c
do 50 i=1,m
xi = i
t = t * (a-b+xi)*x/((1.0d0-b+xi)*xi)
sum = sum + t
50 continue
c
60 sum = dgamma(b-1.0d0) * dgamr(a) * x**(1-n) * xtoeps * sum
c
c next evaluate the infinite sum. ----------------------------------
c
70 istrt = 0
if (n.lt.1) istrt = 1 - n
xi = istrt
c
factor = (-1.0d0)**n * dgamr(1.0d0+a-b) * x**istrt
if (beps.ne.0.0d0) factor = factor * beps*pi/dsin(beps*pi)
c
pochai = dpoch (a, xi)
gamri1 = dgamr (xi+1.0d0)
gamrni = dgamr (aintb+xi)
b0 = factor * dpoch(a,xi-beps) * gamrni * dgamr(xi+1.0d0-beps)
c
if (dabs(xtoeps-1.0d0).gt.0.5d0) go to 90
c
c x**(-beps) is close to 1.0d0, so we must be careful in evaluating the
c differences.
c
pch1ai = dpoch1 (a+xi, -beps)
pch1i = dpoch1 (xi+1.0d0-beps, beps)
c0 = factor * pochai * gamrni * gamri1 * (
1 -dpoch1(b+xi,-beps) + pch1ai - pch1i + beps*pch1ai*pch1i)
c
c xeps1 = (1.0 - x**(-beps))/beps = (x**(-beps) - 1.0)/(-beps)
xeps1 = alnx*dexprl(-beps*alnx)
c
dchu = sum + c0 + xeps1*b0
xn = n
do 80 i=1,1000
xi = istrt + i
xi1 = istrt + i - 1
b0 = (a+xi1-beps)*b0*x/((xn+xi1)*(xi-beps))
c0 = (a+xi1)*c0*x/((b+xi1)*xi)
1 - ((a-1.0d0)*(xn+2.d0*xi-1.0d0) + xi*(xi-beps)) * b0
2 / (xi*(b+xi1)*(a+xi1-beps))
t = c0 + xeps1*b0
dchu = dchu + t
if (dabs(t).lt.eps*dabs(dchu)) go to 130
80 continue
call seteru (
1 60hdchu no convergence in 1000 terms of the ascending series,
2 60, 3, 2)
c
c x**(-beps) is very different from 1.0, so the straightforward
c formulation is stable.
c
90 a0 = factor * pochai * dgamr(b+xi) * gamri1 / beps
b0 = xtoeps * b0 / beps
c
dchu = sum + a0 - b0
do 100 i=1,1000
xi = istrt + i
xi1 = istrt + i - 1
a0 = (a+xi1)*a0*x/((b+xi1)*xi)
b0 = (a+xi1-beps)*b0*x/((aintb+xi1)*(xi-beps))
t = a0 - b0
dchu = dchu + t
if (dabs(t).lt.eps*dabs(dchu)) go to 130
100 continue
call seteru (
1 60hdchu no convergence in 1000 terms of the ascending series,
2 60, 3, 2)
c
c use luke-s rational approximation in the asymptotic region.
c
120 dchu = x**(-a) * d9chu(a,b,x)
c
130 return
end