ランチョス近似

数学において、ランチョス近似は、1964 年にコーネリアス ランチョスによって発表されたガンマ関数を数値的に計算する方法です。これは、固定精度でガンマ関数を計算するための、より一般的なスターリング近似の実用的な代替手段です。

ランチョス近似は次の式で表される。

${\displaystyle \Gamma (z+1)={\sqrt {2\pi }}{\left(z+g+{\tfrac {1}{2}}\right)}^{z+1/2}e ^{-(z+g+1/2)}A_{g}(z)}$

ガンマ関数の場合、

${\displaystyle 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 .}$

ここでgは実定数であり、Re( z + g + ⁠1/2⁠ ) > 0. ^[1] gに依存する係数p の計算はやや難しい（下記参照）。ここで述べた式は右複素半平面の引数に対してのみ有効であるが、反射公式によって複素平面全体に拡張することができる。

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

級数Aは収束し、必要な精度の近似値を得るために切り捨てられる場合があります。適切なg (通常は小さな整数) を選択すると、通常の単精度または倍精度 浮動小数点精度でガンマ関数を計算するのに必要な級数の項は 5 ～ 10 個程度になります。固定のgを選択した場合は、係数を事前に計算することができ、部分分数分解により、合計は次の形式に書き直されます。

${\displaystyle A_{g}(z)=c_{0}+\sum _{k=1}^{N}{\frac {c_{k}}{z+k}}}$

したがって、ガンマ関数の計算は、少数の基本関数を評価し、格納されている定数を掛け合わせるだけの問題になります。Lanczos 近似はNumerical Recipesによって普及しました。これによると、ガンマ関数の計算は「 sin xやe ^xなど、私たちが当然と思っている他の組み込み関数よりもそれほど難しくありません 」。この方法は、GNU Scientific Library、Boost、CPython、muslにも実装されています。

係数は次のように与えられる。

${\displaystyle 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}}$

ここで、はチェビシェフ多項式の係数行列の( n , m )番目の要素を表し、次の恒等式から再帰的に計算できます。 ${\displaystyle C_{n,m}}$

${\displaystyle {\begin{aligned}C_{1,\,1}&=1\\[5px]C_{2,\,2}&=1\\[5px]C_{n+1,\,1}&=-\,C_{n-1,\,1}&{\text{ for }}n&=2,3,4\,\dots \\[5px]C_{n+1,\,n+1}&=2\,C_{n,\,n}&{\text{ for }}n&=2,3,4\,\dots \\[5px]C_{n+1,\,m+1}&=2\,C_{n,\,m}-C_{n-1,\,m+1}&{\text{ for }}n&>m=1,2,3\,\dots \end{aligned}}}$

Godfrey（2001）は係数と切り捨てられた級数Aの値を行列積として得る方法を説明しています。^[2]

ランチョスはレオンハルト・オイラーの積分からこの式を導いた。

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

一連の基本的な操作を実行して

${\displaystyle \Gamma (z+1)=(z+g+1)^{z+1}e^{-(z+g+1)}\int _{0}^{e}{\Big (}v(1-\log v){\Big )}^{z}v^{g}\,dv,}$

そして積分の級数を導出します。

Python プログラミング言語での次の実装は、複素引数に対して機能し、通常は 13 桁の正確な小数点以下桁数を返します。最小の係数を省略すると (たとえば、速度を追求するため)、まったく不正確な結果が得られることに注意してください。項の数が少ない展開では、係数を最初から再計算する必要があります。

 cmath からsin 、sqrt 、pi 、expをインポートします    



コードで使用されている係数は、g = 7、n = 9の場合のもの
です。他のサンプルもいくつかあります。

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.999999999999999298, 
    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.5203​​681218851 、
    - 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 の場合、 
    zを返します。 

def  gamma ( z ): 
    z  =  complex ( z ) 
    if  z . real  <  0.5 : 
        y  =  pi  /  ( sin ( pi  *  z )  *  gamma ( 1  -  z ))   # 反射式
    else : 
        z  -=  1 
        x  =  p [ 0 ] 
        for  i  in  range ( 1 ,  len ( p ) ): 
            x  +=  p [ i ]  /  ( z  +  i ) 
        t  =  z  +  g  +  0.5 
        y  =  sqrt ( 2  *  pi )  *  t  **  ( z  +  0.5 )  *  exp ( - t )  *  x 
    return  drop_imag ( y ) 
"""
上記の反射の使用 (したがって、if-else 構造) は、奇妙に見えるかもしれませんが、
ランチョス法が有効ではない Re(z) < 0.5の z 値に近似を拡張できるため、必要です。 
"""


印刷(ガンマ( 1 ))
印刷(ガンマ( 5 ))   
印刷(ガンマ( 0.5 ))

• スターリング近似
• スポージの近似

• Godfrey, Paul (2001)。「Lanczos ガンマ関数の実装」
• Lanczos, Cornelius (1964). 「ガンマ関数の高精度近似」. Journal of the Society for Industrial and Applied Mathematics, Series B: Numerical Analysis . 1 (1): 86–96. Bibcode :1964SJNA....1...86L. doi :10.1137/0701008. ISSN  0887-459X. JSTOR  2949767.
• Press, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP (2007)、「セクション 6.1. ガンマ関数」、Numerical Recipes: The Art of Scientific Computing (第 3 版)、ニューヨーク: Cambridge University Press、ISBN 978-0-521-88068-8
• Pugh, Glendon (2004). Lanczos ガンマ近似の分析(PDF) (博士論文).
• Toth, Viktor (2005)。「プログラム可能な計算機: Lanczos 近似」。
• ワイスタイン、エリック・W.「ランチョスの近似」。マスワールド。