サイエンスものおきば

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

　･･･　(1)

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

　上記の (1) 式は、４つの部分に分けられます。
　それらを、以下のように、s₁， s₂， s₃， s₄ と表記します。
　s₁ と s₂ については、係数の４と２をひとまず除外しています。
　

　･･･　(2)
　
　いま、円周率πを１６進数で表したときの、小数第ｄ位以下の桁の数字が知りたいとします。
　すぐに分かりますが、これは、π　に 16^dをかけたものの小数部分に等しいです。
　同様に、s₁ の小数第ｄ 位以下の桁の数字は、s₁ に 16^dをかけたものの小数部分に等しいです。
　これは、ｚの小数部分を、frac(ｚ) と表記すると、frac(s₁×16^d ) と書けます。
　frac( s₁×16^d ) は、次のように変形できます。
　

 　･･･　(3)
　
　１行目から２行目は、単にｋ＝ｄで分けただけ。
　２行目から３行目で、余りを求める演算 ｍｏｄ が入っていますが、この式で必要なのは小数部分だけなので、和の結果に影響はありません。

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

　以上のようにして、s₁×16^d の小数部分が求められます。
　同様に s_2， s_3， s₄ についても求め、4*s₁－2*s₂－s₃－s₄ を計算すれば、円周率の小数第ｄ位以下の桁の数字がもとめられます。
　
-----------------------------------------------

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