不完全ガンマ関数の連分数って何項で打ち切ればいいんですか?　前編

はじめに

何項打ち切りシリーズ第二弾です。第一弾は これ です。
普通のプログラミング言語はWhile文という、条件を達成したらすぐにループを抜けるという便利なものがあるのですが、Excelにはない(再帰は制限がきつい、かつ遅いため除外)ので、REDUCE文を用いてループを行います。そのため、十分性のある回数を前もって推定しなければなりません。今回は不完全ガンマ関数編です。(連分数については後方から解くこともあるので、回数の推定が出来ることに越したことはありませんね。)

本題

不完全ガンマ関数の性質は今回は語らずに、いきなり連分数を見せます。
$$\Gamma(a,x) = e^{-x}x^a\dfrac{1}{1+x-a+\dfrac{1\cdot(a-1)}{3+x-a + \dfrac{2\cdot(a-2)}{5+x-a + \dfrac{3\cdot (a-3)}{7+x-a +\cdots}}}}$$
$a_n = 2n-1+x-a\: (n\geq 1),\: b_1 = 1,\: b_n = (n-1)(a-n+1)\: (n\geq 2) $とおくと、
$$\Gamma(a,x) = e^{-x}x^a \dfrac{b_1}{a_1+}\dfrac{b_2}{a_2+}\dfrac{b_3}{a_3+}\cdots$$
と表せます。これの項数を見積もります。特に$x\sim a$(ただし$x\geq a-1, a\to\infty$を仮定)の領域について議論します。

文献によると

Numerical Recipes in C (William H. et al. 1988) にはこの様な記述があります。なお、(6.2.5) は冪級数、(6.2.6)は連分数で、(6.2.7)は収束を改善した連分数で、今回吟味する連分数に一致します。

It turns out that (6.2.5) converges rapidly for x less than about a +1, while (6.2.6) or (6.2.7) convergesrapidly for x greater than about a+1. In these respectiveregimes each requires at most a few times √a terms to converge, and this many only near x = a, where the incomplete gamma functions are varying most rapidly. Thus (6.2.5) and (6.2.7) together allow evaluation of the function for all positive
a and x.

すなわち項数は$O(\sqrt{a})$と言っていますね。
また、 Computing the Incomplete Gamma Function to Arbitrary Precision (Serge W. 2003) には

$$\sqrt{n} = \dfrac{P\ln10 +\ln(4\pi \sqrt{n})+\Re[z+\left(\frac32 -a\right)\ln z-\ln\Gamma(1-a)]}{\sqrt{8|z|+\Re(z)}}$$
To find n, it is enough to perform a few iterations of Eq. (25) starting with n =1. It is clear that at fixed z one needs $ 2n = O(P^2)$ terms of the continued fraction for $P$ digits of precision; the constant of proportionality depends on $a$ and $z$. The required number of terms grows when $z$ approaches the negative real semiaxis or zero, or when $ \Re(a) > 0$and$|a−1|> e|z|$.

これはこの連分数を多倍長演算で、しかも複素数まで拡張しているので今回の議論とは少しずれますが、項数は桁数の2乗の比例するとしか書いておらず、$a$の大きさとの詳しい議論はあまりないことが分かります。(ほかにもあったらコメント等で教えていただけると助かります。)
$$ $$
果たして本当に$O(\sqrt{a})$項なのでしょうか?
確かめていきましょう。

連分数の収束

まず連分数の収束について、 子葉 さんの 連分数１：漸化式と収束性
で以下のような式があります。

これを用います。証明は参照元の記事を確認して下さい。
つまり、これの右辺の減衰オーダーが分かれば、どれくらいで減衰するか分かるし、それから項数を抑えられる...ということですね。

ところで、この式をどう見積もりましょうか。
$\displaystyle \prod_{k=1}^{n+1}b_k$は簡単に分かります。
$$\begin{align*}\prod_{k=1}^{n+1}b_k &= \prod_{k=1}^n k(a-k)\\ &= \dfrac{n!\Gamma(a)}{\Gamma(a-n)}\end{align*}$$
問題は、$n$項目まで進んだ時の分母$q_n$です。ここで、この$q_n$について、以下の漸化式が成り立つことが分かります。

ここで、$a_n,b_n$を代入して整理すると以下の漸化式になります。

この漸化式、どこかで見たことありますね。仮に$q_n =n!L_n$と置くと、両辺$n!$で割れて、

ラゲール陪多項式(ソニン多項式) $L_n^{-a}(-x)$が出てきました。
既存の多項式だと議論しやすいので好都合です。

一旦纏めると、
$$\dfrac{(2n+5)\Gamma(a)}{(n+2)!L_n^{-a}(-x)L_{n+2}^{-a}(-x)\Gamma(a-n)}\leq \left|Z-\dfrac{p_n}{q_n}\right|\leq \dfrac{\Gamma(a)}{(n+1)!L_n^{-a}(-x)L_{n+1}^{-a}(-x)\Gamma(a-n)}$$
が成り立ちます。$L_n^{-a}(-x)$さえ何とかなれば道筋は見えそうです。

漸近展開

ラゲール陪多項式の漸近展開と言えば、Bessel型の漸近展開、特に今回の形で言うと、
$$\begin{align*}L_n^{-a}(-x) &\sim \dfrac{e^{x/2}}{n^a}I_{-a}(2\sqrt{nx})\\ &\sim \dfrac{x^{a/2-1/4}}{n^{a/2+1/4}}\exp\left(\dfrac{x}{2}+2\sqrt{nx}\right) \end{align*}$$
が有名ですが、この公式は$a$を固定して、$n\to\infty$の場合を仮定しているので不適です。
よって、オリジナルの漸近展開を考えなければなりません。
皆さん大好き鞍点法を使いましょう。これも前回のベッセル関数漸近展開の打ち切り項数の議論で出てきました。もはや準レギュラーです。

まず、Schläfliの積分表示を用意します。
$$\begin{align*}L_n^{-a}(-x) &= \dfrac{1}{2\pi i}\oint_{C}\dfrac{e^{xt/(1-t)}(1-t)^{a-1}}{t^{n+1}}\:\text{d}t\\&= \dfrac{1}{2\pi i}\oint \exp(\Phi(t))\: \text{d}t\end{align*}$$
$\left(\Phi(t) = \dfrac{xt}{1-t}+(a-1)\ln(1-t)-(n+1)\ln t\right)$
$\Phi'(t_0) = 0$なる$t_0$を求めましょう。
\begin{align*} \Phi'(t_0) &= \dfrac{x}{(1-t_0)^2}-\dfrac{a-1}{1-t_0}-\dfrac{n+1}{t_0}=0\\ \end{align*}
$t\ll 1$を仮定して、$\dfrac{x}{(1-t_0)^2}\approx x(1+2t_0),\: \dfrac{a-1}{(1-t_0)^2}\approx (a-1)(1+t_0)$と展開したものを代入すると、($x\sim a$より)
\begin{align*} \Phi'(t_0) &\approx x(1+2t_0)-(a-1)(1+t_0) - \dfrac{n+1}{t_0}\\ &\approx at_0^2-n=0 \end{align*}
$t_0 \approx \sqrt{\dfrac{n}{a}}$と求まります。これは$t_0\ll1$を十分満たすと考えてよいでしょう。
では$\Phi(t_0)$を評価しましょう。先ほどと同じように$t_0 \ll 1$として級数展開して近似等を整理すると以下のようになります。
\begin{align*}\Phi(t_0) &\approx (x-a)t_0 +\left(x-\dfrac{a}{2}\right)t_0^2 + \left(x-\dfrac{a}{3}\right)t_0^3+\left(x-\dfrac{a}{4}\right)t_0^4 - n\ln t\\ &\approx \dfrac{n}{2} + \dfrac{2n^{3/2}}{3\sqrt{a}}+\dfrac{3n^2}{4a} - \dfrac n2 \ln\left(\dfrac na\right)\end{align*}
また
\begin{align*} \Phi''(t_0) &= \dfrac{2x}{(1-t_0)^3}-\dfrac{a-1}{(1-t_0)^2}+\dfrac{(n+1)}{t_0^2}\\ &\approx 2x-a+a \approx 2a \end{align*}
最急降下方向に沿って経路を設定すると、(鞍点法から)
\begin{align*} L_{n}^{-a}(-x) &\approx \dfrac{1}{2\pi}\sqrt{\dfrac{2\pi}{|\Phi''(t_0)|}}\exp(\Phi(t_0))\\ &\approx \dfrac{1}{\sqrt{4\pi a}}\exp\left(\dfrac{n}{2}+\dfrac{2n^{3/2}}{3\sqrt{a}}+\dfrac{3n^2}{4a}-\dfrac{n}{2}\ln \left(\dfrac{n}{a}\right)\right) \end{align*}
導出出来ました。これを使って整理しましょう。対数を取っておきましょう。
$$\begin{align*}\ln\left(\dfrac{\Gamma(a)}{(n+1)!L_n^{-a}(-x)L_{n+1}^{-a}(-x)\Gamma(a-n)}\right) &\approx n\ln (a) -\dfrac{n^2}{2a} - \dfrac12 \ln(2\pi)-\left(n+\dfrac32\right)\ln (n)+n+1+\ln(4\pi)+\ln(a)-n-\dfrac{1}{2}-\dfrac{4n^{3/2}}{3\sqrt{a}}-\dfrac{3n^2}{2a}+\left(n+\dfrac12\right)\ln\left(\dfrac{n}{a}\right)\\ &\approx -\dfrac{4n^{3/2}}{3\sqrt{a}}-\dfrac{2n^2}{a}+\ln\left(\dfrac{\sqrt{a}}{n}\right)+\dfrac12 +\ln\left(2\sqrt{2\pi}\right) \end{align*}$$
これが$\ln(10^{-d})$未満になれば条件は達成されるので、
$$-\dfrac{4n^{3/2}}{3\sqrt{a}}-\dfrac{2n^2}{a}+\ln\left(\dfrac{\sqrt{a}}{n}\right)+\dfrac12 +\ln\left(2\sqrt{2\pi}\right) <-d\ln(10)$$
$-\dfrac{4n^{3/2}}{3\sqrt{a}}\sim O(1)$であることを考えると$-\dfrac{2n^2}{a}\to 0$なので無視出来て、
$$\dfrac{4n^{3/2}}{3\sqrt{a}} > d\ln(10)+\ln\left(\dfrac{2\sqrt{2\pi e a}}{n}\right)$$
これはLambertWを用いて解くことが出来ます。冗長なので省略すると
$$n_h > \left[ \frac{\sqrt{a}}{2} W_0\left( 2 (8\pi e)^{3/4} 10^{3d/2} a^{1/4} \right) \right]^{2/3}$$
オーダーを明瞭にするために$W_0(x)\sim \ln(x)$を利用します。
$$\begin{align*}n_h &> \left[ \frac{\sqrt{a}}{2} \left( \frac{1}{4}\ln(a) + \frac{3d}{2}\ln(10) + \frac{5}{2}\ln(2) + \frac{3}{4}\ln(2\pi) + \frac{3}{4} \right) \right]^{\frac{2}{3}}\sim O(a^{1/3}(\ln a)^{2/3})\end{align*}$$
つまり、誤差上界における項数$n_h$は$O(a^{1/3}\ln(a)^{2/3})$で抑えられることが分かりました。
$O(\sqrt{a})$よりもかなり小さいことが分かります。
一応下界の方も求めておきましょう。
$$\begin{align*}\ln\left(\dfrac{(2n+5)\Gamma(a)}{(n+2)!L_n^{-a}(-x)L_{n+2}^{-a}(-x)\Gamma(a-n)}\right) &\approx -\dfrac{4n^{3/2}}{3\sqrt{a}}-\dfrac{1}{2}\ln (n)+\ln(4\sqrt{2\pi})\end{align*}$$
$$n_l > \left[\dfrac{\sqrt{a}}{4}W_0\left(\dfrac{512\pi\sqrt{2\pi}\cdot 10^{3d}}{\sqrt{a}}\right)\right]^{\frac{2}{3}}$$
$$n_l > \left[\dfrac{\sqrt{a}}{4}\left(-\dfrac{1}{2}\ln(a)+3d\ln(10)+\ln(512\pi \sqrt{2\pi})\right)\right]^{\frac{2}{3}} \sim O(a^{1/3}\ln (a)^{-2/3})$$
つまり、必要な項数は、少なくとも$O(a^{1/3}\ln(a)^{-2/3})$です。(あまり意味はありませんが...。)

実際はどうなのか

ということで 数値計算bot さん(a.k.a. YK さん)の全面協力の元、$x=a$におけるこの連分数の収束項数を計算して頂きました。
収束項数(赤が$a$、緑が$d$)
図1のグラフ
項数を$\sqrt{a}$で除した値との関係
項数を$a^{1/3}$で除した値との関係
明らかに$O(a^{1/3})$であることが分かります。なんなら定数に収束しそうな感じもしますね。

どうやら
$$\dfrac{n}{a^{1/3}}\to \left(\dfrac{3}{4}\ln(4\cdot 10^{d})\right)^{2/3}$$に収束しそうです。この件は現在調査中です。
次の記事で多分明らかにします。

$x-a \to ar$ ($r$は正の定数)の場合

今までは$x\approx a$の話をしていました。もし$x-a \sim O(a)$ならばどうなるでしょうか。これは極めて高速で収束することが知られています。

では$L_n^{-a}(-x)$の鞍点法による近似からやっていきましょう。

$\Phi(t) = \dfrac{xt}{1-t}+(a-1)\ln(1-t)-(n+1)\ln (t)$です。(再掲)

$\Phi'(t_0)=0$なる$t_0$を求めます。これも、$t_0\ll 1$であると信じましょう。
$$\begin{align*}\Phi'(t_0) &=\frac{a(1+r)}{(1-t_0)^2} - \frac{a-1}{1-t_0} - \frac{n+1}{t_0}\\&\approx (a-n)t_0^2 + (ar+2n)t_0 - n \end{align*}$$
$$t_0 \approx \frac{-ar + \sqrt{a^2r^2 + 4an}}{2a} \approx \frac{n}{ar}$$
$a\to\infty$なので$t_0\ll 1$は成り立ちそうです。$n\sim o(a)$であると信じましょう。
で、$\Phi(t_0),\: \Phi''(t_0)$を評価します。
$$\begin{align*}\Phi(t_0) &= \frac{a(1+r)t_0}{1-t_0} + (a-1)\ln(1-t_0) - (n+1)\ln (t_0)\\ &\approx n + n \ln\left(\frac{ar}{n}\right) + \frac{n}{ar} + \frac{n^2(1+2r)}{2ar^2} + \ln\left(\frac{ar}{n}\right)\end{align*}$$
\begin{align*} \Phi''(t_0) &= \dfrac{2a(1+r)}{(1-t_0)^3}-\dfrac{a-1}{(1-t_0)^2}+\dfrac{(n+1)}{t_0^2}\\ &\approx \dfrac{n}{t_0^2} \approx \dfrac{a^2r^2}{n} \end{align*}
よって、
$$L_n^{-a}(-x) = \dfrac{\sqrt{n}}{\sqrt{2\pi}ar}\exp\left(n + n \ln\left(\frac{ar}{n}\right) + \frac{n}{ar} + \frac{n^2(1+2r)}{2ar^2} + \ln\left(\frac{ar}{n}\right)\right)$$

$$\ln\left(\dfrac{\Gamma(a)}{(n+1)!L_n^{-a}(-x)L_{n+1}^{-a}(-x)\Gamma(a-n)}\right) \approx n\ln\left(\dfrac{ar^2e}{n}\right)-\dfrac{5}{2}\ln(n)+2\ln(ar)+1+\dfrac{1}{2}\ln(2\pi)$$
です。いつものように
$$n\ln\left(\dfrac{ar^2e}{n}\right)-\dfrac{5}{2}\ln(n)+2\ln(ar)+1+\dfrac{1}{2}\ln(2\pi) <-d\ln(10)$$
を解きたいのですが、$-\dfrac{5}{2}\ln(n)$が邪魔です。今回は惜しいですが消えていただきましょう。
そうすると解けて、
$$n_h > -\dfrac{d\ln(10)+1+\dfrac{1}{2}\ln(2\pi)+2\ln(a)+2\ln(r)}{W_{-1}\left(-\dfrac{d\ln(10)+1+\dfrac{1}{2}\ln(2\pi)+2\ln(a)+2\ln(r)}{ar^2e}\right)}$$
$W_{-1}(-x) \sim \ln(x)$なので、
$$n_h > \frac{d\ln(10) + 1 + \frac{1}{2}\ln(2\pi) + 2\ln(a) + 2\ln(r)}{\ln(a) + 2\ln(r) + 1 - \ln\left( d\ln(10) + 1 + \frac{1}{2}\ln(2\pi) + 2\ln(a) + 2\ln(r) \right)}$$
めっちゃ見づらいのでパターンごとに見ていきましょう。
$d\ln(10)\gg 2\ln(a)$の場合、
$n_h \approx \dfrac{d\ln(10)}{\ln(a)-\ln(d)}$
つまり、$O(d)$となります。
逆に$d\ln(10)\ll 2\ln(a)$の場合
$n_h\approx \dfrac{2\ln(a)}{\ln(a)+2\ln(r)}$
つまり$n\sim O(1)$、$a$が十分大きければ$2$項で収束するという事になります。

まとめ

$x-a\sim O(1)$のときは$n\sim O(a^{1/3}\ln(a)^{2/3})$
$x-a \sim ar,\: (r>0,\text{const.})$のときは$n\sim O(d)\:\text{or}\: O(1)$

おわりに

数値計算に快く協力して頂いた 数値計算bot さん(a.k.a. YK さん)には改めて感謝申し上げます。
では、更に項数上限のオーダーを改善する後編をご期待ください。