Gamma function with the Lanczos approximation (Python)

From LiteratePrograms
Jump to: navigation, search
Other implementations: Perl | Python

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]
        for i in range(1, g+2):
            x += lanczos_coef[i]/(z+i)
        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
hijacker
hijacker
hijacker
hijacker