Lanczos approximation: Difference between revisions

Revision as of 18:01, 23 April 2014

Latest revision as of 05:12, 9 August 2024

(48 intermediate revisions by 34 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 [[Gammagamma function]] numerically, published by [[Cornelius Lanczos]] in 1964. It is a practical alternative to the more popular [[Stirling's approximation]] for calculating the Gammagamma function with fixed precision.

==Introduction==

The Lanczos approximation consists of the formula

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

for the Gammagamma function, with

:<math>A_g(z) = \frac{1}{2}p_0frac12p_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'') > +''g''+{{sfrac|1/|2}})&nbsp;>&nbsp;0.<ref>Pugh{{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-105–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 Gammagamma 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\,2\ell+1)} \frac{left(\sqrtell - \tfrac{2}1}{\pi2} \right)! {\left(a\ell -+ \begin{matrix}g + \fractfrac{1}{2} \endright)}^{matrix-(\ell+1/2)} e^{\right)!ell + g + 1/2 }</math>

{\left(a + g + \begin{matrix} \frac{1}{2} \end{matrix} \right)}^{- \left( a + \frac{1}{2} \right) } e^{a + g + \frac{1}{2} }</math>

withwhere <math>C(iC_{n,j)m}</math> 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}

:{|

|<math>C( C_{1,\,1)} &= 1\,</math> ||\\[5px]

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

|-

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

|<math>C(2,2) = 1\,</math> ||

C_{n+1,\,n+1} &= 2\,C_{n,\,n} & \text{ for } n &= 2, 3, 4\, \dots \\[5px]

|-

|<math>C(iC_{n+1,\,m+1)} &= -C(i-2\,C_{n,\,m} - C_{n-1),\,</math>m+1} ||& \text{ for } n & <math>i m = 31, 42, \dots3\,</math> \dots

\end{align}</math>

|-

|<math>C(i,j) = 2 C(i-1, j-1)\,</math> || <math>i = j = 3, 4, \dots\,</math>

|-

|<math>C(i,j) = 2 C(i-1, j-1) - C(i-2, j)\,</math> || <math>i > j = 2, 3, \dots .</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 ⟶ 43:

performing a sequence of basic manipulations to obtain

:<math>\Gamma(z+1) = (z+g+1)^{z+1} e^{-(z+g+1)} \int_0^e [\Big(v(1-\log v)]\Big)^{z-\frac{1}{2}} v^g\,dv,</math>

and deriving a series for the integral.

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

"""

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

]

"""

# Coefficients used by the GNU Scientific Library

g = 7

n = 9

p = [0.99999999999980993, 676.5203681218851, -1259.1392167224028,

p = [

771.32342877765313, -176.61502916214059, 12.507343278686905,

0.99999999999980993,

-0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7]

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)

# Reflection formula

if z.real < 0.5:

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

else:

z -= 1

x = p[0]

for i in range(1, g+2len(p)):

x += p[i] / (z + i)

t = z + g + 0.5

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

return drop_imag(y)

</source>

"""

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

print(gamma(5))

print(gamma(0.5))

</syntaxhighlight>

==See also==

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

| id = ISSN: 0887459X

| 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 125 ⟶ 192:

[[Category:Gamma and related functions]]

[[Category:Numerical analysis]]

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