docs: function qnbinom(size, prob, lowerTail, logp) {
        var mu, sigma, gamma;
  
        logp = logp === true;
        lowerTail = lowerTail !== false;
  
        mu = size * (1 - prob) / prob;
        sigma = Math.sqrt(size * (1 - prob) / prob / prob);
        gamma = (2 - prob) / prob / sigma;
  
        if (utils.hasNaN(prob, size)) {
           return function() { return NaN; };
        }
        if (prob === 0 && size === 0) {
           return function(p) { return utils.hasNaN(p) ? NaN : 0; };
        }
        if (prob <= 0 || prob > 1 || size < 0) {
           return function() { return NaN; };
        }
        if (prob === 1 || size === 0) {
           return function(p) { return utils.hasNaN(p) ? NaN : 0; };
        }
  
        /* eslint-disable complexity */
        return utils.pWrap(lowerTail, logp, function(ps) {
           var z, ret;
  
           if (ps.p === 0) { return 0; }
           if (ps.q === 0) { return Infinity; }
           // Our addition, R loops instead
           if (size === Infinity) { return Infinity; }
  
           z = ps.p > 0.9 ? qnorm(0, 1, false)(ps.q) : qnorm(0, 1)(ps.p); // initial value
           ret = Math.floor(mu + sigma * (z + gamma * (z * z - 1) / 6) + 0.5);
  
           // Supposed to do: ps.p *= 1 - 64*DBL_EPSILON here
           if (ps.p > 0.9) {
              return utils.discInvCdf2(0, Infinity, ret, ps.q,
                        pnbinom(size, prob, false));
           }
           return utils.discInvCdf(0, Infinity, ret, ps.p, pnbinom(size, prob));
        });
        /* eslint-enable complexity */
     }