gmplogo2.png

Introduction

GMP
GNU Multi Precision = 任意精度数演算ライブラリ

Linuxで数値計算をしたいが、精度がどこまで出せるのか、どこまで大きい・小さい数を扱えるのか? 数十bitじゃあ、だめだめ。そこでGMP(任意精度数演算ライブラリ)です。

普通に標準関数を使ってgccでコンパイルしたプログラムの場合整数 int は32bit(符合つきなら31bit、つまり-2147483648から2147483647)、十進で約10桁。 倍精度浮動小数点数 doubleでは64bitの精度。

  • doubleで扱える最大最小値(絶対値)
    e^{-708}\approx      10^{-308} <|x|<10^{308}\approx e^{710}

通常はこの精度でも問題がない。 しかし、多次元(数万とか)の確率分布を考えると、expの方が数千、数万となり「零」になってしまう。

\log \sum_i \exp( -\chi_i^2/2 )

のような場合(\log[ \exp(-\sum_i \chi_i^2/2) ]ではない)、 一度はものすごく小さい数を扱わなければならない。 このような場合「double」では不十分。 そのような時に任意精度数計算ライブラリの出番である。

GMPを使えば、基本的に任意の桁数の計算ができる。 (GMPはC言語 on Linuxで使うことが想定されている?)

私はGMPを直接使うことはしない。MPFR, MPFRCPPを介してGMPを使っている。

詳細については

  • MPFR :MPFRはCで記述された, 任意精度浮動小数点演算のためのポータブルなライブラリで, ベースにはGNU MPライブラリを用いている。 本ライブラリの目的は,GNU MPライブラリの浮動小数点数クラスを精密なやり方で拡張することである。
  • MPFRCPP :MPFRCPPはMPFRのC++ interfaceを提供する。

を参考に。

MPFR

#include <iostream>
#include <iomanip>
#include <mpfr.h>
using namespace std ;

#define N_M 10000

void get_chisq(const int i,mpfr_t &x) ;

int main()
{
  mpfr_t tmp, x, log_tmp ;
  mpfr_init(tmp) ;
  mpfr_init(x) ;
  mpfr_init(log_tmp) ;

  mpfr_set_ld(tmp,0.0,GMP_RNDN) ;
  mpfr_set_ld(x,0.0,GMP_RNDN) ;
  mpfr_set_ld(log_tmp,0.0,GMP_RNDN) ;

  for( int i=0; i<N_M; i++){
    get_chisq( i, x) ;
    mpfr_add( tmp, tmp, x, GMP_RNDN) ;
  }
  mpfr_log( log_tmp, tmp, GMP_RNDN) ;

  double m = N_M ;
  mpfr_t mm ;
  mpfr_init(mm) ;
  mpfr_set_ld( mm, m, GMP_RNDN) ;

  mpfr_t M ;
  mpfr_init(M) ;
  mpfr_log( M, mm, GMP_RNDN) ;

  mpfr_sub( log_tmp, log_tmp, M, GMP_RNDN) ;

  double result = mpfr_get_d( log_tmp, GMP_RNDN) ;
  cout << result << endl ;

  return 0 ;
}

void get_chisq(const int i,mpfr_t &x)
{
  double chisq = 76000.0 + i%100 ;
  chisq *= -0.5 ;

  mpfr_t chisq2 ;
  mpfr_init(chisq2) ;
  mpfr_set_ld( chisq2, chisq, GMP_RNDN) ;

  mpfr_exp( x, chisq2, GMP_RNDN) ;
}

MPFRCPP

#include <iostream>
#include <mpfrcpp/mpfrcpp.hpp>
#include <mpfrcpp/extra/simple_functions.hpp>

using namespace std ;

#define N_M 100000

void get_chisq(const int i,mpfr::Real &x) ;


int main () {
  mpfr::Real tmp, x, log_tmp ;
  for( int i=0; i<N_M; i++){
    get_chisq(i,x) ;
    tmp += x ;

    cout << "x=" << x << ", tmp=" << tmp << endl ;
  }
  log_tmp = mpfr::log(tmp) ;

  double m = N_M ;
  mpfr::Real M(m) ;

  log_tmp -= mpfr::log(M) ;

  double result = log_tmp ;
  cout << result << endl ;

  return 0;
}

void get_chisq(const int i,mpfr::Real &x)
{
  double Tmp = 76000.0 + i%100 ;
  Tmp *= -0.5 ;

  x = mpfr::exp(Tmp) ;
}

添付ファイル: filempfrcpp.cc 348件 [詳細] filetest.cc 346件 [詳細]