ガンマ関数とは、みなさんおなじみの、
- \(\Gamma(z+1)=z\Gamma(z),\)
- \(\Gamma(1)=1,\)
- \(x>0\) に対し \(\log\Gamma(x)\)は凸
を満たす \(\def\Complex{\mathbf{C}}\Complex\) 上の有理型関数である。\(\DeclareMathOperator\Re{Re}\Re z>0\) に対しては、以下の積分表示がある。\[
\Gamma(z)=\int_0^\infty t^{z-1}e^{-t}dt
\]
階乗 \(n!\) との関係は\[
n!=\Gamma(n+1)
\]となる。
重要な公式としては、反転公式\[
\Gamma(1-z)\Gamma(z)=\frac{\pi}{\sin\pi z}
\]がある。特に、\(z=\frac{1}{2}\) とおけば \(\Gamma\left(\frac{1}{2}\right)=\sqrt{\pi}\) が得られる。
さて、コンピューターでガンマ関数の値を数値的に計算するにはどうすればいいのか。指数関数とか三角関数だとかは分かりやすい冪級数表示があったから良かったが、ガンマ関数にはそういうのはないのか。
いろいろググって調べた結果、ガンマ関数の近似方法として、Lanczos近似というのがあるらしい。が、ググっているだけではいまいちその実体が釈然としないし、導出方法もよくわからない。なので、Lanczosの原論文（1964）を読むことにした。大学院生という身分は便利なもので、大学の図書館でその論文にアクセスできた。
Lanczosの方法では、さっき書いた積分表示\[
\Gamma(z)=\int_0^\infty t^{z-1}e^{-t}dt
\]をいろいろ変数変換とかしていく。その際に、非負のパラメーター \(g\) を導入する（原論文では \(\gamma\) だが、この記号はオイラーの定数と紛らわしいので、ここでは \(g\) とする。Wikipediaの記事でも \(g\) としているようである）。このパラメーターによって、収束の速さとか精度とかそういうのが変わってくる。「非負の」とは言うが、実際に使うのは自然数か、自然数に \(1/2\) を加えたものぐらいだろう。
変数変換の過程で、\[
v(x)(1-\log v(x))=1-x^2,\quad v({-1})=0, v(0)=1, v(1)=e
\]という変換を行う。論文ではこの \(v(x)\) の冪級数展開の数項をいくつか計算している。自分でも係数の満たす漸化式を作ってみたが、論文を最後まで読んだ感じでは、結局使わないっぽい。
論文ではその後「テイラー級数では収束が遅い」とか言ってフーリエ展開を持ち出したり、チェビシェフ多項式（\(\cos 2m\theta\) を \(\cos\theta\) の多項式で表すやつ）を使ってトリックっぽいことをしている。
それで、肝心なところで、\(\cos\) の定積分をガンマ関数を使って表す公式\[
\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} \cos^{2z}\theta \cos 2k\theta d\theta=\sqrt{\pi}\frac{\left(z-\frac{1}{2}\right)!z!}{(z+k)!(z-k)!}
\]を持ってきていて、自分はそういうのに馴染みがなかったので面食らった。これって常識なんだろうか。っていうか、ガンマ関数をガンマ関数で表してどうする、と言いたいのだが、その後\[
\left(z-\frac{1}{2}\right)!=(\text{式変形})=\frac{\left(z-\frac{1}{2}\right)!}{z!}(\text{ほにゃらら})
\]という形に持っていって、式変形の出発点だった \(\left(z-\frac{1}{2}\right)!\) を消してしまうのだった。計算できる式が出てくるのはもともとの \(\left(z-\frac{1}{2}\right)!\) に対してではなくて新参の \(z!\) に対してという。下克上か。
で、最終的なやり方では、計算すべき式からガンマ関数（論文では階乗の表記だが）が消えておらず、おいおい、と思ったが、どうやら、整数 \(\alpha\) に対しての \(\left(\alpha-\frac{1}{2}\right)!\) は \(\left(-\frac{1}{2}\right)!=\Gamma\left(\frac{1}{2}\right)=\sqrt{\pi}\) を使って計算できるから既知でしょ、みたいな感じのようだった。
誤差の評価に関してはなんかフーリエ展開の知識が必要なようで、その辺に詳しくない自分はふ〜んそうなのか〜〜〜と思いながら読み飛ばすほかなかった。（いや、そこが一番重要なところなんだろうけど…）
で、得られる式は\[
z!=\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)
\]という形をしていて、\(A_g(z)\) は \(z\) についての有理関数（の無限級数）なのだが、その \(A_g(z)\) に出てくる係数が割と複雑な形をしている。この係数はあらかじめ計算してソースコードに埋め込んでおくことになるのだろう。
計算の方法は分かっても、
- パラメーター \(g\) を決めないといけない
- 級数をどこで打ち切るか決めないといけない
- 係数の計算が面倒くさい（精度のことを考えると、任意精度ライブラリを使って係数を計算するプログラムを別に書くべきだろう）
ので面倒くさいなあと思った。そのぐらい面倒くさがらずにちゃんとやれという話だが。
