function beta (a, b)
  c june 1977 edition.   w. fullerton, c3, los alamos scientific lab.
        external albeta,  gamma, r1mach
        data xmax, alnsml /0., 0./
  c
        if (xmax.ne.0.) go to 10
        call r9gaml (xmin, xmax)
        alnsml = alog(r1mach(1))
  c
   10   if (a.le.0. .or. b.le.0.) call seteru (
       1  35hbeta    both arguments must be gt 0, 35, 2, 2)
  c
        if (a+b.lt.xmax) beta = gamma(a) * gamma(b) / gamma(a+b)
        if (a+b.lt.xmax) return
  c
        beta = albeta (a, b)
        if (beta.lt.alnsml) call seteru (
       1  41hbeta    a and/or b so big beta underflows, 41, 1, 2)
  c
        beta = exp (beta)
  c
        return
        end