MPFR

mpfrを利用した数値計算例

$$\sum_{k=1}^{n}\frac{1}{k(k+1)}= \frac{n}{n+1}$$

を数値計算し、計算精度を確認する。

プログラム sum.cxx

  1    #include $<$stdio.h$>$
  2    #include $<$mpfrxx.h$>$
  3    
  4    #define LAST 1000
  5    #define PREC 166
  6    
  7    int main(void) {
  8       mpfr_set_default_rounding_mode(GMP_RNDN);
  9       mpfr_set_default_prec(PREC);
 10       int k;
 11       mpfr_class s;
 12       mpfr_class a;
 13       s = 0;
 14       for(k=1; k$\leq$LAST; k++) {
 15          a = 1;
 16          a /= k$\ast$(k+1);
 17          s += a;
 18       }
 19       a = LAST;
 20       a /= (LAST+1);
 21       cout $\ll$ "exact:$\backslash$t" $\ll$ a $\ll$ endl;
 22       cout $\ll$ "sum:$\backslash$t" $\ll$ s $\ll$ endl;
 23       cout $\ll$ "error:$\backslash$t" $\ll$ a-s $\ll$ endl;
 24    }

ここでは、記述の容易さを考慮して、C++言語でプログラムを作成した。

• コンパイルとリンク
  g++ -o sum sum.cxx -lgmp -lmpfr

• 実行と実行結果
  ./sum
  

        exact: 0.999000999000999000999000999000999000999000999000998e0
        sum:   0.999000999000999000999000999000999000999000999000982e0
        error: 0.160365882605531738768432188029088912728911812806303e-49
  

計算精度として、10進数50桁程度を要求すると、

$$50\log_{2}10 = 166.096\cdots$$

であるから、166ビットの多倍長数値計算を行う。

$n=1000$までの総和

$$\sum_{k=1}^{1000}\frac{1}{k(k+1)}$$

と $$\left.\frac{n}{n+1}\right\vert _{n=1000}$$を直接計算した値を 比較すると、実行結果から誤差は

$${\tt error}= 0.16\cdots 10^{-49}$$

であった。

計算順序を入れ換えるなどの工夫をしなくても、簡単に高精度の計算結果を 得ることができる。 また、丸めの方法を変えることで精度保証計算を容易に行うことができる。

GMP:   GNU multiple-precision floating-point computations with exact rounding
MPFR:   multiple-precision floating-roint