code: function bcorr(a, b) {
        if (a > b) { return bcorr(b, a); }
        // a <= b
        return stirlerr(a) + algdiv(a, b) +
                 (a + b - 0.5) * Math.log(1 + a / b) +
                 a * (Math.log(b) - 1);
     }
docs: bcorr Computes stirlerr(a) + stirlerr(b) - stirlerr(a + b) see DiDonato (23)