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