不完全ベータ関数の連分数って何項で打ち切ればいいんですか?　後編

後編というより残業

$s<1$の方が難易度が高いし、その場合を考えていくと$s=1$と綺麗につながるので、それならば一回記事を書きなおした方が良いのではないかと思ったので後編を急遽作りました。

前編のあらすじ

$x = \dfrac{a+1}{a+b+2}$についての不完全ベータ関数の連分数
$$B_x(a,b) = \dfrac{x^a (1-x)^b}{a}\dfrac{1}{\beta_0 + \dfrac{\alpha_1}{\beta_1+\dfrac{\alpha_2}{\beta_2 + \dfrac{\alpha_3}{\beta_3+\cdots}}}}$$
$$\alpha_k = \dfrac{(a+k-1)(a+b+k-1)(b-k)k}{(a+2k-1)^2}x^2$$
$$\beta_k = a + 2k + \left(\dfrac{k(b-k)}{a+2k-1}-\dfrac{(a+k)(a+b+k)}{a+2k+1}\right)x$$
について、
推定項数$m$は
$$\left[\dfrac{3\sqrt{ab(a+b)}}{4(a+2b)}(d\ln 10 + 2\ln 2)\right]^{2/3}$$
と導出しました。今回は、$x= s\dfrac{a+1}{a+b+2},\: 0< s<1$の場合を深堀していきます。

導出

$\alpha_m \approx ams^2\dfrac{t}{t+1}$
$\beta_m \approx a(1-s)+s+2m\left(1+\dfrac{st}{t+1}\right)$
とします。$\beta_m$の近似をより深めました。
ここで、$\nu = a(1-s)+s,\: \kappa = 1+\dfrac{st}{t+1},\: \mu = \dfrac{s^2t}{t+1}$とおくと、
$$\alpha_m\approx a\mu m,\: \beta_m \approx \nu + 2\kappa m$$とすっきりします。
前編の記事と同様に、
$$\lambda_m^{\pm} = \dfrac{\beta_m}{2}\pm\sqrt{\dfrac{\beta_m^2}{4}+\alpha_m} = u_m \pm \sqrt{1+u_m^2}$$
が取れて、更に$u_k = \dfrac{\beta_k}{2\sqrt{\alpha_k}} = \dfrac{\nu+2\kappa k}{2\sqrt{a\mu k}} =\dfrac{p+qk}{\sqrt{k}}$、
$\left(ただし\: p=\dfrac{\nu}{2\sqrt{a\mu}},\: q= \dfrac{\kappa}{\sqrt{a\mu}}\right)$とすると、
$$\rho_k :=\left|\dfrac{\lambda_k^{-}}{\lambda_k^{+}}\right| = \dfrac{\sqrt{1+u_k^2}-u_k}{\sqrt{1+u_k^2}+u_k} $$
$$\ln \rho_k = -2\sinh^{-1}(u_k)$$
が成り立ちます。一方誤差上界の式は
$$\dfrac{|\lambda_m^{-}|}{C^2\sqrt{R_mR_{m+1}}}\prod_{k=1}^{m-1}\rho_k\approx \dfrac{1}{C^2}\exp \left(\int_0^m \ln \rho_k \:\text{d}k\right)$$
と近似できます。
先ほどの$u_k = \dfrac{p+qk}{\sqrt{k}}$を
$\displaystyle F(k) = \int \ln\rho_k\:\text{d}k=\int -2\sinh^{-1}(u_k)\:\text{d}k$に代入して整理すると、
$$F(k) = \dfrac{r(k)}{q}+\dfrac{4pq+1}{q^2}\tanh^{-1}\left(\dfrac{qk}{p-r(k)}\right)-2k\sinh^{-1}\left(\dfrac{p+qk}{\sqrt{k}}\right)$$
ただし、$r(k) = \sqrt{(p+qk)^2+k}$です。
$k\to 0^{+}$の収束値も確認しておくと、
$$F(0^{+}) = \dfrac{p}{q}-\dfrac{(4pq+1)\ln(4pq+1)}{2q^2}$$
となります。
$p=0$、つまり$s=1$の時の特殊形を示しておきます。
$$F_{p=0}(k) = \dfrac{\sqrt{q^2k^2+k}}{q}-\left(\dfrac{1}{q^2}+2k\right)\sinh ^{-1}\left(q\sqrt{k}\right)$$
解くべきは、
$F(m)-F(0^{+}) < -d\ln 10 - 2\ln 2\quad \left(\because C\approx\dfrac12\right)$
なので、ニュートン法を用います。$D = d\ln 10 +2\ln 2$とおくと、
$m_{i+1} = m_i+\dfrac{F(m_i)-F(0^{+})+D}{2\sinh^{-1}\left(\dfrac{p+qm_i}{\sqrt{m_i}}\right)}$
という反復式が得られます。
なお初期値は$\ln \rho_k \approx -2u_k$の近似から三次方程式に帰着させて解いた
$$P =\dfrac{p}{q},\: Q = \dfrac{3D}{8q},\: m_0 = \left(\sqrt[3]{Q + \sqrt{Q^2+P^3}}+\sqrt[3]{Q - \sqrt{Q^2+P^3}}\right)^2$$
を用いて、2回反復すれば十分です。
なお、$m >b$になったときには、$b$を$m$より少し大きい$(1.1m\: 程度)$に、漸化式を用いて引き上げてから、もう一度項数推定をして使うことを推奨します。

精度検証

$a,b,s,d$を様々な条件で試して、実際の項数と比較しました。
$a=2,5,10,20,50,100, 200, 400, 1000,3000,10000,30000,50000,100000,200000,500000$
$b=5.5,10.5,20.5,50.5,100.5,200.5,500.5,1000.5,2000.5,3000.5,10000.5,20000.5,30000.5,50000.5,100000.5,200000.5,500000.5$
$s=0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0$
$d=1, 10, 20,50, 100, 120, 200,1000,2000$
で、すべての組合せを試し、予測項数が$b$を超えたものを除外しました。
結果(差分の小さい順)
結果(差分の大きい順)
今回の結果で最大不足項数は$5$項、
安全側にはみ出した最大超過項数は$511$項となりました。
つまり、精度が$2000$桁までならば、出てきた予測項数に$+5\sim 10$をすれば十分そうです。
で、安全側に超過しやすいケースの特徴もなんとなく分かっており、$b\gg m$が守られない場合によく起こります。とはいえ、項数増加率で言えば最大$20\%$の増加でしかなく、連分数の後ろ向き計算の方が修正Lentz法よりも定数倍で強いこと、更にGPUによる計算と相性が良くなることを考えると、かなり良い結果なのかな、と思います。