Blame view
fvn_fnlib/chu.f
3.81 KB
38581db0c
git-svn-id: https...
|
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 |
function chu (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.
c
external exprel, gamma, gamr, poch, poch1,
1 r1mach, r9chu
data pi / 3.1415926535 8979324 e0 /
data eps / 0.0 /
c
if (eps.eq.0.0) eps = r1mach(3)
c
if (x.lt.0.0) call seteru (31hchu x is negative, use cchu,
1 31, 1, 2)
c
if (x.ne.0.0) go to 10
if (b.ge.1.0) call seteru (
1 35hchu x=0, b ge 1 so chu infinite, 35, 2, 2)
chu = gamma(1.0-b)/gamma(1.0+a-b)
return
c
10 if (amax1(abs(a),1.0)*amax1(abs(1.0+a-b),1.0).lt.0.99*abs(x))
1 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 (abs(1.0+a-b).lt.sqrt(eps)) call seteru (
1 60hchu algorithm is bad when 1+a-b is near zero for small x,
2 60, 10, 2)
c
aintb = aint(b+0.5)
if (b.lt.0.0) aintb = aint(b-0.5)
beps = b - aintb
n = aintb
c
alnx = alog(x)
xtoeps = exp(-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.0
if (n.eq.0) go to 30
c
t = 1.0
m = -n
do 20 i=1,m
xi1 = i - 1
t = t*(a+xi1)*x/((b+xi1)*(xi1+1.0))
sum = sum + t
20 continue
c
30 sum = poch(1.0+a-b, -a) * sum
go to 70
c
c now consider the case b .ge. 1.0.
c
40 sum = 0.0
m = n - 2
if (m.lt.0) go to 70
t = 1.0
sum = 1.0
if (m.eq.0) go to 60
c
do 50 i=1,m
xi = i
t = t * (a-b+xi)*x/((1.0-b+xi)*xi)
sum = sum + t
50 continue
c
60 sum = gamma(b-1.0) * gamr(a) * x**(1-n) * xtoeps * sum
c
c now evaluate the infinite sum. -----------------------------------
c
70 istrt = 0
if (n.lt.1) istrt = 1 - n
xi = istrt
c
factor = (-1.0)**n * gamr(1.0+a-b) * x**istrt
if (beps.ne.0.0) factor = factor * beps*pi/sin(beps*pi)
c
pochai = poch (a, xi)
gamri1 = gamr (xi+1.0)
gamrni = gamr (aintb+xi)
b0 = factor * poch(a,xi-beps) * gamrni * gamr(xi+1.0-beps)
c
if (abs(xtoeps-1.0).gt.0.5) go to 90
c
c x**(-beps) is close to 1.0, so we must be careful in evaluating
c the differences
c
pch1ai = poch1 (a+xi, -beps)
pch1i = poch1 (xi+1.0-beps, beps)
c0 = factor * pochai * gamrni * gamri1 * (
1 -poch1(b+xi, -beps) + pch1ai - pch1i + beps*pch1ai*pch1i )
c
c xeps1 = (1.0 - x**(-beps)) / beps
xeps1 = alnx * exprel(-beps*alnx)
c
chu = 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) - ((a-1.0)*(xn+2.*xi-1.0)
1 + xi*(xi-beps)) * b0/(xi*(b+xi1)*(a+xi1-beps))
t = c0 + xeps1*b0
chu = chu + t
if (abs(t).lt.eps*abs(chu)) go to 130
80 continue
call seteru (
1 60hchu 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 * gamr(b+xi) * gamri1 / beps
b0 = xtoeps*b0/beps
c
chu = 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
chu = chu + t
if (abs(t).lt.eps*abs(chu)) go to 130
100 continue
call seteru (
1 60hchu no convergence in 1000 terms of the ascending series,
2 60, 3, 2)
c
c use luke-s rational approx in the asymptotic region.
c
120 chu = x**(-a) * r9chu(a, b, x)
c
130 return
end
|