前回 boost::numeric::interval(1): C++ Boost 区間演算ライブラリ 四則演算と基本的な関数
に続き Boost 区間演算ライブラリの紹介です。
今日はいくつかの区間用の数学関数を説明します。
min, max, abs
区間の関数は基本的にベースの型の関数を拡張したものになっています。これらは、区間の四則演算と同様「区間内の数に対して関数を適用した結果が 全て 含まれているような区間」を計算するものです。
関数 min は 2 つの区間 [a, b], [c, d] を引数にとって、区間 [min(a, c), min(b, d)] を返します。max も同様で、こちらは区間 [max(a, c), max(b, d)] を返します。例えば、max([1, 3], [2, 4]) は [2, 4] になります（4 や [3, 4] ではない）。
一方、関数 abs が返す区間は abs([a, b]) = [abs(a), abs(b)] では ありません。例えば、abs([-3, 1]) は [0, 3] になります。集合の表記法を使えば、abs([a, b]) = {abs(x) | x≦b} と書くことができます。具体的に書き下すには与えられた区間の下端 a と上端 b の正負と絶対値の大小によって場合分けすることになりますが。
# include <iostream>
# include <boost/numeric/interval.hpp>
# include <boost/numeric/interval/io.hpp>
int main () {
typedef double R;
typedef boost::numeric::interval<R> IR;
using std::cout;
using std::endl;
const IR a(-1, 3);
const IR b(2, 4);
cout << "a = " << a << endl;
cout << "b = " << b << endl;
cout << "min(a, b) = " << min(a, b) << endl;
cout << "max(a, b) = " << max(a, b) << endl;
cout << "abs(a) = " << abs(a) << endl;
return 0;
}
a = [-1,3]
b = [2,4]
min(a, b) = [-1,3]
max(a, b) = [2,4]
abs(a) = [0,3]
square, pow, nth_root
これらの関数も同様ですが、いくつか注意が必要です。
まず、square はその名前の通り、区間の数を平方（2乗）した区間を計算しますが、square(x) と x * x は一般には等しくありません。例えば、square([-3, 1]) は [0, 9] ですが、[-3, 1] * [-3, 1] は [-3, 9] となり負の数も含まれてしまいます。
これは、指数関数 pow についても同様です。与えられた指数 n が偶数の場合には pow(x, n) は非負の区間になります。また一般に、単純に n 回乗算を行うよりも pow を使った方が乗算の回数が少ない分、誤差による区間の広がりを抑えられます。なお指数を 0 にすると、一点区間 [1, 1] が返ります（与えた区間が 0 を含んでいてもそうなるので注意）。
最後に、nth_root は n 乗根を計算する関数ですが、この関数が計算する区間は n によって次のように場合分けできます:
- n が正の偶数の場合:
nth_root(pow(x, n), n)⊇abs(x)をみたすような区間を返す。 - n が正の奇数の場合:
nth_root(pow(x, n), n)⊇xをみたすような区間を返す。 - n が非正整数の場合: 定義されていない。
# include <iostream>
# include <boost/numeric/interval.hpp>
# include <boost/numeric/interval/io.hpp>
int main () {
typedef double R;
typedef boost::numeric::interval<R> IR;
using std::cout;
using std::endl;
cout << "--- interval power functions ---" << endl;
const IR a(-1, 3);
cout << "a = " << a << endl;
cout << "square(a) = " << square(a) << endl;
cout << "pow(a, 3) = " << pow(a, 3) << endl;
const IR b(2, 4);
cout << "b = " << b << endl;
cout << "nth_root(b, 2) = " << nth_root(b, 2) << endl;
cout << "--- checking inclusion property ---" << endl;
const IR c = IR(-6, 8) / 3.0;
cout << "c = " << c << endl;
cout << " k \t d = c^k \t e = d^(1/k) \t inclusion" << endl;
for ( int k = 1; k <= 5; k++ ) {
const IR d = pow(c, k);
const IR e = nth_root(d, k);
cout << k << '\t' << d << '\t' << e << '\t';
if (k % 2 == 1 && subset(c, e)) {
cout << " c is a subset of e." << endl;
} else if (k % 2 == 0 && subset(abs(c), e)) {
cout << "abs(c) is a subset of e." << endl;
} else {
cout << endl;
}
}
return 0;
}
--- interval power functions ---
a = [-1,3]
square(a) = [0,9]
pow(a, 3) = [-1,27]
b = [2,4]
nth_root(b, 2) = [1.41421,2]
--- checking inclusion property ---
c = [-2,2.66667]
k d = c^k e = d^(1/k) inclusion
1 [-2,2.66667] [-2,2.66667] c is a subset of e.
2 [0,7.11111] [0,2.66667] abs(c) is a subset of e.
3 [-8,18.963] [-2,2.66667] c is a subset of e.
4 [0,50.5679] [0,2.66667] abs(c) is a subset of e.
5 [-32,134.848] [-2,2.66667] c is a subset of e.
multiplicative_inverse, division_part1, division_part2
割り算に関しては、operator / や div 以外にもより細かい関数があります。
まず、multiplicative_inverse は逆数 1 / x を計算するための関数です。これについては特に説明することもないでしょう。
さて、前回の四則演算の説明の時にきちんと言及していませんでしたが、Boost の区間演算ライブラリで 0 を含む区間による除算 / を行おうとすると、±∞ が絡んできます。例えば、[1, 2] / [-3, 4] は (-∞, ∞) になります。
Remark. boost/numeric/interval/io.hpp の定義は単純なものなので、無限大を含む区間も [-inf,inf] と同じフォーマットで出力されてしまうので注意してください。
しかし、{ x / y | x∈[1, 2], y∈[-3, 4], y≠0 } を包む最小の閉集合は、実際には全区間ではなく (-∞, -1/3] ∪ [1/4, ∞) という 2 つの交わらない閉区間の和です。これは一つの区間では表せないため、Boost の区間演算 / は最初のように全区間を返す訳です。
そこで、このような場合のために、除算の結果をより小さい区間で得るための関数が division_part1 と division_part2 です。part1 と part2 はそれぞれ、除算の結果が分割される時の負の側の区間と正の側の区間に対応します。
-
division_part1(a, b, parted): a を b で割った区間の負の側を返そうとする。- 区間が分割される場合は参照変数 parted に true がセットされ、分割された負の側の区間が返る。
- 区間が分割されない場合は参照変数 parted に false がセットされ、そのまま a / b と同じ結果が返る。
-
division_part2(a, b, parted): a を b で割った区間の正の側を返そうとする。- parted が省略された場合は true とみなされる。
- parted が true の場合、a を b で割った正の側の区間を返す。
- parted が false の場合や、区間が分割されない場合この関数の結果は undetermined である。(エラーが投げられるなど)
# include <iostream>
# include <boost/numeric/interval.hpp>
# include <boost/numeric/interval/io.hpp>
int main () {
typedef double R;
typedef boost::numeric::interval<R> IR;
using std::cout;
using std::endl;
const IR a(1, 2);
const IR b(-3, 4);
cout << "a = " << a << endl;
cout << "b = " << b << endl;
cout << "a / b = " << a / b << endl;
bool parted = false;
const IR c = division_part1(a, b, parted);
const IR d = division_part2(a, b, parted);
if (parted) {
cout << "negative part of a / b = " << c << endl;
cout << "positive part of a / b = " << d << endl;
}
return 0;
}
a = [1,2]
b = [-3,4]
a / b = [-inf,inf]
negative part of a / b = [-inf,-0.333333]
positive part of a / b = [0.25,inf]
今日もこの辺で。本家より長い説明書いてる気がするのはきっと気のせいです。(自重しないと)
今回のサンプルプログラムと出力例: https://gist.github.com/2954837
次: boost::numeric::interval<int>(3): C++ Boost 区間演算ライブラリ 超越関数