Lanczos approximation: Difference between revisions

Revision as of 15:48, 28 March 2017

Latest revision as of 05:12, 9 August 2024

(33 intermediate revisions by 23 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 (in pursuit of speed, for example) gives totally inaccurate results; the coefficients must be recomputed from scratch for an expansion with fewer terms.

<sourcesyntaxhighlight lang="python">

from cmath import sin, sqrt, pi, exp

def gamma(z):

epsilon = 0.0000001

def withinepsilon(x):

return abs(x) <= epsilon

from cmath import sin,sqrt,pi,exp

"""

p = [ 676.5203681218851, -1259.1392167224028, 771.32342877765313,

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

-176.61502916214059, 12.507343278686905, -0.13857109526572012,

Here are some other samples

9.9843695780195716e-6, 1.5056327351493116e-7]

z = complex(z)

g = 5

# Reflection formula (edit: this use of reflection (thus the if-else structure) seems unnecessary and just adds more code to execute. it calls itself again, so it still needs to execute the same "for" loop yet has an extra calculation at the end)

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,

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:

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

else:

z -= 1

x = p[0.99999999999980993]

for i in range(1, len(p)):

for (i, pval) in enumerate(x += p[i] / (z + i):

t = z + xg += pval/(z+i+1)0.5

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

return t = z + lendrop_imag(py) - 0.5

"""

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

The above use of the reflection (thus the if-else structure) is necessary, even though

it may look strange, as it allows to extend the approximation to values of z where

Re(z) < 0.5, where the Lanczos method is not valid.

"""

print(gamma(1))

if withinepsilon(result.imag):

print(gamma(5))

return result.real

print(gamma(0.5))

return result

</syntaxhighlight>

</source>

==See also==

Line 101 ⟶ 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

| issue = 1

| publisher = Society for Industrial and Applied Mathematics

| issn = 0887-459X

| doi= 10.23071137/29497670701008

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

}}

* {{Citation | last1=Press | first1=W. H. | last2=Teukolsky | first2=S. A. | last3=Vetterling | first3=W. T. | last4=Flannery | first4=B. P. | year=2007 | title=Numerical Recipes: The Art of Scientific Computing | edition=3rd | publisher=Cambridge University Press | publication-place___location=New York | isbn=978-0-521-88068-8 | chapter=Section 6.1. Gamma Function | chapter-url=http://apps.nrbook.com/empanel/index.html?pg=256}}

* {{cite thesis

|last=Pugh

Line 129 ⟶ 192:

[[Category:Gamma and related functions]]

[[Category:Numerical analysis]]

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