MPFR

mpfrを利用した数値計算例

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

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

プログラム sum.cxx

indentation

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 }

• コンパイルとリンク
  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}$$

であった。

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