Nemes.scala

  def gamma(z: Real): Real = {
    // The Nemes approximation to the log Gamma function isn't very accurate for 0 < z < 0.5,
    // so we use a trick taken from Boost's lgamma function and calculate Log(Gamma(z + 1)) - Log(z) instead.
    // See https://www.boost.org/doc/libs/1_50_0/libs/math/doc/sf_and_dist/html/math_toolkit/special/sf_gamma/lgamma.html
    nemesGamma(z + 1) - z.log
  }

  private def nemesGamma(z: Real): Real = {
    // An approximation to Log(Gamma(z)) due to Gergő Nemes. See https://en.wikipedia.org/wiki/Stirling%27s_approximation
    val w = z + (Real.one / ((12 * z) - (Real.one / (10 * z))))
    (Real(Math.PI * 2).log / 2) - (z.log / 2) + (z * (w.log - 1))
  }