円周率の望みの桁を直接計算する法

円周率の計算は、3. 14159265... と、数字の左から順番に始めて、
だんだん小数点以下の桁数の大きい方に進めてゆくのが常套手段です。
ところが、この常識を覆して、円周率の望みの桁数の数字をいきなり直接計算できる方法があるのだそうです。
それは「BBPアルゴリズム」と呼ばれています。
例えば円周率の１００万桁目の数字が知りたかったら、最初の３．１４はすっ飛ばして、
いきなり１００万桁目の数字がわかってしまう、ということなのです。
ただし、この計算は１６進数でないとできません。
なので１００万桁目と言っても、普通に１０進数での１００万桁ではなく、１６進数で書いたときの１００万桁目なのです。
・・・そんな計算、本当にできるのだろうか？ なぜ１６進数なのだろう？
気になったので、調べてみました。

BBPアルゴリズムについては、以下に解説があります。
★ The BBP Algorithm for Pi (David H. Bailey) 8 Sept 2006
http://www.experimentalmath.info/bbp-codes/bbp-alg.pdf
* Pi Directory
>> http://crd.lbl.gov/~dhbailey/pi/
* BBP Code Directory -- C言語とFORTRANのプログラムコードが置かれている。
>> http://www.experimentalmath.info/bbp-codes/
このページには日本語の解説と、十進BASICによるプログラムコードが掲載されています。
* こつこつアルゴリズム -- BBPアルゴリズム
>> http://www14.ocn.ne.jp/~kk62526/pi/BBP.html
このブログにも解説がありました。
* サイエンスものおきば -- 円周率の任意の桁を一発で計算する。
>> http://www.riverplus.net/sci/2005/05/post_83.html
これらの資料、特に最初の★印のPDFを元に“円周率一発計算”の秘密を探りました。

“一発計算”の仕組みを理解するため、まずは円周率ではなく、もう少し分かりやすい数で考えてみます。
★のPDFでは、log 2 (２の自然対数) を使っての説明がありました。
log 2 という数は、次のような足し算を繰り返すことによって計算できます。

※なぜこれで log 2 が計算できるのか、ということについては末尾のメモを参照。
この log 2 という数を２進数で表したときの、小数点以下ｄ桁目の数字を計算することを考えます。
log 2 を 2^d 倍すると、２進数でｄ桁だけ、桁数が上がります。
（１０進数なら 10^d 倍すれば、ｄ桁だけ桁数が上がるよね。それと同じこと。）
log 2 の、小数点以下ｄ桁目の数字を計算するということは、
2^d・log 2 の、小数点以下１桁目を計算するのと同じことです。
ここで、ある数字の少数部分を｛ある数字｝というカッコ記号で書くことにしましょう。
今知りたい小数点以下の数字｛2^d・log 2｝を、最初の足し算の繰り返しに当てはめると、こんな風になります。

この式は、ｄ桁ずらした log 2 の繰り返し足し算を、
最初からｄ番目までの足し算と、ｄ番目より後の足し算に分けて書いたものです。
前半の、最初からｄ番目までの足し算には｛｝が付いています。
それというのも、今、知りたいのは小数点以下の部分だけですから、
整数の部分はバッサリ切り捨てても構わないわけです。
・・・整数の部分は切り捨てても構わない、ということは、前半の部分については
「割り算の商はいらない、余りだけを問題にすればよい」というわけです。
それが式の中に「mod k」という記号を入れた理由です。（赤丸を付けた部分）
（「mod k」というのは、ｋで割ったときの余り、という意味です。）
実はこの「余りだけが分かればよい」というところが、一発計算の秘密なのです。

「余りだけが分かればよい」という計算部分を取り出してみると、つまり
　　(余り) = x ^ d mod k
という計算が高速にできさえすれば、当初の目的が達成できることになるでしょう。
普通に考えると、この計算は xを d回かけ算して、出てきた答を k で割って余りを求める、という手順を踏むことになります。
しかし、d が大きな数になると、d回のかけ算にたいへんな手間がかかります。
そこで、d回のかけ算を、d = (2*2*2*2*・・・) * c といった具合に２のべき乗に分解して、計算ステップを減らすことを考えます。
たとえば１０回のかけ算を、こんな風に２のべき乗に置き換えてみると、、、
・改良前
　　x * x * x * x * x * x * x * x * x * x 　　　（１０回かけ算するところ）
・改良後
　　((((x^2)^2)^2) * x * x　　　（５回のかけ算で済む）
こうやってかけ算の回数が減らせるわけです。
　(x^2)^2)^2)・・・という計算方法によって、x をかける回数は２倍、２倍と指数的に増加してゆきます。
なので、とてつもなく大きな回数のかけ算であっても、この方法によって大幅に計算ステップを減らすことができるわけです。
もう１つ、計算を簡略化するポイントがあります。
x ^ d を先に全部計算して、最後に k で割った余りを求めるのでは、途中に出てくる x ^ d が非常に大きな数になってしまいます。
そこで、最後に k で割るのではなく、かけ算の途中の段階で先に k で割った余りを出しておけば、扱う数が小さくて済みます。
・さらなる改良後
　　((((x^2 mod k)^2 mod k)^2 mod k) * x mod k * x mod k　　　（途中で出てくる数は k より小さくて済む）
以上の考え方をもとにした、x ^ d mod k を高速に計算する方法は「冪乗法」「バイナリ法」などと呼ばれています。
PHPで書くと、こんな風になります。

function BinModExp( $x, $d, $k ){

　　$r = 1;

　　$t = 1;

　　while( $t <= $d ){

　　　　$t = $t * 2;

　　}

　　$t = $t / 2;　　// t に d を超えない最大の 2のべき乗となる数をセットする.

　　

　　while(true){

　　　　if( $d >= $t ){

　　　　　　$r = ($x * $r) % $k;

　　　　　　$d = $d - $t;

　　　　}

　　　　$t = $t / 2;

　　　　if( $t >= 1 ){

　　　　　　$r = ($r * $r) % $k;

　　　　}

　　　　else{

　　　　　　break;

　　　　}

　　}

　　return $r;

}

ここにわかりやすい解説がありました。
* 冪乗法
>> http://www2.cc.niigata-u.ac.jp/~takeuchi/tbasic/BackGround/power.html
* ウィキペディアにも載っていた >> wikipedia:冪剰余

「冪乗法」は、桁数 dが非常に大きくなっても高速という優れた計算方法です。
一発計算の仕組みは、つまるところこの冪乗法によって立っているのです。
もう一度、上の {2^d・log 2} の式を見てみましょう。
赤丸の付いた 2^d-k mod k という部分は、冪乗法によって極めて高速に計算できます。
それでも log 2 の計算には、冪乗法の計算を d回行って足し合わせる必要があるのですが、
実は１回１回の足し算の結果は 1 を超えることがありません。
なぜなら、k で割った余り(mod k) を、さらに分母 k で割っているからです。
d回以降の残りの項は、理屈の上では無限に足し合わせを行っていますが、
いま目的とする数字のためには、必要な精度の項だけを足し合わせれば十分です。
かくして、高速な冪乗法の結果をd回とあともう少し足し合わせることによって、
log 2 の小数点以下特定の数字が計算できるわけなのです。

さて、円周率の特定の数字を計算する方法は、原理的には上の log 2 と同じです。
円周率を計算する方法として、こんな公式が知られています。

なぜこれで円周率が計算できるのかについては、以下のページの「BBP公式の証明」ってところを見てください。(他力本願)
* こつこつアルゴリズム -- BBPアルゴリズム
>> http://www14.ocn.ne.jp/~kk62526/pi/BBP.html
この式の計算に、log 2 と同じように x ^ d mod k の高速計算を適用することを考えます。
計算結果を 16^d 倍すると、それは１６進数 ｄ桁だけ桁数を上げたのと同じことです。
そこで、上の公式を 16^d倍して、その｛小数部分｝の計算を考えてみると・・・

こんな風に、うまいこと冪乗法が適用できるではありませんか！

ここまでわかってしまえば、後はプログラムするだけ。
上のリンクにあったＣ言語版のプログラムを元に、Silverlightに移植してみました。
* BBPアルゴリズムによる円周率の計算 (要Silverlight4)
>> http://brownian.motion.ne.jp/memo/PiBBP/
求めたい桁数を入力し、[計算]ボタンを押してみてください。
桁数は 10^6 = 1000000 程度までは問題なく計算できました。
桁数 10^7 = 10000000 を計算してみると、私の環境(Core i5)では45秒程度で
　7AF5863EFF
という結果になりました。PDF資料にあった正しい値は
　7AF5863EFE
だったので、最後の F <- E が違っていました。
桁数 10^8 = 100000000 は５分程度で
　CB840E2190
正しい値は
　CB840E2192
なので、やはり最後の 0 <- 2 が違っていました。
つまりこのプログラムで正しく計算できたのは 10^6 程度までということです。
プログラム・ソースを公開しておきます。(C#)
>> http://brownian.motion.ne.jp/memo/PiBBP/MainPage.xaml.cs

円周率一発計算の秘密は、冪乗法にありました。
一発と言ってもよく見ると d回の繰り返しがあるため、どんな桁数でも全く同じように計算できるわけではありません。
それでも d回の繰り返しはdに対して一次ですから、まあ許せる程度の時間内に望みの桁数が計算できるわけです。
私は今回、10^8 程度のところでひとまず断念しましたが、チャレンジャーの方、さらに上の桁を目指してはいかがでしょうか？

※ 参考: Log 2 = Σ 1 / (k 2^k) の公式メモ