Lanczos approximation: Difference between revisions

Content deleted Content added
Introduction:restriction on g: it was correct as it was - see cited ref sect 5.3 page 69 - in any case, a restriction on g would have to involve g
Citation bot (talk | contribs)
Added issue. | Use this bot. Report bugs. | Suggested by Abductive | Category:Gamma and related functions | #UCB_Category 18/38
 
(9 intermediate revisions by 6 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 14 ⟶ 15:
:<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>
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 53 ⟶ 54:
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 75 ⟶ 125:
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