double precision function dlbeta (a, b)
  c july 1977 edition.   w. fullerton, c3, los alamos scientific lab.
        double precision a, b, p, q, corr, sq2pil, d9lgmc, dgamma, dlngam,
       1  dlnrel, dlog
        external d9lgmc, dgamma, dlngam, dlnrel
        data sq2pil / 0.9189385332 0467274178 0329736405 62 d0 /
  c
        p = dmin1 (a, b)
        q = dmax1 (a, b)
  c
        if (p.le.0.d0) call seteru (
       1  38hdlbeta  both arguments must be gt zero, 38, 1, 2)
  c
        if (p.ge.10.d0) go to 30
        if (q.ge.10.d0) go to 20
  c
  c p and q are small.
  c
        dlbeta = dlog (dgamma(p) * (dgamma(q)/dgamma(p+q)) )
        return
  c
  c p is small, but q is big.
  c
   20   corr = d9lgmc(q) - d9lgmc(p+q)
        dlbeta = dlngam(p) + corr + p - p*dlog(p+q)
       1  + (q-0.5d0)*dlnrel(-p/(p+q))
        return
  c
  c p and q are big.
  c
   30   corr = d9lgmc(p) + d9lgmc(q) - d9lgmc(p+q)
        dlbeta = -0.5d0*dlog(q) + sq2pil + corr + (p-0.5d0)*dlog(p/(p+q))
       1  + q*dlnrel(-p/(p+q))
        return
  c
        end