CLNというのを知ったので自分が現時点で知っているものをまとめておくことに。(CLNについては後述)

MPFRを直接Cで使うというのは旧ブログで扱いましたが(MPFR を使ってみる)、お世辞にも使いやすいとはいえない(四則演算すら面倒くさい)のでPythonから使うという手もあります。

Pythonでやるにはgmp,mpfr,mpc(これらはCのライブラリ)をインストールしてからpipでgmpy2を入れる。(他の方法もあるはずですが公式のWelcome to gmpy2’s documentation!あたりを参照)

コンパイルも必要ないし、各種の演算子はオーバーロードされているので使いやすい。しかし遅い。

C++,JavaならApfloatがある。

ClojureユーザーならjavaのApfloatをClojureのREPLから使うのが便利。拙作のclojure-apfloatを一応宣伝しておきますが、そもそもClojureユーザー自体が少いので需要はあまりなさそう。

今回知ったのはCLNというC++のライブラリで備忘録をかねて、ちょっとご紹介。

参考: CLN

自分が知らなかっただけでCopyrightをみると1995年からあるようなのですでに広く使われているでしょう。

OSのパッケージシステム(apt,homebrewなど)を使えばインストールは簡単。

サンプルコード:

#include <stdint.h>
#include <iostream>
#include <cln/cln.h>

using namespace cln;

float_format_t PREC;
cl_R EPS;

cl_R f(cl_R x)
{
    return x*x - 2;
}

cl_R findroot(cl_R (*f)(cl_R), cl_R x1, cl_R x2)
{
    cl_R a = x1, b = x2;
    cl_R m,v;
    
    while (true) {
        m = (a+b)/2;
        v = f(m);
        if (b - a < EPS) {
            break;
        }
        // 端点での符号判定、積が負かどうか(あるいは signum を使う)
        if (minusp(f(a)*v)) { 
            b = m;
            continue;
        } else if (minusp(v*f(b))) {
            a = m;
            continue;
        } else {
            std::cout << "No root found in [" << x1 << "," << x2 << "]" << std::endl;
            m = 0;
            break;
        }
    }
    
    return m;
}

int32_t main()
{
    PREC = float_format(100);
    // 初期化には整数か文字列を使うべき
    // EPS = "1e-100_100"; // 文字列の最後に精度を指定する _100 の部分
    EPS = expt(cl_float(10,PREC),-100); // 10^(-100) 

    std::cout << findroot(f, cl_float(1,PREC), cl_float(2,PREC)) << std::endl;
    std::cout << sqrt(cl_float(2,PREC)) << std::endl;
 
    return 0;
}

コードは原始的なfindrootで2の平方根を求めて、組み込みのsqrt(2)と比較するもの。ファイルをhoge.cppとすると、コンパイルは以下のようになる(最適化は好みで):

$ g++ -O3 hoge.cpp -o hoge -lcln

実行すると

$ ./hoge
1.4142135623730950488016887242096980785696718753769480731766797379907324784621070388503875343276415727351010794671223L0
1.4142135623730950488016887242096980785696718753769480731766797379907324784621070388503875343276415727350138462309123L0

となる。最後のほうの桁が一致していないが指定した精度(100桁)はある。単に余分な桁まで表示されてしまうようです。

余分な桁を表示したくないので何かフォーマット指定子のようなものがあるのかドキュメントをあたってみたが分からなかったので、C++の文字列として処理する(substr等)しかないのかもしれない。

まとめ

Pythonからgmpy2を使うのが手軽。それで速度が気になる場合はMPFR,Apfloat,CLNあたりを使ってみる。MPFRをCから使うのはかなり面倒なので選択肢はいくつか知っていてもいいと思う。今のところ、この分野ではC++に軍配があがるような気がします。好き嫌いはともかく。