音楽室と化学室と美術室とPC室の融合部屋 所謂自由室

趣味と気分で適当に色々やります.なんかあるとたまに更新します.

ガンマ関数についてのメモ&プログラミング初心者でもガンマ関数をプロットしたい

どうも数学とプログラミング初心者の☆ありゅ☆さんです.

今回はガンマ関数と積分を組んでみた後,ワイエルシュトラス形式のガンマ関数の対数形の計算結果との実関数としての対比を行った後に複素数平面全体に定義して描画してみたいと思います.理論は調べたので描きますし,自分は組んでないしブログ書いてないのでやります.この記事より性能の良さそうな実装については参考文献のシキノートさんのページを参照.

 

2020/8/10追記 漸近展開の式の途中で\left(\frac{x}{e}\right)^xが抜けていたのを修正しました.

 

 

ガンマ関数ってなーに?

基本事項

ガンマ関数については

aryuaryuaryuryu.hatenablog.com

にも少し出てきたが,階乗の一般化として名高い(?)関数で,基本として

\Gamma(s) = \int_{0}^{\infty}t^{s-1}e^{-t}dt

と表すことができる.

これについては部分積分を行うと

\int_{0}^{\infty}t^{s-1}\left(-e^{-t}\right)'dt

[-t^{s-1}e^{-t}]_{0}^{\infty} + (s-1)\int_{0}^{\infty}t^{s-2}e^{-t}dt

\therefore \Gamma(s) = (s-1)\Gamma(s-1)

となり,また

\Gamma(1) = \int_{0}^{\infty}e^{-t}dt = 1

になるため

\Gamma(4) = 3\Gamma(3) = 3\cdot 2\Gamma(2) = 3\cdot 2 \cdot 1 \Gamma(1) = 3!

といったように階乗を考えることができ,正の整数nに対し

\Gamma(n) = (n-1)!

となる.

また

\int_{0}^{\infty}\left(\frac{1}{s}t^s\right)'e^{-t}dt

として積分をすると

\Gamma(s) = \frac{\Gamma(s+1)}{s}

が得られてこれで\Gamma(s)は0以下の整数以外について定義ができる.

階乗について他にも関数を考えることができるらしいが,対数凸という性質があるものはこのガンマ関数しかないらしく,ボーア・モーレップの定理というらしい.

 

このガンマ関数の有名値としては

\Gamma\left(\frac{1}{2}\right) = \sqrt{\pi}

があり,これはガウス積分

\int_{-\infty}^{\infty}e^{-x^2}dx = \sqrt{\pi}

から示せる.

このガウス積分については

I=\int_{-\infty}^{\infty}e^{-x^2}dx

として二乗すると

I^2=\left(\int_{-\infty}^{\infty}e^{-x^2}dx\right) \cdot \left(\int_{-\infty}^{\infty}e^{-y^2}dy\right)

互いに積分は独立しているので

I^2=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}e^{-(x^2+y^2)}dxdy

と変形ができる.

ここでx = r\cos{\theta}及びy = r\sin{\theta}と置換するとヤコビアンrになるので

I^2=\int_{0}^{\infty}\int_{0}^{2\pi}e^{-r^2}rd\theta dr

I^2=2\pi\int_{0}^{\infty}e^{-r^2}r dr

I^2=-\pi[e^{-r^2}]_0^{\infty}

I^2=\pi

\therefore I=\sqrt{\pi}

からガウス積分が導出でき,結局

 \Gamma\left(\frac{1}{2}\right) = \sqrt{\pi}

\Gamma\left(\frac{1}{2}\right) = \int_{0}^{\infty}t^{-\frac{1}{2}}e^{-t}dt

について\sqrt{t}=vと置換すれば得られることがわかる.

スターリングの公式

ガンマ関数から得られる近似公式というとスターリングの公式.

n! \sim \sqrt{2\pi n}\left(\frac{n}{e}\right)^n

だがこれについてはLaplace's Methodという簡単な方法で得られる.

 

Laplace's Method

Laplace's Methodは

\int_{-\infty}^{\infty}e^{f(t)}dt

という積分について,f(x)x=aで最大値をとり双方ともx=aを中心に減少し続け負の無限に発散していく場合,最大値周辺が一番積分に影響すると考えることができ,他の部分は0に指数関数的に近づいていくことになるので

\int_{-\infty}^{\infty}e^{f(t)}dt \sim \int_{a-\epsilon}^{a+\epsilon}e^{f(t)}dt

と考えることができる.

そこでf(x)x=aを中心としたテイラー展開の先頭三項

f(x)\sim f(a) + \frac{f'(a)}{1!}(x-a) + \frac{f''(a)}{2!}(x-a)^2

と入れ替えるが,最大値を取るということはf'(a) = 0なので実際のところ2項だけになる.

こうして

\int_{-\infty}^{\infty}e^{f(t)}dt \sim \int_{-\infty}^{\infty}e^{f(a) + \frac{f''(a)}{2!}(t-a)^2}dt

と書くことができ,f(a)及びf''(a)は定数なので

e^{f(a)}\int_{-\infty}^{\infty}e^{\frac{f''(a)}{2!}(t-a)^2}dt

になり,これからガウス積分を使うと単調減少なのでこの関数はf''(a) \lt 0であり

e^{f(a)}\int_{-\infty}^{\infty}e^{-\frac{|f''(a)|}{2!}(t-a)^2)}dt

と書き換え

\sqrt{\frac{|f''(a)|}{2}}(t-a) = u

と置換すると

\int_{-\infty}^{\infty}e^{f(t)}dt \sim e^{f(a)}\sqrt{\frac{2}{|f''(a)|}}\int_{-\infty}^{\infty}e^{-u^2}du

になって,あとはガウス積分から

\int_{-\infty}^{\infty}e^{f(t)}dt \sim e^{f(a)}\sqrt{\frac{2\pi}{|f''(a)|}}

と近似結果が得られる.

 

Laplace's Methodによるスターリングの公式の導出

こうしてLaplace's Methodを得られたのでこれをガンマ関数に適用する.

\Gamma(x) = \int_{0}^{\infty}t^{x-1}e^{-t}dt

について指数関数に統一するために

\Gamma(x) = \int_{0}^{\infty}e^{(x-1)\log{t}}e^{-t}dt = \int_{0}^{\infty}e^{(x-1)\log{t}-t}dt

として,対数を消すためにt=e^vと置換すると

\Gamma(x) = \int_{-\infty}^{\infty}e^{xv-e^v}dv

になる.

これでLaplace's Methodを用いると

f(v) = xv-e^v

において

f'(v) = x-e^v

f''(v) = -e^v

になって,f'(v) = x-e^v = 0になるのはv = \log{x}になるときで,常にf''(v) \lt 0なので上に凸.

こうしてf(\log{x}) = x\log{x} - x及びf''(\log{x}) = -xなので

\Gamma(x) \sim e^{x\log{x} - x}\sqrt{\frac{2\pi}{x}}

\Gamma(x) \sim \sqrt{\frac{2\pi}{x}}\left(\frac{x}{e}\right)^x

になる.最後にx=nとして両辺にnをかけると

n! = n\Gamma(n) \sim \sqrt{2\pi n}\left(\frac{n}{e}\right)^n

が得られる.これがスターリングの公式だ.

 ガンマ関数の関連式

ガンマ関数の無限乗積形

ガンマ関数は対数凸という性質があり\log{\Gamma(x)}が下に凸の形をしている.

そこで

0\lt x \leq 1という条件で整数nに対し,始点がそれぞれに右にずれていき終点が同じ点についての傾き,\Gamma(n-1)\Gamma(n)間,\Gamma(n+x)\Gamma(n)間,\Gamma(n+1)\Gamma(n)の三種の傾きの比較をする.

そこで凸関数の性質から

\frac{\log{\Gamma(n-1)} - \log{\Gamma(n)}}{(n-1)-n} \leq \frac{\log{\Gamma(n+x)} - \log{\Gamma(n)}}{(n-x)-n} \leq \frac{\log{\Gamma(n+1)} - \log{\Gamma(n)}}{(n+1)-n}

になる.ここで式を整理すると\Gamma(n) = (n-1)!なので

\log{(n-1)} \leq \frac{\log{\Gamma(n+x)} - \log{\Gamma(n)}}{x} \leq \log{(n)}

と綺麗にまとめることができる.

こうして

\log{(n-1)^x} \leq \log{\Gamma(n+x)} - \log{\Gamma(n)} \leq \log{(n)^x}

\log{(n-1)^x} + \log{\Gamma(n)} \leq \log{\Gamma(n+x)} \leq \log{(n)^x} + \log{\Gamma(n)}

と変形し,\Gamma(n) = (n-1)!なので

\log{(n-1)^x(n-1)!} \leq \log{\Gamma(n+x)} \leq \log{(n)^x(n-1)!}

つまり

(n-1)^x(n-1)! \leq \Gamma(n+x) \leq (n)^x(n-1)!

とまとまる.

ガンマ関数の性質から

\Gamma(x) = \frac{\Gamma(1+x)}{x}

になるので

\Gamma(x) = \frac{\Gamma(1+x)}{x} = \frac{\Gamma(x+2)}{x(x+1)} = \cdots

と変形していくと

\Gamma(x) = \frac{\Gamma(x+n)}{x(x+1)(x+2)\cdots(x+n-1)}

になるので

(n-1)^x(n-1)! \leq \Gamma(n+x) \leq (n)^x(n-1)!

(n-1)^x(n-1)! \leq x(x+1)(x+2)\cdots(x+n-1)\Gamma(x) \leq (n)^x(n-1)!

と書き換えることができる.

\frac{(n-1)^x(n-1)!}{x(x+1)(x+2)\cdots(x+n-1)} \leq \Gamma(x) \leq \frac{(n)^x(n-1)!}{x(x+1)(x+2)\cdots(x+n-1)}

として(n-1)!(x+n-1)など中途半端だから

\frac{n^xn!}{x(x+1)(x+2)\cdots(x+n)}\cdot\frac{x+n}{n}\cdot \frac{(n-1)^x}{n^x}  \leq \Gamma(x) \leq \frac{n^xn!}{x(x+1)(x+2)\cdots(x+n)}\cdot \frac{x+n}{n}

と形式を揃える.

左右の不等式を分割して

\frac{n^xn!}{x(x+1)(x+2)\cdots(x+n)}\cdot\frac{x+n}{n}\cdot \frac{(n-1)^x}{n^x}  \leq \Gamma(x)

\Gamma(x) \leq \frac{n^xn!}{x(x+1)(x+2)\cdots(x+n)}\cdot \frac{x+n}{n}

揃えた分数の部分単体にする

\frac{n^xn!}{x(x+1)(x+2)\cdots(x+n)} \leq \Gamma(x)\cdot \frac{n}{x+n}\cdot \frac{n^x}{(n-1)^x}

\Gamma(x)\frac{n}{x+n} \leq \frac{n^xn!}{x(x+1)(x+2)\cdots(x+n)}

そして不等式をあわせて

\Gamma(x)\frac{n}{x+n} \leq \frac{n^xn!}{x(x+1)(x+2)\cdots(x+n)} \leq \Gamma(x)\cdot \frac{n}{x+n}\cdot \frac{n^x}{(n-1)^x}

この式についてn\to \inftyと極限をとってやると

\Gamma(x) \leq \lim_{n\to \infty}\frac{n^xn!}{x(x+1)(x+2)\cdots(x+n)} \leq \Gamma(x)

これではさみうちの原理から

\Gamma(x) = \lim_{n\to \infty}\frac{n^xn!}{x(x+1)(x+2)\cdots(x+n)}

が得られるが,0\lt x \leq 1の前提があるので定義域を拡張する.

極限を取る前の式を

\Gamma_n(x) = \frac{n^xn!}{x(x+1)(x+2)\cdots(x+n)}

として

\Gamma_n(x+1) = \frac{n^{x+1}n!}{(x+1)(x+2)\cdots(x+1+n)}

になるので

\Gamma_n(x+1) = \frac{n^{x}n!}{x(x+1)(x+2)\cdots(x+n)}\cdot \frac{nx}{x+1+n}

と変形しn\to \inftyになるので

\Gamma_n(x+1) = x\frac{n^{x}n!}{x(x+1)(x+2)\cdots(x+n)}\cdot \frac{n}{x+1+n}

\therefore \Gamma(x+1) = x\Gamma(x)

が得られ,逆に

\Gamma_n(x-1) = \frac{n^{x-1}n!}{(x-1)x(x+1)\cdots(x+n-1)}

とすると

\Gamma_n(x-1) = \frac{n^{x}n!}{x(x+1)\cdots(x+n)}\cdot \frac{(x+n)}{n(x-1)}

\therefore \Gamma_n(x-1) = \frac{\Gamma(x)}{x-1}

が得られるのでこれを繰り返していくと定義域を全ての実数に拡張できて全ての実数xについて

\Gamma(x) = \lim_{n\to \infty}\frac{n^xn!}{\prod_{k=0}^{n}(x+k)}

が説明できる.

 

ガンマ関数のワイエルシュトラス表示

ガンマ関数は0以下の整数について0除算が発生することは

\Gamma(x) = \lim_{n\to \infty}\frac{n^xn!}{\prod_{k=0}^{n}(x+k)}

の式からわかる.これについて0除算を回避するために

\frac{1}{\Gamma(x)} = \lim_{n\to \infty}\frac{\prod_{k=0}^{n}(x+k)}{n^xn!}

と逆数にする.

すると

\frac{1}{\Gamma(x)} = \lim_{n\to \infty}n^{-x}x\frac{(x+1)(x+2)\cdots(x+n)}{n!}

\frac{1}{\Gamma(x)} = \lim_{n\to \infty}x n^{-x}\left(\frac{x}{1}+1\right)\left(\frac{x}{2}+1\right)\cdots\left(\frac{x}{n}+1\right)

\frac{1}{\Gamma(x)} = \lim_{n\to \infty}x n^{-x}\prod_{k=1}^{n}\left(\frac{x}{k}+1\right)

と変形ができるがこれではn^{-x}が邪魔くさいし相乗の部分は絶対収束しない.

ここで

\prod_{k=1}^{n}\left(\frac{x}{k}+1\right)e^{q(x)}e^{-q(x)}

と絶対収束させるためのバイアスをかけてやり

\prod_{k=1}^{n}\left(\frac{x}{k}+1\right)e^{q(x)}

が収束するf(x)を考える.相乗に対数をつけて

\log\left(\prod_{k=1}^{n}\left(\frac{x}{k}+1\right)e^{q(x)}\right) = \sum_{k=1}^{n}\left(\log{\left(\frac{x}{k}+1\right)} + q(x)\right)

 ここで\log{(1-x)}テイラー展開

\log{(1-x)} = -\left(x+\frac{x^2}{2}+\frac{x^3}{3}+\frac{x^4}{4}+\cdots\right)

を使って,\left|\frac{x}{k}\right|\lt 1について

\log{(1-\frac{-x}{k})} + q(x) = -\sum_{l=1}^{\infty}\frac{1}{l}\left(\frac{-x}{k}\right)^l + q(x)

と書き換えることができるので,今後のためにk\gt 2xとなる最小のkrとおこう.

すると1\leq k \lt rの範囲の総和はせいぜい有限和なのでr\lt kについて考えよう.

これでq(x)

q(x) = \sum_{t=1}^{m}\frac{1}{t}\left(\frac{-x}{k}\right)^t

と定義してやると

\log{(1-\frac{-x}{k})} + q(x) = -\sum_{l=1}^{\infty}\frac{1}{l}\left(\frac{-x}{k}\right)^l + \sum_{t=1}^{m}\frac{1}{t}\left(\frac{-x}{k}\right)^t

\log{(1-\frac{-x}{k})} + q(x) = -\sum_{l=m+1}^{\infty}\frac{1}{l}\left(\frac{-x}{k}\right)^l

\log(1-x)テイラー展開の残りの開始位置を調整することで相乗の収束・発散をコントロールすることができそうだ.

この式は絶対値をつけて

\left|\log{(1-\frac{-x}{k})} + q(x)\right| \leq \sum_{l=m+1}^{\infty}\frac{1}{l}\left|\frac{-x}{k}\right|^l

\left|\log{(1-\frac{-x}{k})} + q(x)\right| \leq \sum_{l=m+1}^{\infty}\frac{1}{m+1}\left|\frac{-x}{k}\right|^l

\left|\log{(1-\frac{-x}{k})} + q(x)\right| \leq \frac{1}{m+1}\frac{\left|\frac{-x}{k}\right|^{m+1}}{1-\left|\frac{x}{k}\right|}

となり,現在はk\gt 2xという条件で考えているので

\left|\log{(1-\frac{-x}{k})} + q(x)\right| \leq \frac{2}{m+1}\left|\frac{-x}{k}\right|^{m+1}

と書くことができる.

こうして

\sum_{k=1}^{\infty}\left|\log{(1-\frac{-x}{k})} + q(x)\right|

について

\sum_{k=1}^{r}\left|\log{(1-\frac{-x}{k})} + q(x)\right| + \sum_{k=r+1}^{\infty}\left|\log{(1-\frac{-x}{k})} + q(x)\right|

と分割し,有限和の部分を

\sum_{k=1}^{r}\left|\log{(1-\frac{-x}{k})} + q(x)\right| = A

とすると

\sum_{k=1}^{\infty}\left|\log{(1-\frac{-x}{k})} + q(x)\right| = A+\sum_{k=r+1}^{\infty}\left|\log{(1-\frac{-x}{k})} + q(x)\right|

になる.

これで無限和は

\sum_{k=1}^{\infty}\left|\log{(1-\frac{-x}{k})} + q(x)\right| = A+\sum_{k=r+1}^{\infty}\frac{2}{m+1}\left|\frac{-x}{k}\right|^{m+1}

\sum_{k=1}^{\infty}\left|\log{(1-\frac{-x}{k})} + q(x)\right| = A+\frac{2}{m+1}\sum_{k=r+1}^{\infty}\left|\frac{-x}{k}\right|^{m+1}

\sum_{k=1}^{\infty}\left|\log{(1-\frac{-x}{k})} + q(x)\right| = A+\frac{2|-x|^{m+1}}{m+1}\sum_{k=r+1}^{\infty}\frac{1}{\left|{k}\right|^{m+1}}

これが収束する最小の整数m

\sum_{n=1}^{\infty}\frac{1}{n} \to \infty

及び

\sum_{n=1}^{\infty}\frac{1}{n^2} \lt \infty

なのでm=1のときに収束する.

 ここで収束させるためのバイアスq(x)

q(x) = \sum_{t=1}^{m}\frac{1}{t}\left(\frac{-x}{k}\right)^t

と定義しているので

q(x) = \sum_{t=1}^{1}\frac{1}{t}\left(\frac{-x}{k}\right)^t = -\frac{x}{k}

になる.

こうして

\frac{1}{\Gamma(x)} = \lim_{n\to \infty}x n^{-x}\prod_{k=1}^{n}\left(\left(1 + \frac{x}{k}\right)e^{-\frac{x}{k}}e^{\frac{x}{k}}\right)

\frac{1}{\Gamma(x)} = \lim_{n\to \infty}x \prod_{k=1}^{n}\left(\left(1 + \frac{x}{k}\right)e^{-\frac{x}{k}}e^{\frac{x}{k}}\right) e^{-x\log{n}}

\frac{1}{\Gamma(x)} = \lim_{n\to \infty}x \prod_{k=1}^{n}\left(\left(1 + \frac{x}{k}\right)e^{-\frac{x}{k}}\right) e^{x(\frac{1}{1} +\frac{1}{2} +\frac{1}{3} +\frac{1}{4} + \cdots -\log{n})}

極限を取る際オイラー定数\gamma = \frac{1}{1} +\frac{1}{2} +\frac{1}{3} +\frac{1}{4} + \cdots -\log{n}を使って

\frac{1}{\Gamma(x)} = x \prod_{k=1}^{\infty}\left(\left(1 + \frac{x}{k}\right)e^{-\frac{x}{k}}\right)e^{x\gamma}

\frac{1}{\Gamma(x)} = xe^{x\gamma} \prod_{k=1}^{\infty}\left(1 + \frac{x}{k}\right)e^{-\frac{x}{k}}

といった式が得られる.

これがガンマ関数のワイエルシュトラス表示だ.

 

ガンマ関数の相反公式

ガンマ関数のワイエルシュトラス表示

\frac{1}{\Gamma(x)} = xe^{x\gamma} \prod_{k=1}^{\infty}\left(1 + \frac{x}{k}\right)e^{-\frac{x}{k}}

-xを代入して

\frac{1}{\Gamma(-x)} = -xe^{-x\gamma} \prod_{k=1}^{\infty}\left(1 - \frac{x}{k}\right)e^{\frac{x}{k}}

こうして\frac{1}{\Gamma(x)\Gamma(-x)}を計算すると

\frac{1}{\Gamma(x)\Gamma(-x)} = -x \cdot x \prod_{k=1}^{\infty}\left(1 + \frac{x}{k}\right)e^{-\frac{x}{k}}\left(1 - \frac{x}{k}\right)e^{\frac{x}{k}}

\frac{1}{\Gamma(x)\Gamma(-x)} = -x \cdot x \prod_{k=1}^{\infty}\left(1 + \frac{x}{k}\right)\left(1 - \frac{x}{k}\right)

\frac{1}{\Gamma(x)\Gamma(-x)} = -x \cdot x \prod_{k=1}^{\infty}\left(1 - \frac{x^2}{k^2}\right)

こうして\sin{(x\pi)}の無限積表示

\sin{(x\pi)} = x\pi\prod_{k=1}^{\infty}\left(1 - \frac{x^2}{k^2}\right)

を使うと

\frac{1}{\Gamma(x)\Gamma(-x)} = -x  \frac{\sin{(x\pi)}}{\pi}

\frac{1}{\Gamma(x)\left(-x\Gamma(-x)\right)} =  \frac{\sin{(x\pi)}}{\pi}

\frac{1}{\Gamma(x)\Gamma(-x+1)} =  \frac{\sin{(x\pi)}}{\pi}

\Gamma(x)\Gamma(1-x) =  \frac{\pi}{\sin{(x\pi)}}

が得られる.

ちなみにx=\frac{1}{2}を代入すると

\Gamma\left(\frac{1}{2}\right)^2 =  \pi

\therefore \Gamma\left(\frac{1}{2}\right) =  \sqrt{\pi}

が得られる.

 

漸近展開

ベルヌーイ多項式

ベルヌーイ多項式って何?

ベルヌーイ数によくにた定義で

\frac{te^{xt}}{e^t-1}

t=0を中心としたテイラー展開

\frac{te^{xt}}{e^t-1} = \sum_{n=0}^{\infty}\frac{B_n(x)}{n!}t^n

で定義される関数B_n(x)をベルヌーイ多項式と言い,言われてみるとわかると思うがx=0のときはベルヌーイ数B_nになる.

 

ベルヌーイ多項式の関係式

これは

\frac{t}{e^t-1}e^{xt}

についてテイラー展開に置き換えると

\frac{t}{e^t-1}e^{xt} = \sum_{k=0}^{\infty}\frac{B_k}{k!}t^k \sum_{l=0}^{\infty}\frac{1}{l!}(xt)^l

としてこれを

\frac{t}{e^t-1}e^{xt} = \sum_{k=0}^{\infty}\sum_{l=0}^{\infty}\frac{B_k}{k!}t^{k+l} \frac{1}{l!}x^l

と変形し,k+l = nと置換して

\frac{t}{e^t-1}e^{xt} = \sum_{n=0}^{\infty}\sum_{k=0}^{n}\frac{B_k}{k!} \frac{1}{(n-k)!}x^{n-k} t^{n}

\frac{t}{e^t-1}e^{xt} = \sum_{n=0}^{\infty}\sum_{k=0}^{n}\left(\frac{B_k}{k!(n-k)!} x^{n-k}\right) t^{n}

\frac{t}{e^t-1}e^{xt} = \sum_{n=0}^{\infty}\sum_{k=0}^{n}\left(\frac{n!B_k}{k!(n-k)!} x^{n-k}\right) \frac{t^{n}}{n!}

\therefore \frac{t}{e^t-1}e^{xt} = \sum_{n=0}^{\infty}B_n(x) \frac{t^{n}}{n!}

とつまり

B_n(x) = \sum_{k=0}^{n}\left(\tbinom{n}{k}B_k x^{n-k}\right)

の関係式が綺麗に得られる.

ちなみにn=0,1,2のときは

B_0(x) = 1

B_1(x) = x-\frac{1}{2}

B_2(x) = x^2-x+\frac{1}{6}

といった式になる.

ベルヌーイ多項式微分

得られた式

B_n(x) = \sum_{k=0}^{n}\left(\tbinom{n}{k}B_k x^{n-k}\right)

について微分を行うと,k=nの時定数なので

(B_n(x))' = \sum_{k=0}^{n-1}\left(\frac{n!(n-k)}{k!(n-k)!}B_k x^{n-k-1}\right)

(B_n(x))' = \sum_{k=0}^{n-1}\left(\frac{n!}{k!(n-k-1)!}B_k x^{n-k-1}\right)

(B_n(x))' = n\sum_{k=0}^{n-1}\left(\frac{(n-1)!}{k!( (n-1)-k)!}B_k x^{n-k-1}\right)

(B_n(x))' = nB_{n-1}(x)

微分結果が得られる.

 

オイラー・マクローリンの定理

スターリングの公式はLaplace's Methodを使って出したが,ちゃんとした漸近展開を出すならオイラー・マクローリンの和公式を使うのが比較的楽.

オイラーマクローリンの和公式の最初の段階は,積分と幅1の長方形で積分を近似した場合(右側の区分求積法みたいな感じ)の誤差についての式で,これはまず整数rに対して

\int_{r-1}^{r}(t-[t])f'(t)dt

についての積分を計算する.r-1\leq t \lt rの間は[t]=r-1なので

\int_{r-1}^{r}(t-r+1)f'(t)dt

と書き換えて部分積分を行う.

\int_{r-1}^{r}(t-[t])f'(t)dt = [(t-r+1)f(t)]_{r-1}^{r} - \int_{r-1}^{r}f(t)dt

\int_{r-1}^{r}(t-[t])f'(t)dt = f(r) - \int_{r-1}^{r}f(t)dt

これをm+1からnまで足していくと

\int_{m}^{n}(t-[t])f'(t)dt = \sum_{k=m+1}^{n}f(k) - \int_{m}^{n}f(t)dt

が得られる.

次にこの左辺の積分を少し変えて

\int_{m}^{n}\left(t-[t]-\frac{1}{2}\right)f'(t)dt

として

\int_{m}^{n}\left(t-[t]-\frac{1}{2}\right)f'(t)dt = \int_{m}^{n}(t-[t])f'(t)dt - \frac{1}{2}\int_{m}^{n}f'(t)dt

と展開してさっきの結果を入れる

\int_{m}^{n}\left(t-[t]-\frac{1}{2}\right)f'(t)dt = \sum_{k=m+1}^{n}f(k) - \int_{m}^{n}f(t)dt - \frac{1}{2}(f(n)-f(m))

式を整理すると

\sum_{k=m+1}^{n}f(k) - \int_{m}^{n}f(t)dt = \int_{m}^{n}\left(t-[t]-\frac{1}{2}\right)f'(t)dt + \frac{1}{2}(f(n)-f(m))

この結果はベルヌーイ数とベルヌーイ多項式を用いて

\sum_{k=m+1}^{n}f(k) - \int_{m}^{n}f(t)dt = - B_{1}\cdot (f(n)-f(m)) + \int_{m}^{n}B_1(t-[t])f'(t)dt

と表すことができ,ここで右辺の積分について(B_n(x))' = nB_{n-1}(x)及びB_n(0) = B_nなので

\int_{m}^{n}B_1(t-[t])f'(t)dt = \int_{m}^{n}\frac{1}{2}(B_2(t-[t]))'f'(t)dt

\int_{m}^{n}B_1(t-[t])f'(t)dt = \frac{1}{2}[B_2(t-[t])f'(t)dt]_{m}^{n} - \frac{1}{2}\int_{m}^{n}B_2(t-[t])f''(t)dt

\int_{m}^{n}B_1(t-[t])f'(t)dt = \frac{1}{2}B_2(f'(n)-f'(m)) - \frac{1}{2}\int_{m}^{n}B_2(t-[t])f''(t)dt

となる.これを繰り返すと正の整数aに対し

\int_{m}^{n}B_{a-1}(t-[t])f'(t)dt = \frac{1}{a!}B_a(f'(n)-f'(m)) - \frac{1}{a!}\int_{m}^{n}B_a(t-[t])f^{(a)}(t)dt

が得られるので

\sum_{k=m+1}^{n}f(k) - \int_{m}^{n}f(t)dt = \sum_{k=1}^{N}(-1)^k \frac{B_{k}}{k!}\cdot (f^{(k-1)}(n)-f^{(k-1)}(m)) + \frac{(-1)^{N-1}}{k!}\int_{m}^{n}B_N(t-[t])f^{(N)}(t)dt

になる.

 

漸近展開の計算

このオイラー・マクローリン定理を

 \Gamma(x) = \lim_{n\to \infty}\frac{n^xn!}{\prod_{k=0}^{n}(x+k)}

に用いるのが良いらしい.

 \Gamma(x) = \lim_{n\to \infty}\left(\frac{n^x}{x+n}\prod_{k=1}^{n}\frac{k}{x+k-1}\right)

 と変形して

\log{\Gamma(x)} = \lim_{n\to \infty}\left( (x-1)\log{n} + \log{\frac{1}{\frac{x}{n}+1}} - \sum_{k=1}^{n}\frac{x+k-1}{k}\right)

 と対数で分解する.\log{\frac{1}{\frac{x}{n}+1}}n\to \inftyのとき0になるから除外して考えて

\log{\Gamma(x)} = \lim_{n\to \infty}\left( (x-1)\log{n} - \sum_{k=1}^{n}\log{\frac{x+k-1}{k}}\right)

になる.

これで

\sum_{k=1}^{n}\log{\frac{x+k-1}{k}}

についてオイラーマクローリンの定理

\sum_{k=m+1}^{n}f(k) - \int_{m}^{n}f(t)dt = \sum_{k=1}^{N}(-1)^k \frac{B_{k}}{k!}\cdot (f^{(k-1)}(n)-f^{(k-1)}(m)) + \frac{(-1)^{N-1}}{k!}\int_{m}^{n}B_N(t-[t])f^{(N)}(t)dt

を用いる.

\sum_{k=1}^{n}\log{\frac{x+k-1}{k}} = \log{x} + \sum_{l=2}^{n}\log{\frac{l+s-1}{l}}

と少しだけ総和を分けて

\sum_{k=1}^{n}\log{\frac{x+k-1}{k}} = \log{x} + \sum_{l=2}^{n}\log{\frac{l+s-1}{l}}

ここで

(\log{x})^{k} = (k-1)!x^{-k}(-1)^{k-1}

になるので,m=1とすると

\sum_{k=1}^{n}\log{\frac{x+k-1}{k}} = \log{x} + \int_{1}^{n}\log{\frac{t+x-1}{t}}dt - B_1(\log{(n+x-1)} - \log{n} - \log{x} + \log{1}) \\\\+ \sum_{k=2}^{N}(-1)^k\frac{B_k}{k!}\left( (-1)^{k-2}(k-2)!( (n+x-1)^{-(k-1)} - n^{-(k-1)}) - (-1)^{k-2}(k-2)!( (1+x-1)^{-(k-1)} - 1^{-(k-1)})\right) \\\\+ \frac{(-1)^{N-1}}{N!}\int_{1}^{n}B_N(t-[t])f^{(N)}(t)dt 

になる.めっちゃ長い・・・これを整理して

\sum_{k=1}^{n}\log{\frac{x+k-1}{k}} = \log{x} + \int_{1}^{n}\log{\frac{t+x-1}{t}}dt - B_1(\log{(n+x-1)} - \log{n} - \log{x}) \\\\+ \sum_{k=2}^{N}\frac{B_k}{k(k-1)}\left( (n+x-1)^{-(k-1)} - n^{-(k-1)} - x^{-(k-1)} + 1\right) \\\\+ \frac{(-1)^{N-1}}{N!}\int_{1}^{n}B_N(t-[t])f^{(N)}(t)dt

ここで右辺について

\int_{1}^{n}\log{\frac{t+x-1}{t}}dt = \int_1^{n}\left(\log{(t+x-1)} - \log{t}\right)dt

\int_{1}^{n}\log{\frac{t+x-1}{t}}dt = \int_1^{n}(t+x-1)'\log{(t+x-1)}dt - \int_{1}^{n}t'\log{t}dt

\int_{1}^{n}\log{\frac{t+x-1}{t}}dt = \int_1^{n}(t+x-1)'\log{(t+x-1)}dt - \int_{1}^{n}t'\log{t}dt

\int_{1}^{n}\log{\frac{t+x-1}{t}}dt = [(t+x-1)\log{(t+x-1)}]_1^n - [t\log{t}]_1^n

と変形することで

\int_{1}^{n}\log{\frac{t+x-1}{t}}dt = (n+s-1)\log{(n+x-1)} - x\log{x} - n\log{n}

になり,またベルヌーイ数B_1 = -\frac{1}{2}を用いると

\sum_{k=1}^{n}\log{\frac{x+k-1}{k}} = \log{x} + (n+x-1)\log{(n+x-1)} - x\log{x} - n\log{n} + \frac{\log{(n+x-1)} - \log{n} - \log{x}}{2} \\\\+ \sum_{k=2}^{N}\frac{B_k}{k(k-1)}\left( (n+x-1)^{-(k-1)} - n^{-(k-1)} - x^{-(k-1)} + 1\right) \\\\+ \frac{(-1)^{N-1}}{N!}\int_{1}^{n}B_N(t-[t])f^{(N)}(t)dt

整理したもの

\sum_{k=1}^{n}\log{\frac{x+k-1}{k}} = \left(\frac{1}{2}-x\right)\log{x} + \left(n+x-\frac{1}{2}\right)\log{(n+x-1)} - \left(n+\frac{1}{2}\right)\log{n} \\\\+ \sum_{k=2}^{N}\frac{B_k}{k(k-1)}\left( (n+x-1)^{-(k-1)} - n^{-(k-1)} - x^{-(k-1)} + 1\right) + \frac{(-1)^{N-1}}{N!}\int_{1}^{n}B_N(t-[t])f^{(N)}(t)dt

\log{\Gamma(x)} = \lim_{n\to \infty}\left( (x-1)\log{n} - \sum_{k=1}^{n}\log{\frac{x+k-1}{k}}\right)

に入れる.

\log{\Gamma(x)} = \lim_{n\to \infty}\left( (x-1)\log{n} - \left(\frac{1}{2}-x\right)\log{x} - \left(n+x-\frac{1}{2}\right)\log{(n+x-1)} + \left(n+\frac{1}{2}\right)\log{n} \\\\- \sum_{k=2}^{N}\frac{B_k}{k(k-1)}\left( (n+x-1)^{-(k-1)} - n^{-(k-1)} - x^{-(k-1)} + 1\right) - \frac{(-1)^{N-1}}{N!}\int_{1}^{n}B_N(t-[t])f^{(N)}(t)dt\right)

式を整理して

\log{\Gamma(x)} = \lim_{n\to \infty}\left(\left(x-\frac{1}{2}\right)\log{x} - \left(n+x-\frac{1}{2}\right)\log{(n+x-1)} + \left(n + x-\frac{1}{2}\right)\log{n} \\\\- \sum_{k=2}^{N}\frac{B_k}{k(k-1)}\left( (n+x-1)^{-(k-1)} - n^{-(k-1)} - x^{-(k-1)} + 1\right) - \frac{(-1)^{N-1}}{N!}\int_{1}^{n}B_N(t-[t])f^{(N)}(t)dt\right)

\log{\Gamma(x)} = \lim_{n\to \infty}\left(\left(x-\frac{1}{2}\right)\log{x} + \left(n+x-\frac{1}{2}\right)\log{\frac{n}{n+x-1}} \\\\- \sum_{k=2}^{N}\frac{B_k}{k(k-1)}\left( (n+x-1)^{-(k-1)} - n^{-(k-1)} - x^{-(k-1)} + 1\right) - \frac{(-1)^{N-1}}{N!}\int_{1}^{n}B_N(t-[t])f^{(N)}(t)dt\right)

最後にn\to \inftyを取る際

\left(n+x-\frac{1}{2}\right)\log{\frac{n}{n+x-1}} = -\log{\left(\frac{n+x-1}{n} \right) ^ {n+x-\frac{1}{2}}}

\left(n+x-\frac{1}{2}\right)\log{\frac{n}{n+x-1}} = -\log{\left(1 + \frac{x-1}{n} \right) ^ {n+x-\frac{1}{2}}}

\left(n+x-\frac{1}{2}\right)\log{\frac{n}{n+x-1}} = -\log{\left(1 + \frac{x-1}{n} \right) ^ {n} \left(1 + \frac{x-1}{n} \right)^{x-\frac{1}{2}}}

\left(n+x-\frac{1}{2}\right)\log{\frac{n}{n+x-1}} = -\log{\left(\left(1 + \frac{x-1}{n} \right) ^ {\frac{n}{x-1}}\right)^{x-1} \left(1 + \frac{x-1}{n} \right)^{x-\frac{1}{2}}}

n\to \inftyのとき

\left(n+x-\frac{1}{2}\right)\log{\frac{n}{n+x-1}} = -\log{e^{x-1}}

\left(n+x-\frac{1}{2}\right)\log{\frac{n}{n+x-1}} = 1-x

に注意すると

\log{\Gamma(x)} = \left(x-\frac{1}{2}\right)\log{x} + 1 - x - \sum_{k=2}^{N}\frac{B_k}{k(k-1)}\left(1 - x^{1-k}\right) - \frac{(-1)^{N-1}}{N!}\int_{1}^{\infty}B_N(t-[t])f^{(N)}(t)dt

が得られる.

これは

\log{\Gamma(x)} = \left(x-\frac{1}{2}\right)\log{x} - x + \sum_{k=2}^{N}\frac{B_k}{k(k-1)}x^{1-k} + \left(1-\sum_{k=2}^{N}\frac{B_k}{k(k-1)} - \frac{(-1)^{N-1}}{N!}\int_{1}^{\infty}B_N(t-[t])f^{(N)}(t)dt\right)

\Gamma(x) = \frac{1}{\sqrt{x}}\left(\frac{x}{e}\right)^x \cdot \exp\left(1-\sum_{k=2}^{N}\frac{B_k}{k(k-1)} - \frac{(-1)^{N-1}}{N!}\int_{1}^{\infty}B_N(t-[t])f^{(N)}(t)dt\right) \cdot  \exp\left(\sum_{k=2}^{N}\frac{B_k}{k(k-1)}x^{1-k}\right)

として定数項

\left(1-\sum_{k=2}^{N}\frac{B_k}{k(k-1)} - \frac{(-1)^{N-1}}{N!}\int_{1}^{\infty}B_N(t-[t])f^{(N)}(t)dt\right)

についてはスターリングの公式

\Gamma(x) \sim \sqrt{\frac{2\pi}{x}}\left(\frac{x}{e}\right)^x

よりN\to \inftyのとき

\left(1-\sum_{k=2}^{N}\frac{B_k}{k(k-1)} - \frac{(-1)^{N-1}}{N!}\int_{1}^{\infty}B_N(t-[t])f^{(N)}(t)dt\right) = \frac{1}{2}\log{(2\pi)}

がわかるので

\log{\Gamma(x)} = \left(x-\frac{1}{2}\right)\log{x} + \frac{1}{2}\log{(2\pi)} - x + \sum_{k=2}^{\infty}\frac{B_k}{k(k-1)}x^{1-k}

が得られ,

\Gamma(x) = \exp\left(\left(x-\frac{1}{2}\right)\log{x} + \frac{1}{2}\log{(2\pi)} - x + \sum_{k=2}^{\infty}\frac{B_k}{k(k-1)}x^{1-k}\right)

\Gamma(x) = \sqrt{\frac{2\pi}{x}}\left(\frac{x}{e}\right)^x\exp\left(\sum_{k=2}^{\infty}\frac{B_k}{k(k-1)}x^{1-k}\right)

とまとめることができる.これを展開すると以前求めたベルヌーイ数

aryuaryuaryuryu.hatenablog.com

を先頭10個ほど

B_2 = \frac{1}{6}
B_3 = 0
B_4 = -\frac{1}{30}
B_5 = 0
B_6 = \frac{1}{42}
B_7 = 0
B_8 = -\frac{1}{30}
B_9 = 0
B_{10} = \frac{5}{66}

B_{11} = 0
B_{12} = -\frac{691}{2730}
を使うと

\Gamma(x) = \sqrt{\frac{2\pi}{x}}\left(\frac{x}{e}\right)^x\exp\left(\frac{1}{12x} - \frac{1}{360x^3} + \frac{1}{1260x^5} - \frac{1}{1680x^7} + \frac{1}{1188x^9} - \frac{691}{360360x^{11}} + \cdots \right)

になる.

プログラムを組んでみる

積分形式

ガンマ関数について色々見てきたからそろそろプログラムを組んでみたい.

最初は簡単なガンマ関数の形の積分表示

\Gamma(s) = \int_{0}^{\infty}t^{s-1}e^{-t}dt

を組んでみる.

 

積分の実装というと,直感的な方法は長方形近似

f:id:Aryuaryuaryuryu:20200211181854p:plain

長方形近似

\sum f(t)\cdot \Delta t

で,これだと図を見ての通り幅のとり方によっては誤差が大きくなるが,台形近似

f:id:Aryuaryuaryuryu:20200211182145p:plain

台形近似

\sum \frac{(f(t) + f(t+\Delta t))\cdot \Delta t}{2}

を用いると少しはマシになる.

ちなみにこの刻み幅\Delta tは幅がh倍になると誤差はh^2倍になる.

 

今後なんか積分計算したいときに使い回せるようにもしたいからそれも考えて実装すると

/*積分開始,積分終了,被積分関数,被積分関数の引数,積分間隔*/

double Integral(double Min, double Max, std::function<double(double,std::vector<double>)> fn, std::vector<double> v, double interval = 0.01) {

double result = 0, i;

for (i = Min; i < Max; i += interval)

result += ( (fn(i, v) + fn(i + interval, v)) * interval) / 2.;

result += ( (fn(i - interval, v) + fn(Max, v)) * (Max - i + interval)) / 2.;

return result;

}

として

\int_{Min}^{Max}f(t,v)dt

についてf(t,v)を実装すると良いようにした.

これで

#define M_E (2.71828182845904)

double Gamma1(double t,std::vector<double> v) {

double X = v[0];

return pow(t, X - 1.) * pow(M_E, -t);

}

として計算を行うと恐ろしいほど計算時間が遅いガンマ関数が実装できた.

この実装の性能としては

積分開始位置0.0001,

積分終了位置10000

■台形の幅0.001

\Gamma(0.5)を計算するのにおおよそ1分弱,計算精度は\sqrt{\pi} = 1.7724538\cdotに対し計算結果は1.7723とまぁまぁボチボチな結果が得られた.

わかってはいたけど控えめに言ってクソ.

そこでs=0.5のときの被積分関数のグラフを見てみると

f:id:Aryuaryuaryuryu:20200211222123p:plain

x=0.5のときの被積分関数

になる.t=0\to 0.5の範囲はすごく急に変化しているのに対しその後はゆったりしてt=4の時点でめっちゃゼロ.

 じゃぁ積分範囲を小さくして

積分開始位置0.0001,

積分終了位置5

■台形の幅0.00001

とすると結果は2.5秒ほどかかって1.76348・・・うーん・・・

 

ただs=5のときは

積分開始位置0.0001,

積分終了位置100

■台形の幅0.001

の設定では2msで4! = 24が得られた.

これは多分

f:id:Aryuaryuaryuryu:20200211224357p:plain

s=5のときの被積分関数

のようなゆったりとしたカーブを描いているのと,台形近似による誤差で少し多めに計算されたのが積分を途中で打ち切っているのを埋めているのかなと思う.

 

そこでs=7の場合は

f:id:Aryuaryuaryuryu:20200211224629p:plain

s=7のときの被積分関数

という結構急なカーブを描いているのに720と意外とちゃんとした値が得られている.

じゃぁ多分t=0~1の範囲が上に跳ね上がるのがいけないのかな?というと

f:id:Aryuaryuaryuryu:20200211225112p:plain

x=1のときの被積分関数

今度は結果は0.99999と結構1に近い結果が得られた.やっぱり0.5のときのスタート地点がいけなかったのかな?

といろいろスタート地点等をいじったけど最終的にはやっぱ刻み幅によりそうなので断念.ただある程度大きい数値に対してはそれなりにアバウトな設定でも速くていい感じの値が得られた.

 

無限乗積形式

ガンマ関数の無限乗積表示

\Gamma(x) = \lim_{n\to \infty}\frac{n^xn!}{\prod_{k=0}^{n}(x+k)}

について,これを実装してみると

double Gamma3(double x,double Max = 50000) {

double result = pow(Max, x) / x;

for (int i = 1; i < 50000; i++)result *= (double)i / (i + x);

return result;

}

となる.

正直ネタのつもりだったからあらかたこの記事を書いた後で実装したら,対数形式とあまり性能は変わらなかったと思う.

正直対数と実装はあまり変わってないからそりゃそうだよね・・・

 

対数形式

ガンマ関数のワイエルシュトラス表示形式

\frac{1}{\Gamma(x)} = xe^{x\gamma} \prod_{k=1}^{\infty}\left(1 + \frac{x}{k}\right)e^{-\frac{x}{k}}

に対数をつけて

\log{\Gamma(x)} = -\log{x} - x\gamma -  \sum_{k=1}^{\infty}\log{\left(1 + \frac{x}{k}\right)} -\frac{x}{k}

となるので,\log{\Gamma(x)}を計算して後で\expにぶちこめば良いのでは?

ということで

double Gamma2(double x) {

double result = -log(x) - x * M_GAMMA;

for (int i = 1; i < 50000; i++)result -= log(1. + x / i) - x / i;

 return pow(M_E,result);

}

と実装した.

同様に\Gamma(0.5)を計算するのにおおよそ0.002秒,計算精度は\sqrt{\pi} = 1.7724538\cdotに対し計算結果は1.77245が得られた.さっきの積分形式に比べると恐ろしく性能がよい.

 

そこでx=5を代入した場合,理論値は24だが23.994が出てきて誤差0.006に.

x=9を代入した結果理論値は362880だが結果は40287.4にと代入値が大きくなると誤差が大きくなっていった.

これについて考えてみると

\log{\Gamma(x)} = -\log{x} - x\gamma -  \sum_{k=1}^{\infty}\left(\log{\left(1 + \frac{x}{k}\right)} -\frac{x}{k}\right)

についてx\lt\lt 1のとき\log{(1+x)}テイラー展開

\log{(1+x)} = x-\frac{x^2}{2}+\frac{x^3}{3}-\frac{x^4}{4}+\frac{x^5}{5}+\cdots

を利用して

\log{\Gamma(x)} = -\log{x} - x\gamma -  \sum_{k=1}^{\infty}\left(-\frac{\left(\frac{x}{k}\right)^2}{2} + \frac{\left(\frac{x}{k}\right)^3}{3} - \frac{\left(\frac{x}{k}\right)^4}{4} + \cdots\right)

と変形して

\Gamma(x) = x^{-1}e^{-x\gamma}\cdot \prod_{k=1}^{\infty}\exp\left(\frac{\left(\frac{x}{k}\right)^2}{2} - \frac{\left(\frac{x}{k}\right)^3}{3} + \frac{\left(\frac{x}{k}\right)^4}{4} - \cdots\right)

と書き換えることができる.

xが1未満だとこの\expの内部の総和が早く収束するので,引数が小さい場合は結構いい感じの値が得やすいのかな?と思う.

そこで1未満の値の場合はそのまま出力して,1超過の場合は

\Gamma(s+1) = s\Gamma(s)

函数等式を使うのが良さ気で,0未満の場合は

\Gamma(x)\Gamma(1-x) = \frac{\pi}{\sin{\pi x}}

を使えそう.

ただこれは実軸方向への拡張なので,虚数方向に移動するとなると別の方法が必要になる.

 

漸近展開による方法

ベルヌーイ数は予め求めてたからいきなり値ぶっこんでおくw

そこで漸近展開を用いて計算するなら簡単に書くと

double Gamma4(double x) {

std::vector<double> Bernoulli;

Bernoulli.push_back(1./12.); Bernoulli.push_back(0);

Bernoulli.push_back(-1./360.); Bernoulli.push_back(0);

Bernoulli.push_back(1./1260.); Bernoulli.push_back(0);

Bernoulli.push_back(-1./1680.); Bernoulli.push_back(0);

Bernoulli.push_back(1./1188.); Bernoulli.push_back(0);

Bernoulli.push_back(-691./360360.); Bernoulli.push_back(0);

Bernoulli.push_back(0); Bernoulli.push_back(1./156.);

double Result = sqrt(2. * M_PI / x) * pow(x / M_E, x);

double tmp = 0;

for (int i = 1; i <= Bernoulli.size(); i++) {

double temp=1;

for (int j = 0; j < i; j++)temp *= x;

tmp += Bernoulli[i - 1]/temp;

}

return Result * exp(tmp);

}

になる. 

引数が1以上なら結構いい具合の値が出てくる.1以上に対しての計算はもうこれでいいのではないか

 絶対値が1未満なら対数形式で,1超過ならこれで十分かも.

 

複素数平面全体について定義する

基本的に今まで作った型の引数等を『std::complex<double>』に変えてやるだけ.

0以下の整数はとりあえず計算最大値INFTY(適当に決める)を返して,正の整数は階乗返しとけば良いかも.

引数の絶対値が1未満では対数形式,1超過だと漸近展開形式が良さ気.絶対値が1の場合は(1+i)を代入して比較した感じ漸近展開形式がマシだろうはと思う.

 ゼータ関数と同じようにDxlibで描画.|\Gamma(s)|は0~4の範囲で,\Re(\Gamma(s))及び\Im(\Gamma(s))は-4~4の範囲で描画してみた.

あとは偏角も描画してみた.こっちはHSVのHを偏角として,S=V=1で描画した.

f:id:Aryuaryuaryuryu:20200214003946p:plain

|Γ(s)|の描画

f:id:Aryuaryuaryuryu:20200214004103p:plain

Re(Γ(s))の描画

f:id:Aryuaryuaryuryu:20200214004156p:plain

Im(Γ(s))の描画

f:id:Aryuaryuaryuryu:20200214040133p:plain

Γ(s)の偏角

 

参考文献

スターリングの公式の簡単な導出法

ガンマ関数の数値計算 | シキノート

素数とゼータ関数

複素解析