« 科学シンポジウムに寄ってみた。 | メイン | 暗号パズル。 »

円周率の任意の桁を一発で計算する。

 先日の記事で紹介した、円周率の任意の桁の数字をもとめる方法について。 
 16進数での表現をもとめるには、以下の公式を紹介しました。
 この公式を使うことで、きわめて小さいメモリ使用量での計算が可能になります。

20050521a
 ・・・ (1)

 ですが、ぱっと見ただけでは、実際にこの式どう使うのかがよく分かりませんでした。
 適当にグーグルしたものの、適切な解説は見つからず。
 そこで、この公式を使った、円周率の任意の桁の導出アルゴリズムについて、以下に解説を提供。
 というか、むしろ、ここのページの日本語訳です。
 
 
【BBP公式による円周率の任意の桁の導出アルゴリズム】

 上記の (1) 式は、4つの部分に分けられます。
 それらを、以下のように、s1, s2, s3, s4 と表記します。
 s1 と s2 については、係数の4と2をひとまず除外しています。
 
20050521b
 ・・・ (2)
 
 いま、円周率πを16進数で表したときの、小数第d位以下の桁の数字が知りたいとします。
 すぐに分かりますが、これは、π に 16d をかけたものの小数部分に等しいです。
 同様に、s1 の小数第d 位以下の桁の数字は、s1 に 16d をかけたものの小数部分に等しいです。
 これは、zの小数部分を、frac(z) と表記すると、frac(s1×16d ) と書けます。
 frac( s1×16d ) は、次のように変形できます。
 
20050521c
  ・・・ (3)
 
 1行目から2行目は、単にk=dで分けただけ。
 2行目から3行目で、余りを求める演算 mod が入っていますが、この式で必要なのは小数部分だけなので、和の結果に影響はありません。

 さて、この式の前半部分の計算をします。
 ここに登場する、べき乗した数から余りをもとめる計算というのは、それほど大変な作業ではないんですね。
 16k の値を計算して格納する必要はなく、1回16をかけるごとに mod 演算をしても答えは変わりないので、メモリ使用量を節約して計算することができます。
 3兆桁程度であれば、値を格納するのに64ビットのメモリがあれば十分だそうです。
 
 こうして、式の前半部分の計算が完了します。
 
 次に、後半部分を計算します。
 後半部分には無限個の項が含まれています。
 ですが、この級数はゼロに収束するので、ある程度のkのところで計算を打ち切ってかまいません。
 ここの計算では、前半部分と違い、べき乗の計算を自力でやらなくてはいけないのですが、この収束のスピードは非常に速いため、それほど大きな計算量というわけではありません。

 以上のようにして、s1×16d の小数部分が求められます。
 同様に s2, s3, s4 についても求め、4*s1-2*s2-s3-s4 を計算すれば、円周率の小数第d位以下の桁の数字がもとめられます。
 
-----------------------------------------------

 以上です。
 さて、やっぱり気になるのは、任意の基数での表現方法ですね。。
 Simon Plouffe により、この方法が見つけられたそうなのですが(ここ)、どういう方法なのかなんとなくは分かったものの、実際にどうやって計算できるのか、未だ不十分な理解で困っています。
  
 と思ったら、上記の方法のコードを書いてる方がいらっしゃいました。ここ Home>今日の一行>2005年5月9日)。
 いったいどういう仕組みなんだろう、これ。そもそも全く見たことない言語だ;

トラックバック

このエントリーのトラックバックURL:
http://www.riverplus.net/cgi/mt/mt-tb.cgi/325

コメント (1)

匿名:

もっと簡単にお願いします

コメントを投稿