円周率を10000桁まで求める

マチンの公式で円周率を10000桁まで求めるプログラムです。ページ中ほどで実際に動かしてみることができます。 

マチンの公式を使うと、

π= ( arctan(1/5) - arctan(1/239) ) * 4

 で、円周率が求められます。 これを変形すると、

π= ( arccot(5) - arccot(239) ) * 4	   

となります。
そのため、円周率は、arccotを求める多倍長の演算関数を作れば、求められることになります。
arccotは、

arccot(k) = (1/k) - (1/3)(1/k)^3 + (1/5)(1/k)^5 - (1/7)(1/k)^7 + ... 

ですから、当然のことながら、実数の計算が必要になるわけですが、
ここでは、 整数部が一桁で足りることを使って、多倍長整数演算で、arccotを求めています。

例えば、2/7 を小数点8桁まで求めるとすれば、正数 9桁で以下の ように演算すれば、
2/7 の結果である 0.28571428 が求められることになります。

2/7 = 2 * 100000000 / 7
      = 200000000 / 7
      = 028571428

上位の 1桁目と2桁目の間に、小数点があると仮定すれば、答えの 0.28571428 が得られるわけです。

つまり、多倍長整数による、加減乗除ができれば、円周率を求めることができます。
.NET Framework4 では、BigIntという任意の桁数の演算が行えるクラスが用意されているようですが、
ここでは、.NET Framework3.5以前でも利用できるようにBigNumber というクラスを自作しています。

このBigNumberクラスは、円周率を求めるための最低限の機能のみを実装していますので、
他に流用するには、もう少し汎用性を持たせる必要があると思います。

ソースコードのコメントにも書いていますが、内部データとしては、intの配列として値を 保持しています。


value[] の各要素には4桁の整数を格納するようにしています。
value[0] が最下位の桁、value[Length-1]が最上位の桁です。
たとえば、配列のサイズが 2 のとき、1.23456 という数値は、

 
と格納されています。

なぜ、4桁なのかですが、4桁にすることで、9999 * 9999 < 2^31 になり、
int どうしの演算が可能に なります。
ところが、5桁の場合、99999 * 99999 > n^31 となり、オーバーフローしてしまい、
桁上がりの処理ができなくなってしまいます。
これは、見方を変えると、10000進数を表現しているとも言えます。

なお、下位の桁では誤差が生じますので、余分に配列を確保しています。

Silverlightアプリとして作成しましたので、実際に動かして確かめることができます。
初期値は、1000桁として有りますが、10000桁まで指定できます。


以下、C#+Silverlightのコードです。Solver.csは、Silverlightに依存していません。
XAMLでのEllipseとTextBlockの配置がかなり力技であるのは、笑って許してください。

■PICalculator.cs


■MainPage.xaml.cs

■MainPage.xaml