docs: function bratioLog(a, b, x) {
        var p, q, y, bint, bbar, lambda;
  
        p = a / (a + b);
        q = b / (a + b);
        y = 1 - x;
        bint = Math.floor(b);
        bbar = b - bint;
        if (bbar === 0) {
           bbar += 1;
           bint -= 1;
        }
        lambda = a > b ? (a + b) * y - b
                       : a - (a + b) * x;
  
        // By convention, if 0 < x < 1 then a == 0 -> result = 1
        // and b == 0 -> result = 0. Can't have a == b == 0.
        if (a === 0) { return { lower: 0, upper: -Infinity }; }
        if (b === 0) { return { lower: -Infinity, upper: 0 }; }
  
        // R code handles specially the case where max(a, b) < eps * .001
        // but we will not (for now, anyway)
  
        // This logic follows pages 368-369 of algorithm 708
        /*  CASE where min(a, b) <= 1 */
        if (Math.min(a, b) <= 1) {
           if (x > 0.5) { return tailFlip(bratioLog(b, a, y)); }
           // we are skipping the case, in R code, for b < min(eps, a*eps)
           // also skipping a < min(eps, b*eps)
           if (Math.max(a, b) > 1) {
              // either a <= 1 < b or b <= 1 < a
              if (b <= 1) { return fromLower(bpser(a, b, x)); } // (12c)
              // a <= 1 < b
              if (x >= 0.3) { return fromUpper(bpser(b, a, y)); } // (13b)
              // x < 0.3
              if (x < 0.1 && Math.pow(x * b, a) <= 0.7) { return fromLower(bpser(a, b, x)); } // (12d)
              // x >= 0.1 or (bx)^a > 0.7
              // for now, assume bgrat will start with w0 = 0
              if (b > 15) { return fromUpper(bgrat(b, a, y)); } // (14a, 14b)
              // b <= 15
              return fromUpper(logspaceAdd(bup(b, a, y, x, 20),
                                           bgrat(b + 20, a, y))); // (15a, 15b)
           }
           // max(a, b) <= 1
           if (a >= Math.min(0.2, b)) { return fromLower(bpser(a, b, x)); } // (12a)
           // a < min(0.2, b)
           if (Math.pow(x, a) <= 0.9) { return fromLower(bpser(a, b, x)); } // (12b)
           // x^a > 0.9
           if (x >= 0.3) { return fromUpper(bpser(b, a, y)); } // (13a)
           // x < 0.3
           return fromUpper(logspaceAdd(bup(b, a, y, x, 20),
                                           bgrat(b + 20, a, y)));  // 15c
        }
        // min(a, b) > 1
        if (x > p) { return tailFlip(bratioLog(b, a, y)); }
        if (b < 40) {
           if (b * x <= 0.7 || lambda > 650) { return fromLower(bpser(a, b, x)); } // (16a)
           // b * x > 0.7
           if (x <= 0.7) {
              return fromLower(logspaceAdd(bup(bbar, a, y, x, bint),
                                           bpser(a, bbar, x)));
           } // (17a)
           // x > 0.7
           if (a > 15) {
              return fromLower(logspaceAdd(bup(bbar, a, y, x, bint),
                                           bgrat(a, bbar, x)));
           } // (18a)
           // a <= 15
           return fromLower(logspaceAdd(
                                logspaceAdd(bup(bbar, a, y, x, bint),
                                            bup(a, bbar, x, y, 20)),
                                bgrat(a + 20, bbar, x))); // (19a)
        }
        // b >= 40
        if (a <= b) {
           if (a <= 100) { return fromLower(bfrac(a, b, x, y)); } // (20a)
           // a > 100
           if (x < 0.97 * p) { return fromLower(bfrac(a, b, x, y)); } // (20b)
           // x >= 0.97 * p
           return fromLower(basym(a, b, x, y)); // (21a)
        }
        // a > b
        if (b <= 100) { return fromLower(bfrac(a, b, x, y)); } // (20c)
        // b > 100
        if (y > 1.03 * q) { return fromLower(bfrac(a, b, x, y)); } // (20d)
        // y <= 1.03 * q
        return fromLower(basym(a, b, x, y)); // (21c)
     }