10.1 MPFR C++ by Pavel Holoborodko

GNU MPFR を使ってみよう、と考えた。 C++ のインターフェイスはどうなっているんだっけ？ (本当は、MPFR 自体が C++ インターフェイス持っていれば良いんだよな…) The GNU MPFR Library を見たら、 百花繚乱になっているらしいことが分かった。 うーん。 まずは ``MPFR C++'' の検索でトップに出て来た MPFR C++ - Pavel Holoborodko をやってみる。

ライセンスは GPL である。と書いてあるけれど、

MPFR C++ is free for usage in free projects. If you intend to use it in commercial application please contact author for permission.

とも書いてある。はてな。

boost で使われていると書いてある (本当？)。

Eigen で使うには、

#include <unsupported/Eigen/MPRealSupport>

だそうな。

とにかくやってみよう。mpfrc++-3.6.2.zip (2015年7月に出来た？) を入手して

mkdir mpfrc++-3.6.2
cd mpfrc++-3.6.2
unzip ../mpfrc++-3.6.2.zip

実にあっさりしている。 実体は、mpreal.h というインクルード・ファイル一つ (3100行くらい)。

MPFR 自体は、MacPorts でインストールされているので、 /opt/local/include, /opt/local/lib に置かれている。

mpreal.h はどこに置くのが良いか。 自分で導入したものなので /usr/local/include が自然か。 ここに置いておくと、LLVM の g++ は無指定で見てくれる。 一方で、MacPorts で入れた GCC の g++ は、 /opt/local/include を無指定で見てくれる。 そこで

sudo cp -p mpreal.h /usr/local/include
sudo ln -s /usr/local/include/mpreal.h /opt/local/include

とした。これで mpreal.h を見るために -I は必要ない。

サンプル・プログラムは唯一 example.cpp がついている (107行)。それはあっさり動いた。

$a_1=1$ , $a_{n+1}=3a_n+2$ ($n\in\N$ ) で定義される 数列の第$100$ 項を計算してみる (すごい勢いで大きくなるので、double では正確に計算できない)。

/*
 * test-mpfrcc++.cpp
 * 準備
 *   (1) http://www.holoborodko.com/pavel/mpfr/ から mpfrc++-3.6.2.zip 入手
 *       unzip mpfrc++-3.6.2.zip として出て来る mpreal.h を適当な場所におく。
 *   (2) MPFR を用意する。MacPorts を使っていると、/opt/local にあるだろう。
 * コンパイル
 *   g++ -I/opt/local/include test-mpfrc++.cpp -L/opt/local/lib -lmpfr
 */

#include <iostream>
#include "mpreal.h"

using namespace std;

int main(void)
{
  using mpfr::mpreal;
  const int digits = 50;
  mpreal::set_default_prec(mpfr::digits2bits(digits));
  mpreal a;
  int i;
  cout.precision(digits);
  a = 1;
  for (i = 2; i <= 100; i++)
    a = 3 * a + 2;
  cout << a << endl;

  return 0;
}

g++-mp-6 test-mpfrc++.cpp -L/opt/local/lib -lmpfr

あるいは

g++-mp-5 -std=c+11 test-mpfrc++.cpp -L/opt/local/lib -lmpfr

あるいは

g++ -I/opt/local/include test-mpfrc++.cpp -L/opt/local/lib -lmpfr

でコンパイルする。 これであっさり動いた (343585013821340887357640753177080848468071681333 という結果が表示される)。

次は、 自然対数の底 $e$ を 1000 桁計算してみよう。 画面に表示する前に文字列にして、50桁ずつわけて表示する。

e1000.cpp

/*
 *   g++ -I/opt/local/include e1000.cpp -L/opt/local/lib -lmpfr
 */

#include <iostream>
#include <string>
#include "mpreal.h"

using namespace std;

int main(void)
{
  int n, maxn = 20;
  using mpfr::mpreal;
  const int digits = 1000 + 10; // 10は余裕
  mpreal::set_default_prec(mpfr::digits2bits(digits));
  mpreal t, s, news;
  string str;
  t = 1; s = t;
  for (n = 1; n < 1000; n++) {
    t /= n;
    news = t + s;
    if (news == s) {
       break;
    }
    s = news;
  }
  cout.precision(digits);
  // cout << "e=" << s << endl;
  str = s.toString();
  cout << "e=" << str.substr(0, 2) << endl;
  for (n = 0; n < digits; n += 50) {
    cout << "    " << str.substr(n + 2, 50) << endl;
  }
  return 0;
}

e1000の実行結果

[chronos:math/new-sotsuken-text/on-multiprecision] mk% ./a.out
e=2.
    71828182845904523536028747135266249775724709369995
    95749669676277240766303535475945713821785251664274
    27466391932003059921817413596629043572900334295260
    59563073813232862794349076323382988075319525101901
    15738341879307021540891499348841675092447614606680
    82264800168477411853742345442437107539077744992069
    55170276183860626133138458300075204493382656029760
    67371132007093287091274437470472306969772093101416
    92836819025515108657463772111252389784425056953696
    77078544996996794686445490598793163688923009879312
    77361782154249992295763514822082698951936680331825
    28869398496465105820939239829488793320362509443117
    30123819706841614039701983767932068328237646480429
    53118023287825098194558153017567173613320698112509
    96181881593041690351598888519345807273866738589422
    87922849989208680582574927961048419844436346324496
    84875602336248270419786232090021609902353043699418
    49146314093431738143640546253152096183690888707016
    76839642437814059271456354906130310720851038375051
    01157477041718986106873969655212671546889570350354
    0212340769
[chronos:math/new-sotsuken-text/on-multiprecision] mk%

Mathematica で N[E,1001] として表示したものと比較して、 一致することを確認出来た。