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