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