4 MPFR

MPFR については、 MPFR というサイトがある。

Laurent Fousse, Guillaume Hanrot, Vincent Lefèvre, Patrick Pélissier, ``MPFR: A Multiple-Precision Binary Floating-Point Library With Correct-Rounding'' という INRIA のレポートがある。 その abstract によると、 MPFR は、 C言語で記述され、GMP に基いた、多倍長浮動小数点演算ライブラリィである。 正しい丸めと例外を用意することで、 IEEE 754 標準の考え方を多倍長演算に拡張した点に特徴があるとしている。

また C99 に含まれるすべての数学関数を多倍長でサポートしている (底が $e$ , $2$ , $10$ の対数関数と指数関数と log1p() ( $=\log(1+x)$ ) と expm1() ( $=\exp(x)-1$ ), 三角関数、逆三角函数、双曲線函数、逆双曲線関数、 ガンマ関数、ゼータ関数、誤差関数 (?erfc())、 算術幾何平均、冪乗 pow() ($=x^y$ )。 すべて「正しく丸められる」そうで、 精度保証数値計算家にとってはありがたいだろう …と思ったら既に MPFR に基づく区間演算ライブラリィがいくつかあるそうだ。

使用しているアルゴリズムについては、 THE MPFR TEAM によって ``THE MPFR LIBRARY: ALGORITHMS AND PROOFS'' という文書が用意されている。 ``proofs'' が入っているところが流石。

普段使っている MacBook Air には、気がついたときには既に MPFR が入っていた。 おいおい、どういうことだ (調査中)。 GMP に MPFR が取り込まれたという説もある。 どこかまっさらなマシンに自分でインストールしてみないと分からない。

napier-MPFR.c

/*
 * napier-MPFR.c --- MPFRライブラリィを使って自然対数の底を求める。
 *  http://www.mpfr.org/sample.html
 */

#include <stdio.h>

#include <gmp.h>
#include <mpfr.h>

int main (void)
{
  unsigned int i;
  mpfr_t s, t, u;

  mpfr_init2 (t, 200);
  mpfr_set_d (t, 1.0, GMP_RNDD);
  mpfr_init2 (s, 200);
  mpfr_set_d (s, 1.0, GMP_RNDD);
  mpfr_init2 (u, 200);
  for (i = 1; i <= 100; i++)
    {
      mpfr_mul_ui (t, t, i, GMP_RNDU);
      mpfr_set_d (u, 1.0, GMP_RNDD);
      mpfr_div (u, u, t, GMP_RNDD);
      mpfr_add (s, s, u, GMP_RNDD);
    }
  printf ("Sum is ");
  mpfr_out_str (stdout, 10, 0, s, GMP_RNDD);
  putchar ('\n');
  mpfr_clear (s);
  mpfr_clear (t);
  mpfr_clear (u);
  return 0;
}

MacBook Airでのコンパイルと実行

katsurada:~/study-c-c++-2010 mk$ gcc napier-MPFR.c -lmpfr -lgmp
katsurada:~/study-c-c++-2010 mk$ ./a.out 
Sum is 2.7182818284590452353602874713526624977572470936999595749669131
katsurada:~/study-c-c++-2010 mk$

• mpfr_t 変数名リスト; で変数を定義出来る。
• mpfr_init2(変数名, ビット数) で変数の精度を指定できる。
• mpfr_set_d(変数名, double型の数値, GMP_RNDD) で、 double 型の数値を変数に代入する。_d は double を意味する。
• GMP_RNDN, GMP_RNDZ, GMP_RNDU, GMP_RNDD, GMP_RND_MAX という丸めモードがある。
• mpfr_mul(z,x,y,GMP_RNDD), mpfr_mul_d(z,x,y,GMP_RNDD), mpfr_mul_ui(z,x,y,GMP_RNDD) はいずれも、 z=x*y という意味だが y の型が違う。
• mpfr_out_str(stdout, 基数, 0, 数値, GMP_RNDN)

(2016/5/10記) kv の MPFR ラッパーというのを試してみた (http://nalab.mind.meiji.ac.jp/~mk/labo/text/studying-kv/node15.html)。