boost::numeric::interval<int>(3): C++ Boost 区間演算ライブラリ 超越関数

大分間が空いてしまいましたが、前回の boost::numeric::interval<int>(2): C++ Boost 区間演算ライブラリ 数学関数 に続き Boost 区間演算ライブラリの紹介です。

さて、前回「数学関数」を紹介したあとだいぶ放置した訳ですが… おや何ですか？数学関数を扱うとか言いながら三角関数すら見ていない？はい、そうですね〜 じゃぁ今日は三角関数のような「超越関数」を扱う方法について見ていきましょう。

default_policies

実は、通常推奨されている boost::numeric::interval<double> 型をそのまま使っただけでは、sin, cos などの三角関数や exp などの超越関数（以下、まとめて超越関数と呼びます）は使うことができないようになっています。

これは、boost::numeric::interval<double> におけるこれら区間値超越関数の実装が、rounding mode を切り替えつつ cmath の超越関数を呼び出すだけの単純なものになっているからです。このような実装では、区間演算の包含単調性が保証されません。即ち、区間値関数の計算結果が、（数学的に厳密な） 真の値を含んでいる保証がありません！

こういった理由により、boost::numeric::interval<double, Policies> の 2 つ目のテンプレート引数を省略して default_policies をそのまま使うと、これらの関数が使えないようになっています。使おうとしてもエラーが出ます。

test-interval-transc-error.cpp

#include <iostream>
#include <boost/numeric/interval.hpp>
#include <boost/numeric/interval/io.hpp>

int main () {
    using std::cout;
    using std::endl;
    using namespace boost::numeric;
    typedef double R;
    typedef interval<R> IR;
    const IR x = 1;
    cout << "exp(1) = " << exp(x) << endl;
    return 0;
}

$ g++ test-interval-transc-error.cpp
In file included from /usr/include/boost/numeric/interval.hpp:30:0,
                 from test-interval-transc-error.cpp:2:
/usr/include/boost/numeric/interval/transc.hpp: In function ‘boost::numeric::interval<T, Policies> boost::numeric::exp(const boost::numeric::interval<T, Policies>&) [with T = double, Policies = boost::numeric::interval_lib::policies<boost::numeric::interval_lib::rounded_math<double>, boost::numeric::interval_lib::checking_strict<double> >]’:
test-interval-transc-error.cpp:12:30:   instantiated from here
/usr/include/boost/numeric/interval/transc.hpp:34:64: error: ‘boost::numeric::interval_lib::policies<boost::numeric::interval_lib::rounded_math<double>, boost::numeric::interval_lib::checking_strict<double> >::rounding’ has no member named ‘exp_down’
/usr/include/boost/numeric/interval/transc.hpp:34:64: error: ‘boost::numeric::interval_lib::policies<boost::numeric::interval_lib::rounded_math<double>, boost::numeric::interval_lib::checking_strict<double> >::rounding’ has no member named ‘exp_up’

コンパイルしたらそんな関数は知らんと怒られてますね。

Remark. 他の区間演算ライブラリでは、包含単調性を保証しつつこれらの関数を提供してくれるものもあります（例えば MatLab の IntLab など）。ただ、IntLab の実装を読んでみた限りでは結構打ち切り方とかハードコーディングっぽかったので、自分が必要とする精度で計算したい場合は自分で実装してしまう方がいいのかもしれません。Boost 区間演算ライブラリが標準で超越関数を提供していないのはきっとその辺の思惑もあるのでしょう。

change_rounding

しかし、テンプレート引数 Policies をきちんと指定してやれば、超越関数も 無理矢理 使えるようになります。もちろん、計算の包含単調性は保証されませんから、そこは自己責任でということになりますが。

Remark. 包含単調性が保証されないだけで、おおよそは正しい区間を計算してくれる訳ですから、数値計算以外の用途、例えばゲームの衝突判定に使うと言った用途ならば 無理矢理 使っても何も問題はないでしょう。

さて、じゃぁ Policies に何を書けばいいの？って話ですが、これを直接書くよりも、change_rounding を使った方が（多少）楽です。なので、これを使って超越関数を使えるようにしたサンプルを見てみましょう。

test-interval-transc-change_rounding.cpp

#include <iostream>
#include <boost/numeric/interval.hpp>
#include <boost/numeric/interval/io.hpp>

int main () {
    using std::cout;
    using std::endl;
    using namespace boost::numeric;
    typedef double R;
    typedef interval_lib::change_rounding<
        interval<R>,
        interval_lib::save_state<
            interval_lib::rounded_transc_opp<
                R
            >
        >
    >::type TranscIR;
    const TranscIR x = 1;
    cout << "exp(1) = " << exp(x) << endl;
    return 0;
} 

$ g++ test-interval-transc-change_rounding.cpp -o test-interval-transc-change_rounding.out
$ ./test-interval-transc-change_rounding.out
exp(1) = [2.71828,2.71828]

型の辺りがエラい激しくなりましたがコンパイルも通り期待通り動いています。分かりやすいよう敢えて名前空間を省略してないので激しさマシマシですが、一応、軽く意味を説明していきましょう。

• interval_lib::change_rounding<...>::type
  • 重要なポイントその 1 です。区間型の rounding policies を他のものに変えた新しい区間型を作ってくれます。
  • ::type よく忘れるので注意。
• interval<R>
  • 元になる型です。
• interval_lib::save_state<...>
  • おまじないです。rounding mode save についてはまた後日扱う予定です。
• interval_lib::rounded_transc_opp<...>
  • 重要なポイントその 2 です。"transc" が超越関数、"opp" は rounding mode の切り替え方を表しています。
  • 要するに、超越関数を計算する時に "opp" を使うよってことです。
  • 他にも rounded_transc_std などがあります。
• interval_lib::rounded_arith_opp<R>
  • "arith" が基本的な四則演算を表しています。
  • 要するに、四則演算をする時に "opp" を使うよってことです。
  • デフォルトではこの rounded_arith_opp までしか使われていません。
  • 他にも rounded_arith_std などがあります。

今回はとにかく無理矢理でもいいから使えるようにすることが目的ですから、詳しくは、本家ドキュメントとソースを参照してください。

---

まとめ。

• 正しく精度保証できないため、推奨デフォルトの区間型では超越関数は使えない。
• 正しく精度保証したいなら自分で実装しよう！（頑張ります……）
• 超越関数を使えるようにした新しい区間型は
  • change_rounding<IR, save_state<rounded_transc_opp<R, rounded_arith_opp<R> > > >::type

今回のサンプルプログラムと入出力例: https://gist.github.com/3330726