円周率の計算 - Gushwell's C# Programming Page

円周率を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

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

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

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

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

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

value[0] : 5600 value[1] : 1234
と格納されています。

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

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

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

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

■PICalculator.cs

using System;
using System.Text;

namespace CalcPI10000 {
    // 円周率の計算
    // Silverlight 非依存
    static class PICalculator {
        // 小数部figure桁までの円周率を求める
        public static string Calculate(int figure) {
            BigNumber.Length = figure / 4 + 2;    // 最下位の桁には誤差が含まれるので、余分に領域を確保する
            // マチンの公式により、円周率を求める
            BigNumber a = BigNumber.Atan(5).Multi(4);
            BigNumber b = BigNumber.Atan(239);
            BigNumber pi = a.Sub(b).Multi(4);
            //string pi = Regex.Replace(PI, @"\s", "", RegexOptions.Multiline);
            return pi.ToString().Substring(0, figure + 2);   // +2 は、整数部とピリオドの表示用
        }

// Calculateで求めた文字列を、整形する。
        public static string Format(string s) {
            StringBuilder sb = new StringBuilder();
            sb.Append(s.Substring(0, 2));
            int j = 0;
            for (int i = 2; i < s.Length; i += 10) {
                if (j % 10 == 0 && j > 0)
                    sb.Append("  ");
                sb.Append(s.Substring(i, Math.Min(10, s.Length - i)));
                sb.Append(" ");
                if (++j % 10 == 0)
                    sb.AppendLine();
            }
            return sb.ToString();
        }
    }

// 高精度演算を行うクラス 最低限の機能のみ実装 （汎用的には使えない）
    // Silverlight 非依存
    class BigNumber {
        // Length * 4 が有効桁数 
        // 1000は、デフォルト値
        public static int Length = 1000;   
        // value[n] には、4桁の数が格納される。いわゆる10000進数
        // value[0] が最下位の桁、value[Length-1]が最上位の桁
        private int[] value;

// コンストラクタ
        public BigNumber() {
            value = new int[Length];
            for (int i = 0; i < value.Length; i++)
                value[i] = 0;
        }

// 1をセット 
        public void SetOne() {
            // 整数部は最上位の一桁だけであり、それ以外はすべて小数部
            // 例えば、value[Length - 1] = 1234; ならば、1.234 をセットしていることになる。
            value[Length - 1] = 1000;
        }

// ATan(1/m) を求める
        public static BigNumber Atan(int m) {
            BigNumber one = new BigNumber();
            one.SetOne();
            int mm = m * m;
            int sign = -1;
            BigNumber item = one.Divide(m);   // (1/m)
            BigNumber at = item.Multi(1);

for (int step = 3; ; step += 2) {
                item = item.Divide(mm);
                var a = item.Divide(step);
                if (a.IsZero())
                    break;
                if (sign == 1)
                    at = at.Add(a);
                else
                    at = at.Sub(a);
                sign = sign * -1;
            }
            return at;
        }

// ゼロかどうか判断
        private bool IsZero() {
            for (int i = 0; i < Length; i++)
                if (value[i] != 0)
                    return false;
            return true;
        }

// 加算  a + b   ただし (a >= 0 , b >= 0)
        public BigNumber Add(BigNumber num) {
            int[] result = new int[Length];
            int carry = 0;
            for (int i = 0; i < Length; i++) {
                int n = value[i] + num.value[i] + carry;
                if (n >= 10000) {
                    result[i] = n % 10000;
                    carry = n / 10000;
                } else {
                    result[i] = n;
                    carry = 0;
                }
            }
            return new BigNumber() { value = result, };
        }

// 減算 a - b  (ただし a > b > 0)
        public BigNumber Sub(BigNumber num) {
            int[] result = new int[Length];
            int carry = 0;
            int fig = 1;
            for (int i = 0; i < Length; i++) {
                int n = value[i] - num.value[i] - carry;
                if (n < 0) {
                    result[i] = n + 10000;
                    carry = 1;
                } else {
                    result[i] = n;
                    carry = 0;
                    if (n != 0)
                        fig++;
                }
            }
            return new BigNumber() { value = result };
        }

// 乗算  
        public BigNumber Multi(int num) {
            int[] result = new int[Length];
            int carry = 0;
            for (int i = 0; i < Length; i++) {
                int n = value[i] * num + carry;
                if (n >= 10000) {
                    result[i] = n % 10000;
                    carry = n / 10000;
                } else {
                    result[i] = n;
                    carry = 0;
                }
            }
            return new BigNumber() { value = result };
        }

// 除算  
        public BigNumber Divide(int num) {
            int[] result = new int[Length];
            int carry = 0;
            for (int i = Length - 1; i >= 0; i--) {
                int a = value[i] + carry * 10000;
                int c = a / num;
                carry = a - c * num;
                result[i] = c;
            }
            return new BigNumber() { value = result };
        }

public override string ToString() {
            int ix = Length - 1;
            string s = value[ix].ToString("0000");
            s = s.Insert(1, ".");
            for (int i = ix - 1; i >= 0; i--) {
                s += value[i].ToString("0000");
            }
            return s;
        }
    }
}

■MainPage.xaml.cs

■MainPage.xaml

</ScrollViewer>
        <Button Content="計算" Height="23" HorizontalAlignment="Left" Margin="205,14,0,0" Name="button1" VerticalAlignment="Top" Width="75" Click="button1_Click" />
        <TextBox Height="24" Width="50" Margin="70,13,0,0" Name="textBox1" HorizontalAlignment="Left" VerticalAlignment="Top" Text="1000" />
        <TextBlock Height="23" HorizontalAlignment="Left" Margin="124,17,0,0" Name="textBlock2" Text="桁まで求める" VerticalAlignment="Top" Width="78" />
        <TextBlock Height="23" HorizontalAlignment="Left" Margin="19,17,0,0" Name="textBlock3" Text="円周率を" VerticalAlignment="Top" Width="50" />
    </Grid>
</UserControl>

Gushwell's C# Programming Page

円周率を10000桁まで求める

円周率を10000桁まで求める