fetburner.core

コアダンプ

ラマヌジャンの公式で円周率を計算する

要旨

ラマヌジャンの公式$$ \frac{4}{\pi}=\sum_{n=0}^{\infty} \frac{(-1)^n (4n)! (1123+21460n)}{882^{2n+1}(4^n n!)^4} $$を用いて,円周率を計算するプログラムをHaskellで書きました. メモリ8GBのPCで動作させたところ,2000万桁程度まで実用的な速度で計算できました.

動機

この間ARCの過去問を解いている時\sum_{n=0}^{L-1} a \cdot 10^{d} \quad \pmod{b} *1 をダブリングで解く問題に遭遇しまして,その解法が円周率計算の文脈でbinary splittingと呼ばれているものに似ていると思ったので,こちらも実装したくなったのが動機です.

*2 あとまぁ本日3月14日はπの日なのでちょうど良いでしょう.

Binary splittingとは?

$$ \frac{p_{1}}{q_{1}}\left(a_{1}+\frac{p_{2}}{q_{2}}\left(a_{2} + \cdots \frac{p_{n-1}}{q_{n-1}}\left(a_{n-1}+\frac{p_{n}}{q_{n}}a_{n} \right) \right) \right) $$の形をした式を高速に計算するアルゴリズムです. 今回円周率の計算に用いたラマヌジャンの公式に限らず,多くの級数はこの形に帰着できます.

ハードウェア組み込みの,いわゆる固定長で計算できる程度の精度であれば,この式を素朴に計算するだけでも十分高速に動作します. しかし円周率の計算のように,普通ハードウェアがサポートしない極端な精度を求める場合,一工夫が必要です.

何しろハードウェアがサポートしないので,筆算の要領で計算する多倍長計算のプログラムを書かなくてはなりません. N桁の数とM桁の数の和を求めたい場合,1桁目の数同士で和を取って繰り上がりを処理して,2桁目の数同士で和を取って繰り上がりを処理して…といった計算を行わなければなりません.最悪時間計算量はO(max{N, M})ですね. 固定長の数とN桁の数の積も,同様にO(N)で計算できそうです.

除算の事は今は考えず,分数の形で始めの式を計算してみましょう. a_{n}, p_{n}, p_{n}は固定長に収まると仮定しても,N桁の数とM桁の数の積はM+N桁ぐらいになるので計算の途中結果の桁数は伸びていきますね.円周率をN桁計算したい場合,Nに比例した数の項の和を取る必要がありますから,時間計算量はおよそO(n^2)になると考えられます.*3

これは桁数を増やそうとすると支障が出るぐらい重い計算量です.どうにかして計算量を減らせないでしょうか? 幸いにしてN桁の数とN桁の数の積はFFTを使えばO(N\log N)で計算できることが知られています.従って同じような桁数同士の和や積を取るように,例えばa_{1}, a_{2}, \cdots, a_{n}の積を計算したい場合,{(\cdots( ( (a_{1} \times a_{2}) \times a_{3}) \times a_{4}) \times \cdots ) \times a_{n} }と順に計算するのではなく,{(\cdots( ( (a_{1} \times a_{2}) \times (a_{3} \times a_{4})) \times \cdots ) }とトーナメント状になるように計算を進めてやれば良さそうです.

この要領で初めの式を計算するのがbinary splittingと呼ばれるアルゴリズムです. 初めの式をn項からi項だけ計算した結果の分子をx_{i},分母をy_{i}とすると,i+1項まで計算した結果は$$\frac{x_{i+1}}{y_{i+1}}=\frac{p_{n-i}}{q_{n-i}}\left(a_{n-i}+\frac{x_{i}}{y_{i}}\right)=\frac{p_{n-i}(a_{n-i}y_{i}+x_{i})}{q_{n-i}y_{i}}$$ になりますが,この式を行列を使って表現すると察しは付くでしょう. $$ \left(\begin{array}{c} x_{i+1} \\ y_{i+1} \\ \end{array}\right) = \left(\begin{array}{cc} p_{n-i} & p_{n-i}a_{n-i} \\ 0 & q_{n-i} \\ \end{array}\right) \left(\begin{array}{c} x_{i} \\ y_{i} \\ \end{array}\right) $$ 初めに挙げた式は行列の積として表現でき,それをトーナメント状に計算することで高速化するのがbinary splittingというアルゴリズムです.*4

a_{i}, p_{i}, q_{i}がうまく1桁に収まっていると仮定すると,トーナメントのi段目はi桁同士の乗算がN/2^iに比例した回数行われ,トーナメントの段数は\log N段ありますから,時間計算量はO(N (\log N)^2)まで落ちます. もっとも,ラマヌジャンの公式を用いるとa_{i}, p_{i}, q_{i}\log iに比例した桁数になってしまうため,時間計算量はO(N (\log N)^3)程度まで増えてしまいますが…

なぜHaskellを用いたのか?

最大の理由は多倍長整数を使いやすかったからです. Haskellでは多倍長整数を表すInteger型がサポートされており,Haskeller達はこれをあたかも固定長整数かのごとく気楽に使っている印象があります. それぐらい使いやすいのです.

残りの理由は僕が比較的馴染んでいたから,ぐらいでしょうか.

公式の選定

円周率の計算にはラマヌジャンの公式 $$ \frac{4}{\pi}=\sum_{n=0}^{\infty} \frac{(-1)^n (4n)! (1123+21460n)}{882^{2n+1}(4^n n!)^4} $$ を用いましたが,これは

  • 収束がそこそこ速い
  • 平方根の計算が不要

といった点が理由として挙げられます. この級数には項を1つ加えるごとに有効桁数が5桁程度ずつ増えていく性質があり,逆三角関数級数展開に比べると格段に計算が速いと言えるでしょう. もっとも,平方根の計算が必要な公式に目を向けると,有名な方のラマヌジャンの公式 $$ \frac{1}{\pi}=\frac{2\sqrt{2}}{99^2}\sum_{n=0}^{\infty} \frac{(4n)! (1103+26390n)}{(4^n 99^{n} n!)^4} $$ は1項あたり10桁程度,チュドノフスキーの公式 $$ \frac{1}{12\pi} = \sum_{n=0}^{\infty} \frac{(-1)^n(6n)!(13591409+545140134n)}{(3n)!(n!)^3 640320^{3\left(n + \frac{1}{2}\right)}} $$ に至っては14桁程度も増えるのだから恐れ入ります. 今回これらを用いなかったのは,HaskellのIntegerには平方根を任意精度で計算する関数なんて当然用意されていないので,ニュートン法なんかで再実装する手間がかかるからです.*5

実装

今回実装した円周率の計算を行うプログラムのソースコードを掲載します.

{-# LANGUAGE Strict #-}
{-# LANGUAGE StrictData #-}

binarySplitting p q r i 0 = (p i, q i * r i, r i)
binarySplitting p q r i n =
  let
     (a,  b,  c)  = binarySplitting p q r (2 * i - 1) (n - 1)
     (a', b', c') = binarySplitting p q r (2 * i) (n - 1)
  in (a * a', a * b' + b * c', c * c')

main :: IO ()
main = do
  s <- getLine
  let
    n = read s
    (_, x, y) =
      binarySplitting
        (\i -> (1 - 4 * i) * (2 * i - 1) * (4 * i - 3))
        (\i -> (21460 * i - 20337))
        (\i -> 24893568 * i * i * i) 1 $
      max 0 $ ceiling $ logBase 2 (fromIntegral n) - 1.55849716603186539
  putStrLn $ show $ 3528 * y * (10 ^ n) `div` x

binarySplitting p q r i nはその名の通りbinary splittingによって行列の積の係数を求める関数です.説明と微妙に変数名が食い違ってますが,pqrは積を求める行列の係数を表しています.そしてniはピラミッドの下から何段目の,左から何番目の行列の係数を求めるのかを表しています.

mainに目を向けると,標準入力から整数を受け取り,binarySplitting p q r i nによって円周率を計算して表示する処理を行なっています. ここではa_{i}, p_{i}, q_{i}の桁数が計算時間を左右すると言っても過言ではないので,念入りに約分を施して計算した式を与えています. 分数として出てきた計算結果をN桁の小数に直すところは手抜きですね.分子に適当に10の冪を掛けてから分母で割ってます.

性能

CPUがHaswell世代のCore i5,メモリが8GBの僕のPCに計算させたところ,1000万桁の計算は67秒程度,2000万桁の計算は153秒でした. 2000万桁を計算させた時点でメモリがカツカツで,これ以上の桁数を計算させようとすると途端に遅くなるので,ここで打ち止めにしました.

反省

一応2000万桁までは計算できましたが,tanakhさんは同じくHaskellを用いて1億桁を計算しているので物足りない結果と言えるでしょう.何がいけなかったんでしょうね〜?

どうやらtanakhさんはチュドノフスキーの公式を使っているようなので,そちらで差がついたのかもしれません. とすると,使用する公式によってどのようにして定数倍の差が付くのか分析するのも面白そうです.

それでは皆さん,良い円周率の日を.

*1:b素数とは限りません

*2:昔のC問題なので今のE問題相当?

*3:AGM系の公式は\log Nに比例した回数だけ計算を反復すれば十分ですが,あれは級数ではないので今回の話とは無関係です

*4:左下の係数は常に0なので,計算を省略できます

*5:a.k.a. 手抜き