double precision function dbeta (a, b)
  c june 1977 edition.   w. fullerton, c3, los alamos scientific lab.
        double precision a, b, alnsml, xmax, xmin,
       1  d1mach, dexp, dgamma, dlbeta, dlog
        external d1mach, dgamma, dlbeta
        data xmax, alnsml / 0.d0, 0.d0 /
  c
        if (xmax.ne.0.d0) go to 10
        call d9gaml (xmin, xmax)
        alnsml = dlog (d1mach(1))
  c
   10   if (a.le.0.d0 .or. b.le.0.d0) call seteru (
       1  35hdbeta   both arguments must be gt 0, 35, 2, 2)
  c
        if (a+b.lt.xmax) dbeta = dgamma(a)*dgamma(b)/dgamma(a+b)
        if (a+b.lt.xmax) return
  c
        dbeta = dlbeta (a, b)
        if (dbeta.lt.alnsml) go to 20
        dbeta = dexp (dbeta)
        return
  c
   20   dbeta = 0.d0
        call seteru (41hdbeta   a and/or b so big beta underflows, 41, 1,
       1  1)
        return
  c
        end