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