Lanczos approximation: Difference between revisions

Revision as of 19:18, 26 January 2019

Latest revision as of 05:12, 9 August 2024

(26 intermediate revisions by 18 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 4 ⟶ 5:

The Lanczos approximation consists of the formula

:<math>\Gamma(z+1) = \sqrt{2\pi} {\left( z + g + \tfrac12 \right)}^{z + \frac121/2 } e^{-\left(z+g+\frac12\right1/2)} A_g(z)</math>

for the gamma function, with

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, thesis]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 Gammagamma 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 Gammagamma 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==

The coefficients are given by

:<math>p_k(g) = \frac{\sqrt{2\,}}{\pi} \sum_{a\ell=0}^k C(C_{2k+1, 2a+1) \frac{\sqrt{,2}}{\piell+1} \left(a\ell - \tfrac12tfrac{1}{2} \right)! {\left(a\ell + g + \tfrac12tfrac{1}{2} \right)}^{- \left( a \ell+ \frac12 \right1/2) } e^{a\ell + g + \frac121/2 }</math>

withwhere ''C''(''i''<math>C_{n,m}</math> ''j'') denotingrepresents the (''in'', ''jm'')th element of the [[Chebyshevmatrix polynomial(mathematics)|matrix]] coefficientof coefficients for the [[matrixChebyshev (mathematics)|matrixpolynomials]], which can be calculated [[recursion|recursively]] from thethese identities:

:<math>\begin{align}

C( C_{1,\,1)} &= 1 \\[5px]

C( C_{2,\,2)} &= 1 \\[5px]

C(i C_{n+1,\,1)} &= -C(i\,C_{n-21,\, 1)} & i\text{ for } n &= 2, 3, 4\, \dots \\[5px]

C(iC_{n+1,j)\,n+1} &= 2 C(i-1\,C_{n,\,n} j-1) & i &= j\text{ for } n &= 2, 3, 4\, \dots \\[5px]

C(iC_{n+1,j)\,m+1} &= 2\,C_{n,\,m} C(i-1, jC_{n-1),\,m+1} -& C(i-2,\text{ j)for & i} n & > jm = 1, 2, 3\, \dots

\end{align}</math>

PaulGodfrey Godfrey(2001) describes how to obtain the coefficients and also the value of the truncated series ''A'' as a [[matrix multiplication|matrix product]].<ref name=Godfrey2001>{{cite web |last=Godfrey |first=Paul |year=2001 |title=Lanczos implementation of the gamma function |url=http://www.numericana.com/answer/info/godfrey.htm |website=Numericana}}</ref>

==Derivation==

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.

<sourcesyntaxhighlight lang="python">

from cmath import sin, sqrt, pi, exp

p = [676.5203681218851

,-1259.1392167224028

,771.32342877765313

,-176.61502916214059

,12.507343278686905

,-0.13857109526572012

,9.9843695780195716e-6

,1.5056327351493116e-7

]

"""

EPSILON = 1e-07

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

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

]

"""

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

def drop_imag(z):

if abs(z.imag) <= EPSILON:

z = z.real

return z

def gamma(z):

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 87 ⟶ 137:

"""

print (gamma(1))

print (gamma(5) )

print (gamma(0.5))

</source>

</syntaxhighlight>

==See also==

Line 112 ⟶ 163:

| title = A Precision Approximation of the Gamma Function

| jstor = 2949767

| journal =Journal of the Society for Industrial and Applied Mathematics, Series B: Numerical Analysis

| journal = [http://www.siam.org/journals/sinum.php SIAM Journal on Numerical Analysis series B]

| volume = 1

| 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]]