gamma function

double gamma(
  1. double z
)

Computes 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 should achieve 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).

Our unit tests suggest that the accuracy of the Gamma function is correct up to 13 floating point digits.

Implementation

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

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

    return s *
        constants.twoSqrtEOverPi *
        stdmath.pow((z - 0.5 + _gammaR) / constants.e, z - 0.5);
  }
}