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