利用者:Flightbridge/sandbox/ランチョス近似

（en:Lanczos approximation oldid=713051180）

ランチョス近似（英: Lanczos approximation）とは、1964年にコルネリウス・ランチョスにより発表された^[1]、ガンマ関数を数値計算する方法である。これは、ガンマ関数を一定の精度で計算する方法として有名なスターリングの公式の実用的な代替となっている。

導入

ランチョス近似は、次のガンマ関数についての式からなる。

\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)\qquad ({\mathrm {Re} }(z)>0)

ここで

A_{g}(z)={\frac {1}{2}}p_{0}+{\frac {z}{z+1}}p_{1}+{\frac {z(z-1)}{(z+1)(z+2)}}p_{2}+\cdots

である。 $g \geq 0$ は任意に選ばれる定数^[2]であり、係数 $p 0, p 1, ...$ は $g$ によって決まる値である（求め方は後述）。上記の式は複素右半平面（英語版）でしか成り立たないが、相補公式（英語版）を用いることで複素平面全体へと拡張できる。

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

級数 $A g$ は収束するので、近似値を求めるには必要な精度に応じて計算を打ち切ればよい。一般的な単精度あるいは倍精度浮動小数点を用いる場合、適切な $g$ （一般に小さな整数）を選べば 5 から 10 項ほどでガンマ関数を計算できる。級数 $A g$ の第 $N + 1$ 項での打ち切りを考えると、打ち切り誤差 $ε$ は $g = 5, N = 6$ のときで $| ε | < 2 \times 10 -10$ よりも小さくなる^[3]。

$N = 2$ のとき、各項を部分分数分解して整理すると

{\begin{aligned}&{\frac {1}{2}}p_{0}+{\frac {z}{z+1}}p_{1}+{\frac {z(z-1)}{(z+1)(z+2)}}p_{2}\\&={\biggl (}{\frac {1}{2}}p_{0}+p_{1}+p_{2}{\biggr )}+{\frac {1}{z+1}}{\biggl (}-p_{1}+2p_{2}{\biggr )}+{\frac {1}{z+2}}{\biggl (}-6p_{2}{\biggr )}\end{aligned}}

となる。同様に、第 $N + 1$ 項で打ち切った $A g$ は適当な $c 0, ..., c N$ を用いて

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

と書ける。このとき各 $c 0, ..., c N$ は $p 0, ..., p N$ の式であるから、 $g$ を固定すれば定数として予め計算しておくことが可能になる。

以上よりガンマ関数の計算は、ごく少数の初等関数の評価と予め求めておいた定数との積へと帰着される。ランチョス近似は『ニューメリカルレシピ（英語版）』において「こうすれば，ガンマ関数の計算は，当然初めから用意されている $sin(x)$ や $e x$ のような関数の計算と比べて，ことさら難しい計算になることはない^[3]」として紹介され、知られるようになった。この方法はGNU Scientific Libraryにおいて実装されている^[4]。

係数

係数 $p 0, p 1, ...$ は次の式で与えられる。

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

ここで $C (i, j)$ は、 $i$ 次のチェビシェフ多項式における $x j$ の係数を表す^[2]。この値は次の式を用いて再帰的に求められる。

{\begin{array}{lll}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 \end{array}}

係数 $c 0, ..., c N$ は次の式で与えられる^[2]。

{\begin{aligned}c_{0}&={\frac {1}{2}}p_{0}+\sum _{k=1}^{N}p_{k}\\c_{j}&=\sum _{k=j}^{N}p_{k}(-1)^{k-j+1}{\frac {(k+j-1)!}{(k-j)!((j-1)!)^{2}}}\end{aligned}}

ポール・ゴッドフリーは係数 $c 0, ..., c N$ とその級数の値を行列の積によって得る方法を述べた^[5]^[6]。

導出

簡単な実装

def gamma(z):
    
    from cmath import sin,sqrt,pi,exp
    
    # GNU Scientific Library の値を使用
    g = 7
    N = 8
    c = [ 0.99999999999980993,     676.5203681218851,   -1259.1392167224028,
           771.32342877765313,   -176.61502916214059,    12.507343278686905,
         -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7]
    
    z = complex(z)
    
    if z.real < 0.5:
        # 相補公式
        return pi / (sin(pi * z) * gamma(1 - z))
    
    z -= 1
    x  = c[0]
    
    for k in range(1, N+1):
        x += c[k]/(z + k)
    
    t = z + g + 0.5
    return sqrt(2 * pi) * t**(z + 0.5) * exp(-t) * x

脚注

^ 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.
^ ^a ^b ^c Pugh, Glendon (2004). An analysis of the Lanczos Gamma approximation (PDF) (PhD thesis).
^ ^a ^b William H. Press、Saul A. Teukolsky、William T. Vetterling、Brian P. Flannery「6.1　ガンマ関数，ベータ関数，階乗関数，２項係数」『ニューメリカルレシピ・イン・シー日本語版―C言語による数値計算のレシピ』丹慶勝市、佐藤俊郎、奥村晴彦、小林誠、169–170頁。
^ “GNU Scientific Library – Reference Manual: Gamma Functions”. 2016年10月21日閲覧。
^ “Lanczos Implementation of the Gamma Function” (2001年). 2005年6月3日閲覧。
^ “Programmable Calculators: The Lanczos Approximation” (2005年). 2006年1月11日閲覧。

外部リンク

Weisstein, Eric W. "Lanczos Approximation". mathworld.wolfram.com (英語).
GSL - GNU Scientific Library

{@{デフォルトソート:らんちよすきんし}} [@[Category:数学に関する記事]] [@[Category:数値解析]]

[1] 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.

[pughthesis-2] Pugh, Glendon (2004). An analysis of the Lanczos Gamma approximation (PDF) (PhD thesis).

[numerical_recipe_japanese-3] William H. Press、Saul A. Teukolsky、William T. Vetterling、Brian P. Flannery「6.1　ガンマ関数，ベータ関数，階乗関数，２項係数」『ニューメリカルレシピ・イン・シー日本語版―C言語による数値計算のレシピ』丹慶勝市、佐藤俊郎、奥村晴彦、小林誠、169–170頁。

[4] “GNU Scientific Library – Reference Manual: Gamma Functions”. 2016年10月21日閲覧。

[5] “Lanczos Implementation of the Gamma Function” (2001年). 2005年6月3日閲覧。

[6] “Programmable Calculators: The Lanczos Approximation” (2005年). 2006年1月11日閲覧。

[1]

[2]

[3]

[4]

[5]

[6]

導入

係数

導出

簡単な実装

関連項目

脚注

外部リンク