function albeta (a, b)
  c july 1977 edition.   w. fullerton, c3, los alamos scientific lab.
        external alngam, alnrel, gamma, r9lgmc
        data sq2pil / 0.9189385332 0467274 e0 /
  c
        p = amin1 (a, b)
        q = amax1 (a, b)
  c
        if (p.le.0.0) call seteru (
       1  38halbeta  both arguments must be gt zero, 38, 1, 2)
        if (p.ge.10.0) go to 30
        if (q.ge.10.0) go to 20
  c
  c p and q are small.
  c
        albeta = alog(gamma(p) * (gamma(q)/gamma(p+q)) )
        return
  c
  c p is small, but q is big.
  c
   20   corr = r9lgmc(q) - r9lgmc(p+q)
        albeta = alngam(p) + corr + p - p*alog(p+q) +
       1  (q-0.5)*alnrel(-p/(p+q))
        return
  c
  c p and q are big.
  c
   30   corr = r9lgmc(p) + r9lgmc(q) - r9lgmc(p+q)
        albeta = -0.5*alog(q) + sq2pil + corr + (p-0.5)*alog(p/(p+q))
       1  + q*alnrel(-p/(p+q))
        return
  c
        end