gammaLn function

double gammaLn (double z)

Computes the logarithm of the Gamma function.

This implementation of the computation of the gamma and logarithm of the gamma function follows the derivation in "An Analysis Of The Lanczos Gamma Approximation", Glendon Ralph Pugh, 2004. We use the implementation listed on p. 116 which achieves an accuracy of 16 floating point digits. Although 16 digit accuracy should be sufficient for double values, improving accuracy is possible (see p. 126 in Pugh).

Unit tests suggest that the accuracy of the Gamma function is correct up to 14 floating point digits.

Implementation

double gammaLn(double z) {
  if (z < 0.5) {
    double s = _gammaDk[0];
    for (int i = 1; i <= _gammaN; i++) {
      s += _gammaDk[i] / (i - z);
    }

    return constants.lnPi -
        stdmath.log(trig.sin(constants.pi * z)) -
        stdmath.log(s) -
        constants.logTwoSqrtEOverPi -
        ((0.5 - z) * stdmath.log((0.5 - z + _gammaR) / constants.e));
  } else {
    double s = _gammaDk[0];
    for (int i = 1; i <= _gammaN; i++) {
      s += _gammaDk[i] / (z + i - 1.0);
    }

    return stdmath.log(s) +
        constants.logTwoSqrtEOverPi +
        ((z - 0.5) * stdmath.log((z - 0.5 + _gammaR) / constants.e));
  }
}