マチンの公式で円周率を100万桁まで求めるプログラムを書いてみました。

マチンの公式

円周率を求める公式の一つに、マチンの公式というものがあります。
マチンの公式は、以下の式で表すことができます。

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

これを変形すると、

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

となります。

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

C#で実装するには

arccotは、

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

ですから、当然のことながら、実数の計算が必要になるわけですが、ここでは、整数部が一桁で足りることを使って、.NET FrameworkのBigIntegerクラスを使ってarccotを求めています。

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

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

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

なお、最下位の桁では誤差が生じますので、100万桁よりも余分に計算しています。

C#のコード

以下そのコードです。

using System;
using System.Numerics;
using System.Text;

namespace PICalc {
    class Program {
        static void Main(string[] args) {
            string pi = PICalculator.Calculate(1000000);
            Console.WriteLine(PICalculator.Format(pi));
        }
    }

    // 円周率の計算
    static class PICalculator {

        // 小数部figure桁までの円周率を求める
        public static string Calculate(int figure) {
            // マチンの公式により、円周率を求める
            BigInteger a = ArcCot(5, figure) * 4;
            BigInteger b = ArcCot(239, figure);
            BigInteger pi = (a - b) * 4;
            return pi.ToString().Substring(0, figure + 1);   // +1 は、整数部 表示用
        }

        // Calculateで求めた文字列を整形する。
        public static string Format(string s) {
            StringBuilder sb = new StringBuilder();
            sb.Append(s.Substring(0, 1));
            sb.AppendLine(".");

            int j = 0;
            for (int i = 1; i < s.Length; i += 10) {
                sb.Append(s.Substring(i, Math.Min(10, s.Length - i)));
                sb.Append(" ");
                if (++j % 10 == 0)
                    sb.AppendLine();
            }
            return sb.ToString();
        }

        public static BigInteger ArcCot(int m, int figure) {
            var one = BigInteger.Pow(10, figure + 9);
            var mm = m * m;
            int sign = -1;
            BigInteger item = one / m;   // (1/m)
            BigInteger at = item;

            for (int step = 3; ; step += 2) {
                item = item / mm;
                var a = item / step;
                if (a == 0)
                    break;
                if (sign == 1)
                    at = at + a;
                else
                    at = at - a;
                sign = sign * -1;  // signをトグルさせる
            }
            return at;
        }
    }
}

実行してみる

さすがに、100万桁だとすぐに終わってくれません。時間がかかります。
気長に終わるまでまってください。正確に測ってなかったのですが、僕の環境だと30分以上かかりました。

ただ、コンソールにそのまま表示してもコンソールのバッファーを超えてしまうので、以下のようにファイルにリダイレクトしました。

PICalc > piresult.txt

実行した結果をここに載せるには、100万桁はあまりにも大きいので、最初の3000桁だけ掲載します。

3.
1415926535 8979323846 2643383279 5028841971 6939937510 5820974944 5923078164 0628620899 8628034825 3421170679 
8214808651 3282306647 0938446095 5058223172 5359408128 4811174502 8410270193 8521105559 6446229489 5493038196 
4428810975 6659334461 2847564823 3786783165 2712019091 4564856692 3460348610 4543266482 1339360726 0249141273 
7245870066 0631558817 4881520920 9628292540 9171536436 7892590360 0113305305 4882046652 1384146951 9415116094 
3305727036 5759591953 0921861173 8193261179 3105118548 0744623799 6274956735 1885752724 8912279381 8301194912 
9833673362 4406566430 8602139494 6395224737 1907021798 6094370277 0539217176 2931767523 8467481846 7669405132 
0005681271 4526356082 7785771342 7577896091 7363717872 1468440901 2249534301 4654958537 1050792279 6892589235 
4201995611 2129021960 8640344181 5981362977 4771309960 5187072113 4999999837 2978049951 0597317328 1609631859 
5024459455 3469083026 4252230825 3344685035 2619311881 7101000313 7838752886 5875332083 8142061717 7669147303 
5982534904 2875546873 1159562863 8823537875 9375195778 1857780532 1712268066 1300192787 6611195909 2164201989 
3809525720 1065485863 2788659361 5338182796 8230301952 0353018529 6899577362 2599413891 2497217752 8347913151 
5574857242 4541506959 5082953311 6861727855 8890750983 8175463746 4939319255 0604009277 0167113900 9848824012 
8583616035 6370766010 4710181942 9555961989 4676783744 9448255379 7747268471 0404753464 6208046684 2590694912 
9331367702 8989152104 7521620569 6602405803 8150193511 2533824300 3558764024 7496473263 9141992726 0426992279 
6782354781 6360093417 2164121992 4586315030 2861829745 5570674983 8505494588 5869269956 9092721079 7509302955 
3211653449 8720275596 0236480665 4991198818 3479775356 6369807426 5425278625 5181841757 4672890977 7727938000 
8164706001 6145249192 1732172147 7235014144 1973568548 1613611573 5255213347 5741849468 4385233239 0739414333 
4547762416 8625189835 6948556209 9219222184 2725502542 5688767179 0494601653 4668049886 2723279178 6085784383 
8279679766 8145410095 3883786360 9506800642 2512520511 7392984896 0841284886 2694560424 1965285022 2106611863 
0674427862 2039194945 0471237137 8696095636 4371917287 4677646575 7396241389 0865832645 9958133904 7802759009 
9465764078 9512694683 9835259570 9825822620 5224894077 2671947826 8482601476 9909026401 3639443745 5305068203 
4962524517 4939965143 1429809190 6592509372 2169646151 5709858387 4105978859 5977297549 8930161753 9284681382 
6868386894 2774155991 8559252459 5395943104 9972524680 8459872736 4469584865 3836736222 6260991246 0805124388 
4390451244 1365497627 8079771569 1435997700 1296160894 4169486855 5848406353 4220722258 2848864815 8456028506 
0168427394 5226746767 8895252138 5225499546 6672782398 6456596116 3548862305 7745649803 5593634568 1743241125 
1507606947 9451096596 0940252288 7971089314 5669136867 2287489405 6010150330 8617928680 9208747609 1782493858 
9009714909 6759852613 6554978189 3129784821 6829989487 2265880485 7564014270 4775551323 7964145152 3746234364 
5428584447 9526586782 1051141354 7357395231 1342716610 2135969536 2314429524 8493718711 0145765403 5902799344 
0374200731 0578539062 1983874478 0847848968 3321445713 8687519435 0643021845 3191048481 0053706146 8067491927 
8191197939 9520614196 6342875444 0643745123 7181921799 9839101591 9561814675 1426912397 4894090718 6494231961 

以下省略