Lanczos approximation: Difference between revisions

Revision as of 18:03, 11 May 2020

Latest revision as of 05:12, 9 August 2024

(17 intermediate revisions by 13 users not shown)

Line 1:

{{Short description|Numerical method for calculating the gamma function}}

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.

Line 10 ⟶ 11:

:<math>A_g(z) = \frac12p_0(g) + p_1(g) \frac{z}{z+1} + p_2(g) \frac{z(z-1)}{(z+1)(z+2)} + \cdots.</math>

Here ''g'' is a real [[Constant (mathematics)|constant]] that may be chosen arbitrarily subject to the restriction that Re(''z'')&nbsp;>&nbsp;+''g''+{{sfrac|1|2}})&nbsp;>&nbsp;0.<ref>{{cite thesis |url=https://web.viu.ca/pughg/phdThesis/phdThesis.pdf#110 |author=Pugh, Glendon |year=2004 |title=An analysis of the Lanczos Gamma approximation |type=Ph.D.}}</ref> The coefficients ''p'', which depend on ''g'', 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]],

:<math>\Gamma(1-z) \; \Gamma(z) = {\pi \over \sin \pi z}.</math>

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

:<math>A_g(z) = c_0 + \sum_{k=1}^{N} \frac{c_k}{z+k}</math>

Thus computing the gamma function becomes a matter of evaluating only a small number of [[elementary function]]s 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 sin&nbsp;''x'' or ''e''^''x''"." The method is also implemented in the [[GNU Scientific Library]], [[Boost (C++ libraries)|Boost]], [[CPython]] and [[musl]].

==Coefficients==

Line 24 ⟶ 25:

:<math>p_k(g) = \frac{\sqrt{2\,}}{\pi} \sum_{\ell=0}^k C_{2k+1,\,2\ell+1} \left(\ell - \tfrac{1}{2} \right)! {\left(\ell + g + \tfrac{1}{2} \right)}^{-(\ell+1/2)} e^{\ell + g + 1/2 }</math>

where <math>C_{n,m}</math> represents the (''n'', ''m'')th element of the [[matrix (mathematics)|matrix]] of coefficients for the [[Chebyshev polynomialpolynomials]]s, which can be calculated [[recursion|recursively]] from these identities:

:<math>\begin{align}

Line 47 ⟶ 48:

==Simple implementation==

The following implementation in the [[Python (programming language)|Python programming language]] works for complex arguments and typically gives 1513 correct decimal places.

Note that omitting the smallest coefficients does(in notpursuit result inof aspeed, fasterfor butexample) slightlygives lesstotally accurateinaccurate implementationresults; the coefficients must be recomputed from scratch for an expansion with fewer terms.

<syntaxhighlight lang="python">

from cmath import sin, sqrt, pi, exp

p = [676.5203681218851

"""

,-1259.1392167224028

The coefficients used in the code are for when g = 7 and n = 9

,771.32342877765313

Here are some other samples

,-176.61502916214059

,12.507343278686905

g = 5

,-0.13857109526572012

n = 5

,9.9843695780195716e-6

p = [

,1.5056327351493116e-7

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

]

"""

g = 7

n = 9

p = [

0.99999999999980993,

p = [ 676.5203681218851,

,-1259.1392167224028,

,771.32342877765313,

,-176.61502916214059,

,12.507343278686905,

,-0.13857109526572012,

,9.9843695780195716e-6,

,1.5056327351493116e-7

]

EPSILON = 1e-07

Line 72 ⟶ 122:

z = complex(z)

if z.real < 0.5:

y = pi / (sin(pi * z) * gamma(1 - z)) # Reflection formula

else:

z -= 1

x = p[0.99999999999980993]

for (i, pval) in enumeraterange(1, len(p)):

x += pvalp[i] / (z + i+1)

t = z + len(p)g -+ 0.5

y = sqrt(2 * pi) * t ** (z + 0.5) * exp(-t) * x

return drop_imag(y)

"""

Line 90 ⟶ 140:

print(gamma(5))

print(gamma(0.5))

</syntaxhighlight>

Line 116 ⟶ 167:

| pages = 86&ndash;96

| year = 1964

| issnissue = 0887-459X1

| issn = 0887-459X

| doi= 10.1137/0701008

| bibcode = 1964SJNA....1...86L

Line 140 ⟶ 192:

[[Category:Gamma and related functions]]

[[Category:Numerical analysis]]

[[Category:Articles with example Python (programming language) code]]