MPFUN90 で pi 千桁

コンピュータによる多倍長計算を手段とする実験数学を提唱して、最近気炎を上げている D. H. Bailey という米国人がおります。そのべ氏の製作による多倍長計算パッケージを MPFUN90 といいます。High-Precision Software Directory

このパッケージは、Fortran90 で書かれておりまして Fortran90 の演算子のオーバーロード機能等を使用して、かなり自然な形で多倍長計算プログラムが書けるようになっています。またパッケージ自体は FORTRAN77 時代から開発が続いていて、信頼度の高いものとなっております。

最近流行の PSLQ アルゴリズムも、べ氏自身によるこのパッケージの成果のひとつです。

*

ここでは、MPFUN90 を利用して Pi を千桁まで求めてみました。用いたアルゴリズムはガウスの算術幾何平均によるものです。１２反復で値が収束しています。プログラム中 mppic というのは、パッケージにあらかじめ定義されている Pi の値です。

ソース

PROGRAM gauss_pi
    USE mpmodule
    IMPLICIT NONE
    INTEGER, PARAMETER :: n = 25
    INTEGER :: i, k
    TYPE (mp_real) :: a(n), b(n), c(n), s, pi
   
    CALL mpinit(1000)
    k = 12
    MPOUD = 1000
    A(1) = '1.0'
    B(1) = 1 / SQRT( MPREAL('2.0 e0') )
    C(1) = B(1)

    DO i = 2, k
     A(I) = 0.5d0 * ( A(i - 1) + B(i - 1) )
     B(I) =     SQRT( A(i - 1) * B(i - 1) )   
     C(I) = 0.5d0 * ( A(i - 1) - B(i - 1) )
    END DO
    
    s = 0.0d0
    DO i = 1, k
     s = s + 2**(i - 1) * C(i) * C(i)
    END DO
    pi = 2.0d0 * a(k) * a(k) / ( 1.0d0 - s )
    
    PRINT *, 'Pi ='
    CALL mpwrite (6, pi)
    PRINT *
    CALL mpwrite (6, mppic)
    
    STOP
END PROGRAM gauss_pi

実行結果