Lanczos approximation

Template:Short description In mathematics, the Lanczos approximation is a method for computing the gamma function numerically, published by Cornelius Lanczos in 1964. It is a practical alternative to the more popular Stirling's approximation for calculating the gamma function with fixed precision.

Introduction

The Lanczos approximation consists of the formula

$\mathrm{\Gamma }(z+1)=\sqrt{2\pi }{(z+g+\frac{1}{2})}^{z+\frac{1}{2}}{e}^{-(z+g+\frac{1}{2})}{A}_g(z)$

for the gamma function, with

$A_g(z)=\frac{1}{2}{p}_0(g)+p_1(g)\frac{z}{z+1}+p_2(g)\frac{z(z-1)}{(z+1)(z+2)}+\cdots .$

Here Template:Mvar is a real constant that may be chosen arbitrarily subject to the restriction that Template:Math. The coefficients Template:Mvar, which depend on Template:Mvar, are slightly more difficult to calculate (see below). Although the formula as stated here is only valid for arguments in the right complex half-plane, it can be extended to the entire complex plane by the reflection formula,

$\mathrm{\Gamma }(1-z)\mathrm{\Gamma }(z)=\frac{\pi }{\sin \pi z}.$

The series Template:Mvar is convergent, and may be truncated to obtain an approximation with the desired precision. By choosing an appropriate Template:Mvar (typically a small integer), only some 5–10 terms of the series are needed to compute the gamma function with typical single or double floating-point precision. If a fixed Template:Mvar is chosen, the coefficients can be calculated in advance and, thanks to partial fraction decomposition, the sum is recast into the following form:

$A_g(z)=c_0+\sum _{k=1}^N\frac{c_k}{z+k}$

Thus computing the gamma function becomes a matter of evaluating only a small number of elementary functions and multiplying by stored constants. The Lanczos approximation was popularized by Numerical Recipes, according to which computing the gamma function becomes "not much more difficult than other built-in functions that we take for granted, such as Template:Math or Template:Math." The method is also implemented in the GNU Scientific Library, Boost, CPython and musl.

Coefficients

The coefficients are given by

$p_k(g)=\frac{\sqrt{2}}{\pi }\sum _{l=0}^k{C}_{2k+1,2l+1}(l-\frac{1}{2})!{(l+g+\frac{1}{2})}^{-(l+\frac{1}{2})}{e}^{l+g+\frac{1}{2}}$

where Template:Math represents the Template:Mathth element of the matrix of coefficients for the Chebyshev polynomials, which can be calculated recursively from these identities:Template:Cn

$\begin{array}{cc}{C}_{1,1}& =1\\ C_{2,2}& =1\\ C_{n+1,1}& =-C_{n-1,1}& \text{ for }n& =2,3,4\dots \\ C_{n+1,n+1}& =2C_{n,n}& \text{ for }n& =2,3,4\dots \\ C_{n+1,m+1}& =2C_{n,m}-C_{n-1,m+1}& \text{ for }n& >m=1,2,3\dots \end{array}$

Godfrey (2001) describes how to obtain the coefficients and also the value of the truncated series Template:Mvar as a matrix product with good precision in floating-point.

Derivation

Lanczos derived the formula from Leonhard Euler's integral

$\mathrm{\Gamma }(z+1)={\int }_0^{\mathrm{\infty }}{t}^z{e}^{-t}dt,$

performing a sequence of basic manipulations to obtain

$\mathrm{\Gamma }(z+1)={(z+g+1)}^{z+1}{e}^{-(z+g+1)}{\int }_0^e(v(1-\log v){)}^z{v}^{g}dv,$

and deriving a series for the integral.

Simple implementation

The following implementation in the Python programming language works for complex arguments and typically gives 13 correct decimal places. Note that omitting the smallest coefficients (in pursuit of speed, for example) gives totally inaccurate results; the coefficients must be recomputed from scratch for an expansion with fewer terms.

from cmath import sin, sqrt, pi, exp

# Constant g and Coefficients p_k(g)
G = 7
P = [
    0.99999999999980993,
    676.5203681218851,
    -1259.1392167224028,
    771.32342877765313,
    -176.61502916214059,
    12.507343278686905,
    -0.13857109526572012,
    9.9843695780195716e-6,
    1.5056327351493116e-7
]

EPSILON = 1e-07
def drop_imag(z: complex) -> float | complex:
    """Turn (a+bj) into a when |b| <= EPSILON"""
    if abs(z.imag) <= EPSILON:
        z = z.real
    return z


def gamma(z: int | float | complex) -> float | complex:
    z = complex(z)
    if z.real < 0.5:
        # Lanczos is nor valid for this z.real. Use reflection formula.
        return pi / (sin(pi * z) * gamma(1 - z))  
    else:
        z -= 1
        x = P[0] + sum(P[i] / (z + i) for i in range(1, len(P)))
        t = z + G + 0.5
        return sqrt(2 * pi) * t ** (z + 0.5) * exp(-t) * x

print(gamma(1))
print(gamma(5))   
print(gamma(0.5))

"""
The coefficients used in the code are for when g = 7 and n = len(p) = 9
For calculation of p_k, see Godfrey (2001)

Here are some other samples

g = 5
n = 5
p = [
    1.0000018972739440364,
    76.180082222642137322,
    -86.505092037054859197,
    24.012898581922685900,
    -1.2296028490285820771
]

g = 5
n = 7
p = [
    1.0000000001900148240,
    76.180091729471463483,
    -86.505320329416767652,
    24.014098240830910490,
    -1.2317395724501553875,
    0.0012086509738661785061,
    -5.3952393849531283785e-6
]

g = 8
n = 12
p = [
    0.9999999999999999298,
    1975.3739023578852322,
    -4397.3823927922428918,
    3462.6328459862717019,
    -1156.9851431631167820,
    154.53815050252775060,
    -6.2536716123689161798,
    0.034642762454736807441,
    -7.4776171974442977377e-7,
    6.3041253821852264261e-8,
    -2.7405717035683877489e-8,
    4.0486948817567609101e-9
]
"""

See also

• Stirling's approximation
• Spouge's approximation

References

• Template:Cite web
• Template:Cite journal
• Template:Citation
• Template:Cite thesis
• Template:Cite web
• Template:MathWorld