docs: function basym(a, b, x, y) {
        var z, p, q, bg, bgn, Ls, as, bs, cs, es, r, ab, sgn, abnp1,
           root2z, root2znm1, ser;
  
        p = a / (a + b);
        q = b / (a + b);
        bg = a < b ? Math.sqrt(q / a) : Math.sqrt(p / b); // betaGamma (11.3)
        bgn = 1 / bg;                        // betaGamma^n
        z = Math.sqrt(a * phi(x / p) + b * phi(y / q));   // (11.12)
        Ls = [];
        as = [];
        cs = [];
        es = [];
        ab = a <= b ? a / b : b / a;
        abnp1 = 1;  // (-1)^n * (a/b)^(n+1) or (b/a)^(n+1) for (11.4 - 11.5)
        root2z = Math.sqrt(2) * z;
        root2znm1 = 1; // (sqrt(2) * z) ^ (n - 1) for (11.21)
        sgn = -1; // (-1)^n
        ser = series(function(n) {
                 var i, seriesTerms;
  
                 bgn *= bg;
                 abnp1 = -abnp1 * ab;
                 sgn = -sgn;
                 if (n === 0) {
                    es[0] = 1;
                    Ls[0] = Math.sqrt(Math.PI) / 4 * Math.exp(z * z) *
                             erfc(z);
                 } else if (n === 1) {
                    as[1] = 2 / 3 * (1 - abnp1) * (a <= b ? q : sgn * p);
                    cs[1] = -0.5 * as[1];
                    es[1] = -cs[1];
                    Ls[1] = Math.pow(2, -1.5);
                 } else {
                    root2znm1 *= root2z;
                    Ls[n] = Math.pow(2, -1.5) * root2znm1 +
                             (n - 1) * Ls[n - 2];     // (11.21)
                    as[n] = 2 / (n + 2) * (1 - abnp1) * (a <= b ? q : sgn * p);
                    r = (-n - 1) / 2;
                    seriesTerms = function(j) {
                          if (j === 0) { return 0; }
                          return ((i - j) * r - j) * bs[j] * as[i - j];
                       };
                    bs = [1, r * as[1]]; // first two bs for current r
                    for (i = 2; i <= n; i += 1) {
                       // add the series for b_n^r
                       bs[i] = r * as[i] + 1 / i * series(seriesTerms, i);
                    }
                    cs[n] = bs[n] / (n + 1); // (11.6)
                    es[n] = -series(function(k) {
                       return es[k] * cs[n - k];
                    }, n);
                 }
                 return es[n] * (Ls[n] * bgn);
               }, 100);
        return Math.log(2 / Math.sqrt(Math.PI)) -
               bcorr(a, b) -   // (11.20)
               z * z + Math.log(ser);
     }