matsukawar’s blog

個人的な技術ブログです。SAPネタを充実したい。Twitter : https://twitter.com/matsukawar

.Net 4.0で実装されたBigIntegerクラスを使って円周率を計算する

たぶん、あまり話題になっていない?みたいなのですが、.NetFramework4.0から正式搭載されたSystem.Numerics.BigIntegerクラスで遊んでみたいと思います。前回のバージョンでは搭載が断念された機能でいよいよ搭載された科学技術演算用を想定した仕組みです。

Javaには以前からこのような仕組みがあったみたいなのですが、.Netにはいままで無く、多倍長浮動小数点演算や多倍長整数演算のような大規模な演算は、自力で計算処理を実装する必要がありました。独自に実装したほうが、早いと仰る人もいると思いますが、数学が得意じゃない私からすればありがたい機能です。

BigIntegerクラスでは、四則演算はもちろんのこと、比較演算子がoverrideされてるため、使い方がめちゃくちゃ簡単です。これで、科学技術演算の分野でも.Netが活躍する機会が整った?といえるのではないでしょうか??笑

BigInteger one = new BigInteger(1);
BigInteger big =  BigInteger.Parse("12345678912345678912345678912345678912345678912345・・・・");

big += one;

上の通り、
既存のBIT数で表現しきれない膨大な数値を取り扱うことができます。
通常ですと、規定のbit数(16, 32, 64)を超える数値を計算しようとすると、オーバーフローしてしまいますが、BigIntegerを使うことで、これは回避できます。

たとえば、以下は、ニュートン法による平方根を求める関数です。
※残念ながら、平方根または、逆平方根については、BigIntegerクラスでは標準装備されていないみたいなので、自力でやる必要がありそうです。

/// <summary>
/// ニュートン法による平方根の算出
/// </summary>
/// <param name="x">入力値</param>
/// <param name="x0">初期値</param>
/// <returns>平方根(X)</returns>
/// <remarks>X0として近似値を与えると高速で算出できる</remarks>
private BigInteger SqrtNewton(BigInteger x, BigInteger x0)
{
    BigInteger xnn = x0;                //X(n+1)
    BigInteger xn = new BigInteger(0);  //X(n)
    while (xnn != xn)
    {
        xn = xnn;
        xnn = (xnn + x / xnn) >> 1;
    }
    return xnn;
}

いろいろ調べたところ、ニュートン法は逆平方根の場合に有利で、1/root(X)を算出してから、Xを掛けてroot(X)を求める方法が一般的なようです。ニュートン法の欠点は、割り算を多様する点です。割り算は四則演算の中でも計算量が多く、時間がかかります。この解決策としては、初期値X(0)に、近似値を入れることにより、より早く収束させることです。X(0)の値の求め方についてはいろいろ工夫がありそうです。

本題ですが、円周率を計算してみましたので、
プログラムを紹介します。

アルゴリズムは、サラミン・ブレント法を使用しました。

//桁数
int iSize = Convert.ToInt32(<円周率の桁数>);

//各種変数の初期化            
BigInteger oOne = new BigInteger(1);
for (int i = 0; i < iSize; i++)
{
    oOne = oOne * 10;
}

BigInteger x = oOne;
BigInteger y = oOne >> 1;
BigInteger z = y;
BigInteger sq = 0;

//サラミンブレント法でPIを算出
int iLoop = (int)Math.Log((double)iSize, 2.0) + 1;
for (int i = 1; i < iLoop; i++)
{
    BigInteger pb = (x + y) >> 1;

    //ニュートン法による平方根の算出
    sq = this.SqrtNewton(x * y, pb);

    x = (pb + sq) >> 1;
    y = sq;
    z = z - ((x - sq) << i);
}
BigInteger oPi = ((x + sq) * oOne) / z;

実行結果

左が計算にかかった時間、右が表示にかかった時間です
CPUはCorei7を使用してますorz

たとえば、BOINCのような科学技術計算の体系を.Netで実現してみたり、
計算処理のマルチスレッド化、分散コンピューティングもできるでしょうから(行列演算とか)、
計算処理をWindows Azure上でやるようなツールなんかおもしろそうです。

ref

MSDN : BigInteger 構造体
http://msdn.microsoft.com/ja-jp/library/system.numerics.biginteger.aspx

津留和生先生
アセンブラを使わないC++による多倍長計算ライブラリ
http://hp.vector.co.jp/authors/VA018507/japanese.html