Chudnovsky の公式を用いた円周率の計算用メモ

公式

Chudnovsky の公式は

\frac{1}{\pi} = \frac{12}{\sqrt{C^3}} \sum_{k=0}^{\infty} \frac{(-1)^k(6k)!(A+Bk)}{(3k)!(k!)^3 C^{3k}} = \frac{12}{\sqrt{C^3}} \sum_{k=0}^{\infty} P_k

\begin{cases}
A = 13591409\\
B = 545140134\\
C = 640320
\end{cases}

となっている。この等式が成り立つ証明は面倒臭いというか理解していないので省く。モジュラーだか保型形式だかよくわからない。
が、この公式の証明はされているし、これを使って計算するのが現在知られている最も効率的な方法であるのは間違いないので、Chudnovsky の公式を計算プログラムに落としこむために必要なことをメモしていこう。

計算項数の制限

まずは ∞ 項の計算はできないので、計算項数を制限するところから。この式は1項進むごとに

\frac{P_{k+1}}{P_k} = \frac{-24(6k+1)(2k+1)(6k+5)(A+Bk+B)}{C^3(k+1)^3(A+Bk)}

くらいのスケールで小さくなっていく。もう少し具体的な数値として見るために $k\rightarrow\infty$ として考えると

\lim_{k\rightarrow\infty}\frac{P_{k+1}}{P_k} = \frac{-24\cdot6^2\cdot2}{C^3}
=\frac{-1}{151931373056000}

となってる。この分母の数の常用対数値 $\log_{10}151931373056000 = 14.1816474627$ がよくサンプルプログラムの冒頭くらいに書かれていて、"１ 項進むごとに約 14 桁ずつ精度が上がる"の根拠はここにある。
ちなみにこの値(の絶対値の常用対数)は $P_1/P_0$ の 13.72607 （つまり １ 項計算するだけで 13 桁分正しいことが分かる）に始まり、 $k$ が大きくなるごとに小さくなっていき(正確になる桁数が増えていき)、上記の 14.1816474 桁/項の速さに収束する。つまり求めたい桁数 $N$ に対しての必要項数 $n$ を

n = \frac{N}{14.1816474627}

で計算していると、特に $N$ が小さいうちは足りなくなる可能性が一応ある。ので適当な数(例えば10)の追加項を含めて計算することをすすめる。

Binary Splitting Method

さて次にアルゴリズム的な話に移ろう。Chudnovsky の公式の計算では Binary Splitting Method を用いると効率よく計算できることがよく知られているというか使う前提になっている。
Binary Splitting Method というのは ３ つの整数列 $\{x_i\}$、$\{y_i\}$、$\{z_i\}$ が定義されてるときに

\frac{1}{x_0}\left(y_0+\frac{z_0}{x_1}\left(y_1+\frac{z_1}{x_2}\left(\cdots +\frac{z_{n-2}}{x_{n-1}}(y_{n-1})\right)\right)\cdots\right)

という数式の計算コストを $O(n^2)$ から $O(n(\log n)^3)$ くらいまで下げるアルゴリズムである。($m$桁同士の乗算が$O(m\log m)$でできる前提で)。 後の進行のために、単純なインデックスではなく $X(k,k+1)=x_k$ みたいな対応で範囲を使うように書き直させてもらうと、計算する式は

\frac{1}{X(0,1)} \left(Y(0,1)+\frac{Z(0,1)}{X(1,2)}\left(Y(1,2)+\frac{Z(1,2)}{X(2,3)}\left(\cdots +\frac{Z(n-2,n-1)}{X(n-1,n)}(Y(n-1,n)) \right)\right)\cdots\right)

となる。置き換えただけ。 さて、この中の一部を取り出したとき、他の部分に影響を与えないまま括弧を 1 つ無くす変形ができる。具体的には

\cdots +\frac{Z(?,l)}{X(l,m)}\left(Y(l,m)+\frac{Z(l,m)}{X(m,u)}\left(Y(m,u)+\frac{Z(m,u)}{X(u,?)}\Bigl(\cdots\right.\right.

という部分について、

\begin{align}
X(l,u) &= X(l,m)X(m,u)\\
Y(l,u) &= X(m,u)Y(l,m)+Y(m,u)Z(l,m)\\
Z(l,u) &= Z(l,m)Z(m,u)\\
\end{align}

というような計算をすると $m$ に関わる項を取り除けてこんなふうになる。

\cdots +\frac{Z(?,l)}{X(l,u)}\left(Y(l,u)+\frac{Z(l,u)}{X(u,?)}\Bigl(\right.\cdots

これを繰り返し適用していくことで、元の式を

\frac{Y(0,n-1)}{X(0,n-1)}

というシンプルな分数にするわけだが、処理自体の考え方としては逆方向に $[0,n)$ のタスクを $m=(l+u)/2$ として分割統治していく形にもとれる。ちなみに数式として直接つながりがあるわけではないが、計算方法としては

M(l,u) = \begin{pmatrix} Z(l,u) & Y(l,u) \\ 0 & X(l,u) \end{pmatrix}

と置いたときの行列積

M(0,n-1)=M(0,1)M(1,2)M(2,3)\cdots M(n-2,n-1)

と対応している。

公式の変形

さて、元々の Chudnovsky の公式に戻って、Binary Splitting Method を適用するため以下のように書き換えよう。

\begin{eqnarray}
\frac{1}{\pi} &=& \frac{12}{\sqrt{C^3}}\sum_{k=0}^{\infty} \frac{(-1)^k(6k)!(A+Bk)}{(3k)!(k!)^3 C^{3k}}\\
&=& \frac{12}{\sqrt{C^3}} \left(
\frac{A}{1} + \frac{(-1)6!(A+B)}{3!(1!)^3 C^3} +
\frac{(-1)^2 12! (A+2B)}{6!(2!)^3 C^6}\cdots\right)\\
&=& \frac{12}{\sqrt{C^3}} \left(
A +
\frac{(-1)(1\cdot2\cdot3\cdot4\cdot5\cdot6)}{(1\cdot2\cdot3)1^3 C^3}\left((A+B) + \frac{(-1)(7\cdots12)}{(4\cdots6) 2^3 C^3}
\left((A+2B) + \frac{(-1)(13\cdots18)}{(7\cdots9) \cdot 3^3 C^3}
\right)\right)\cdots\right)\\
\end{eqnarray}

対応する項を見比べると、$X$、$Y$、$Z$ を以下のように定義すればとりあえず計算ができそうである。

\begin{align}
x_k &= (3k-2)(3k-1)(3k)k^3C^3\quad(ただし\ x_0=1)\\
y_k &= A+Bk\\
z_k &= (-1)(6k+1)(6k+2)(6k+3)(6k+4)(6k+5)(6k+6)\quad(ただし\ z_{n-1}=0)\\
\end{align}

ただ、実際そのまま当てはめてみると、数式のままでも約分ができることに気づく。

\begin{eqnarray}
&& \cdots +\frac{Z(k-1,k)}{X(k,k+1)}\left(Y(k,k+1)+\frac{Z(k,k+1)}{X(k+1,k+2)}\cdots\right.\\
&=& \cdots + \frac{(-1)(6k-5)(6k-4)(6k-3)(6k-2)(6k-1)(6k)}{(3k-2)(3k-1)(3k)k^3C^3}\left(A+Bk+\frac{(-1)(6k+1)(6k+2)(6k+3)(6k+4)(6k+5)(6k+6)}{(3k+1)(3k+2)(3k+3)(k+1)^3C^3}\cdots\right.\\
&=& \cdots + \frac{-(6k-5)(2k-1)(6k-1)}{k^3C^3/24}\left(A+Bk+\frac{-(6k+1)(2k+1)(6k+5)}{(k+1)^3C^3/24}\cdots\right.\\
\end{eqnarray}

ということで、

\begin{align}
x_k &= k^3C^3/24\quad(ただし\ x_0=1)\\
y_k &= A+Bk\\
z_k &= (-1)(6k+1)(2k+1)(6k+5)\quad(ただし\ ｚ_{n-1}=0)\\
\end{align}

という形で定義すると $k$ の次数が少なくなってコンピュータ的に嬉しい。
ちなみに、$z_k$ の定義から $(-1)$ の部分を無くす代わりに

\begin{align}
X(2k,2k+2) &= x_{2k}x_{2k+1}\\
Y(2k,2k+2) &= x_{2k+1}y_{2k}-y_{2k+1}z_{2k}\\
Z(2k,2k+2) &= z_{2k}z_{2k+1}\\
\end{align}

という処理を必須で行うようにすれば $Y$ も $Z$ も常に正の値になることが保証され、負の数の処理に関わる実装が必要なくなるのでさらに嬉しい。

最終処理

さて Binary Splitting Method で分数をまとめると、最終的に円周率(有限桁)を計算するための式が

\begin{eqnarray}
\frac{1}{\pi} &=& \frac{12}{\sqrt{C^3}}\frac{Y(0,n)}{X(0,n)}\\
\therefore \pi &=& \frac{\sqrt{C^3}}{12}\frac{X(0,n)}{Y(0,n)}\\
\end{eqnarray}

もしくは

\begin{eqnarray}
\frac{1}{\pi} &=& \frac{12}{\sqrt{C^3}}\left(A-\frac{5Y(1,n)}{X(1,n)}\right)\\
\therefore \pi &=& \frac{\sqrt{C^3}}{12}\frac{X(1,n)}{AX(1,n)-5Y(1,n)}\\
\end{eqnarray}

となることが分かる。(後者のパターンだと $x_0$ を使わないので分岐をなくすことができる。) $X(0,n)$ と $Y(0,n)$ についてはまぁコンピュータに頑張ってもらうとして、残りの部分の式を展開してみる。先にネタばらししておくと平方根は逆数のままの方が計算コストが(ちょっとだけ)少ないので平方部分を取り出した後は分子の有理化処理をする。

\begin{eqnarray}
\pi &=& \frac{\sqrt{C^3}}{12}\frac{X(0,n)}{Y(0,n)}\\
&=& 53360\sqrt{640320}\frac{X(0,n)}{Y(0,n)}\\
&=& 426880 \sqrt{10005}\frac{X(0,n)}{Y(0,n)}\\
&=& \frac{4270934400 X(0,n)}{\sqrt{10005}Y(0,n)}
\end{eqnarray}

もしくは

\begin{eqnarray}
\pi &=& \frac{\sqrt{C^3}}{12}\frac{X(1,n)}{AX(1,n)-5Y(1,n)}\\
&=& \frac{53360 X(1,n)\sqrt{640320}}{AX(1,n)-5Y(1,n)}\\
&=& \frac{426880 X(1,n)\sqrt{10005}}{AX(1,n)-5Y(1,n)}\\
&=& \frac{4270934400 X(1,n)}{(AX(1,n)-5Y(1,n))\sqrt{10005}}
\end{eqnarray}

ちなみに 分子部分の 4270934400 は符号なし 32bit 整数の範囲にギリギリ収まる。この数は $X$ にとりあえずそのままかけておこう。分母部分については逆数平方根 ($1/\sqrt{10005}$) と逆数 ($1/Y$ もしくは $1/(AX-5Y)$) を求めることになる。そしてその 3 つを全部掛け合わせることで円周率の値が求まることになる。はず。

ニュートン法

ではどうやって逆数平方根や長い数の逆数を求めるかというと、ニュートン法を使う。ニュートン法は、一般的なことで言うと、方程式の根を数値的に求める方法で、適当な点から開始して、関数を一次関数で近似して解く、という処理をひたすら繰り返す処理である。うん、わけがわからない。

逆数の計算

具体的な計算として逆数を求めることを考えよう。今回の $Y(0,n)$ みたいに長い桁の数があっても、乗算と加減算だけで逆数(近似)が求められる。一応最初の例なので理論背景もくっつけつつ。
適当に与えられた数 $A$ の逆数を求めたいとき、関数

f(x)=A-\frac{1}{x}

を考える。 ちなみに $x=1/A$ で $f(x)=0$ になるのは容易に分かる。 ステップの途中として、$x_k$ がそれなりに ＄1/A＄ を近似しているとき、点 $(x_k,f(x_k))$ で $y=f(x)$ のグラフを一次展開する。

\begin{eqnarray}
y-f(x_k) &=& f'(x_k)(x-x_k)\\
y-A+\frac{1}{x_k} &=& \frac{x-x_k}{x_k^2}
\end{eqnarray}

この一次関数の根は $x=2x_k+Ax_k^2$ となるので、そのまま

x_{k+1} = 2x_k - Ax_k^2 = x_k + x_k(1 - Ax_k)

という漸化式を作る。すると $x_{k+1}$ は $x_k$ より $1/A$ に近くなっていて、正確な桁が 2 倍くらいになっているはずである。
ちなみに右の式変形をなぜしているかというと、計算途中の計算コストが下がるからである。$x_k$ が $m$ 桁くらい正確なとき、$Ax_k$ は $2m$ 桁の精度があればいい上に上半分くらいは１０進数表記であれば9が埋まっている状態で、 $1-Ax_k$ が $m$ 桁精度くらいで足りることになっている。¹ その結果、$x_k(1-Ax_k)$ の上半分が求まれば、最後の $x_k$ との足し算で本質的に足し算の演算が必要な部分は少しだけ、他は単純に後ろに追記するだけで事足りる。処理速度の向上に役立てよう。

逆数平方根の計算

逆数平方根を求める場合もまぁ似たような流れで、 $1/\sqrt{A}$ を求めたいとき、

f(x) = A - \frac{1}{x^2}

を考える。ちなみに $A$ の桁数はそんなに長くない前提にしよう。後で必要になるかわからない前提だけどとりあえず。そうすると逆数の時と同じような流れで

\begin{align}
y-f(x_k) &= f'(x_k)(x-x_k)\\
y-A+\frac{1}{x_k^2} &= \frac{2(x-x_k)}{x_k^3}\\
x_{k+1} &= x_k + \frac{1}{2}x_k(1-Ax_k^2)
\end{align}

という感じに漸化式が求まる。ちなみに $1/2$ については引数に依存するわけではないのでプログラムは専用ルーチン作ってもいい。あと、最後の部分は逆数と同じような理由で全項を $x_k/2$ でまとめることはしない。

初期値

漸化式的な処理については上記の通りで、じゃぁ最初はどうするんだというと、そこは普通の倍精度演算に任せる。結局のところ、上記の漸化処理で精度はバンバン上がっていくので初期値は素直に分かる範囲で求めるのが現実的である。ただ最初の精度は高い方がいいのは間違いないので、特に逆数の計算をする場合には上位2要素分を一旦 double などに落としこんでから初期値計算することをおすすめする。

1. TODO: 簡単な証明を書く ↩