docs: function lbeta(a, b) {
        if (a > b) { return lbeta(b, a); }
        // a <= b
        if (a < 0) { return NaN; }
        if (a === 0) { return Infinity; }
        if (b === Infinity) { return -Infinity; }
        if (a > 10) {
           return -0.5 * Math.log(b) + Math.log(sqrt2pi) +
                  (stirlerr(a) + stirlerr(b) - stirlerr(a + b)) +
                  (a - 0.5) * Math.log(a / (a + b)) +
                  b * log1p(-a / (a + b));
        }
        if (b >= 10) {
           return lgamma(a) + (stirlerr(b) - stirlerr(a + b)) +
                  a - a * Math.log(a + b) +
                  (b - 0.5) * log1p(-a / (a + b));
        }
        return lgamma(a) + (lgamma(b) - lgamma(a + b));
     }