diGammaInv function

double diGammaInv (double p)

Computes the inverse Digamma function: this is the inverse of the logarithm of the gamma function. This function will only return solutions that are positive.

This implementation is based on the bisection method.

Implementation

double diGammaInv(double p) {
  if (p.isNaN) {
    return double.nan;
  }

  if (p.isInfinite && p.isNegative) {
    return 0.0;
  }

  if (p.isInfinite && !p.isNegative) {
    return double.infinity;
  }

  var x = stdmath.exp(p);
  for (var d = 1.0; d > 1.0e-15; d /= 2.0) {
    x += d * (p - diGamma(x)).sign;
  }

  return x;
}