Gamma function with the Lanczos approximation (Python)

From LiteratePrograms

Other implementations: Perl | Python | Ruby

The Lanczos approximation offers a relatively simple way to calculate the gamma function of a complex argument to any desired precision. The precision is determined by a free constant g and a series whose coefficients depend on g and may be computed in advanced if a fixed precision is desired. The following implementation uses the same coefficients as the GNU Scientific Library's gamma function. The typical relative error is about 10^-15 (nearly full double precision accuracy).

Lanczos' series is only valid for numbers in the right half complex plane, so the reflection formula for the gamma function is employed for negative arguments (and also arguments in (0, 0.5), since this should increase accuracy near the singularity at zero.)

<<gamma_lanczos.py>>=
from cmath import *
g = 7
lanczos_coef = [ \
     0.99999999999980993,
   676.5203681218851,
 -1259.1392167224028,
   771.32342877765313,
  -176.61502916214059,
    12.507343278686905,
    -0.13857109526572012,
     9.9843695780195716e-6,
     1.5056327351493116e-7]
def gamma(z):
    z = complex(z)
    if z.real < 0.5:
        return pi / (sin(pi*z)*gamma(1-z))
    else:
        z -= 1
        x = lanczos_coef[0] + \
            sum(lanczos_coef[i]/(z+i)
                for i in range(1, g+2))
        t = z + g + 0.5
        return sqrt(2*pi) * t**(z+0.5) * exp(-t) * x

Testing against a known value (Γ(10) = 9! = 362880):

>>> gamma(10.0)
(362880.00000000047+0j)

Download code