pub fn ln_beta(a: f64, b: f64) -> f64
ln B(a, b) = ln Gamma(a) + ln Gamma(b) - ln Gamma(a + b), same domain rule per argument.