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 aint, alog, exp, exprel, gamma, gamr, poch, poch1,
     1  r1mach, r9chu, sin, sqrt
      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