Lanczos approximation

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

\Gamma (z+1)={\sqrt {2\pi }}{\left(z+g+{\frac {1}{2}}\right)}^{{z+{\frac {1}{2}}}}e^{{-\left(z+g+{\frac {1}{2}}\right)}}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 g is a constant that may be chosen arbitrarily subject to the restriction that Re(z) > 1/2.^[1] 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,

\Gamma (1-z)\;\Gamma (z)={\pi \over \sin \pi z}.

The series A is 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 or double floating-point precision. If a fixed g is chosen, the coefficients can be calculated in advance and 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 sin x or e^x". The method is also implemented in the GNU Scientific Library.

Coefficients

The coefficients are given by

p_{k}(g)=\sum _{{a=0}}^{k}C(2k+1,2a+1){\frac {{\sqrt {2}}}{\pi }}\left(a-{\begin{matrix}{\frac {1}{2}}\end{matrix}}\right)!{\left(a+g+{\begin{matrix}{\frac {1}{2}}\end{matrix}}\right)}^{{-\left(a+{\frac {1}{2}}\right)}}e^{{a+g+{\frac {1}{2}}}}

with $C(i,j)$ denoting the (i, j)th element of the Chebyshev polynomial coefficient matrix which can be calculated recursively from the identities

$C(1,1)=1\,$
$C(2,2)=1\,$
$C(i,1)=-C(i-2,1)\,$	$i=3,4,\dots \,$
$C(i,j)=2C(i-1,j-1)\,$	$i=j=3,4,\dots \,$
$C(i,j)=2C(i-1,j-1)-C(i-2,j)\,$	$i>j=2,3,\dots .$

Paul Godfrey describes how to obtain the coefficients and also the value of the truncated series A as a matrix product.

Derivation

Lanczos derived the formula from Leonhard Euler's integral

\Gamma (z+1)=\int _{0}^{\infty }t^{{z}}\,e^{{-t}}\,dt,

performing a sequence of basic manipulations to obtain

\Gamma (z+1)=(z+g+1)^{{z+1}}e^{{-(z+g+1)}}\int _{0}^{e}[v(1-\log v)]^{{z-{\frac {1}{2}}}}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 15 correct decimal places:

def gamma(z):             # great function from Wiki, but maybe could use memoization?
    epsilon = 0.0000001
    def withinepsilon(x):
        return abs(x) <= epsilon

    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]
    z = complex(z)

    # 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)
    if z.real < 0.5:
        result = pi / (sin(pi*z) * gamma(1-z))
    else:
        z -= 1
        x = 0.99999999999980993

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

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

    if withinepsilon(result.imag):
        return result.real
    return result

References

↑ Pugh thesis

Godfrey, Paul (2001). "Lanczos Implementation of the Gamma Function".
Lanczos, Cornelius (1964). "A Precision Approximation of the Gamma Function". SIAM Journal on Numerical Analysis series B. Society for Industrial and Applied Mathematics. 1: 86–96. doi:10.2307/2949767. ISSN 0887-459X. JSTOR 2949767. External link in |journal= (help)
Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. (2007), "Section 6.1. Gamma Function", Numerical Recipes: The Art of Scientific Computing (3rd ed.), New York: Cambridge University Press, ISBN 978-0-521-88068-8
Pugh, Glendon (2004). An analysis of the Lanczos Gamma approximation (PDF) (PhD thesis).
Toth, Viktor (2005). "Programmable Calculators: The Lanczos Approximation".
Weisstein, Eric W. "Lanczos Approximation". MathWorld.

This article is issued from Wikipedia - version of the 4/1/2016. The text is available under the Creative Commons Attribution/Share Alike but additional terms may apply for the media files.