complex(kind(1.d0)) function zbeta (a, b)

        implicit none
  c july 1977 edition.  w. fullerton, c3, los alamos scientific lab.

        complex(kind(1.d0)) a, b, zgamma, zlbeta
        real(kind(1.d0)) xmax,xmin

        external  zgamma, zlbeta
        data xmax / 0.0 /
  c
        if (xmax.eq.0.0) call d9gaml (xmin, xmax)
  c
        if (real(a).le.0.0 .or. real(b).le.0.0) call seteru (
       1  48hzbeta   real part of both arguments must be gt 0, 48, 1, 2)
  c
        if (real(a)+real(b).lt.xmax) zbeta = zgamma(a) * (zgamma(b)/
       1  zgamma(a+b) )
        if (real(a)+real(b).lt.xmax) return
  c
        zbeta = exp (zlbeta(a, b))
  c
        return
        end