前回 boost::numeric::interval<int>(0): C++ Boost 区間演算ライブラリ概要の続きです。今日は区間の四則演算と、よく使う便利な関数について簡単に紹介します。

区間同士の四則演算

前回、区間演算における重要なコンセプトと加法の例を既に書きましたが、他の四則演算についても基本的に同様です。浮動小数点誤差まで考慮した上で区間の数全てに対する演算結果を全て含むような区間を計算して返します。

test-interval-arith-IRxIR.cpp

# 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.0, 2.0);
	const IR b(3.0, 4.0);
	cout << "a = " << a << endl;
	cout << "b = " << b << endl;
	cout << "a + b = " << a + b << endl;
	cout << "a - b = " << a - b << endl;
	cout << "a * b = " << a * b << endl;
	cout << "a / b = " << a / b << endl;
	return 0;
}

また、左辺値に対して四則演算を直接行う演算子 +=, -=, *=, /= もあります。

Remark. 一般に、区間の二項演算 + と - は、互いに逆演算の関係には なっていません 。これは、浮動小数点数型の場合に限りません。例えば、元々誤差がなく、+ と - が互いに逆演算の関係にある int をベースの型とする区間型 interval<int> の場合を考えても、[1, 2] - [1, 2] は [0, 0] にはなりません。* と / についても同様です。このように、一般には区間演算を繰り返すこで区間はどんどん大きくなっていきます。浮動小数点数型の場合には、区間型を考えることで、元々逆演算の関係になかった四則演算のもつ「誤差」という性質がより顕著になったともいえるでしょう。

区間と数の四則演算

四則演算の引数の片方は、ベースの型 T の値で与えることもできます。そうした場合、まず引数を interval<T> に変換して演算を行ったのと同じ結果になります（これは以降出てくる他の関数についてもおよそ同様です）。

Remark. interval<T> には、下端と上端の 2 引数をとるコンストラクタ以外にも、ベースの型 T の値 x をとって一点区間 [x, x] を作る 1 引数のコンストラクタがあります。

test-interval-arith-IRxR.cpp

# 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.0, 2.0);
	const R b = 3.0;
	cout << "a = " << a << endl;
	cout << "b = " << b << endl;
	cout << "a + b = " << a + b << endl;
	cout << "a - b = " << a - b << endl;
	cout << "a * b = " << a * b << endl;
	cout << "a / b = " << a / b << endl;
	return 0;
}

数同士の四則演算

では、両方の引数をベースの型で与えて四則を行うにはどうすればよいでしょうか。え？？そのまま float や double の built-in operator を使えばいいじゃない？？ダメダメ！ダメですよ！！それでは浮動小数点誤差が生じてしまいます。こんな時には boost::numeric::interval_lib に定義されている関数 add, sub, mul, div を使います。これらは、ベースの型 T の値を 2 つ引数にとって、四則演算の切り捨て・切り上げをそれぞれ下端・上端とするような interval<T> 型の区間を計算して返します。

test-interval-arith-RxR.cpp

# 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 R a = 1.7;
	const R b = 3.3;
	cout << "a = " << a << endl;
	cout << "b = " << b << endl;
	{
		using namespace boost::numeric::interval_lib;
	 	cout << "a + b = " << add<IR>(a, b) << "w(" << width(add<IR>(a, b)) << ")" << endl;
		cout << "a - b = " << sub<IR>(a, b) << "w(" << width(sub<IR>(a, b)) << ")" << endl;
		cout << "a * b = " << mul<IR>(a, b) << "w(" << width(mul<IR>(a, b)) << ")" << endl;
		cout << "a / b = " << div<IR>(a, b) << "w(" << width(div<IR>(a, b)) << ")" << endl;
	}
	return 0;
}

上のプログラムで正しく区間演算が行われれば、次のような結果が得られます。ここで、 width は区間の幅、つまり上端と下端の差を計算する関数です。浮動小数点演算による浮動小数点誤差が非零の区間幅として表れていることが見てとれます。

output-for-test-interval-arith-RxR.txt

a = 1.7
b = 3.3
a + b = [5,5]w(8.88178e-16)
a - b = [-1.6,-1.6]w(0)
a * b = [5.61,5.61]w(8.88178e-16)
a / b = [0.515152,0.515152]w(1.11022e-16)

Remark. もちろん、1 引数コンストラクタを用いて IR(1) / IR(9) などと書いても同様のことを実現できますが、このような一点区間同士の演算だと同じ計算が余分に行われてしまう分 div<IR>(1, 9) の方が若干効率がよいです。

下端・上端・中点・ノルム

区間の幅を計算する関数 width のように、interval<T> 型で与えられた区間に対してなんらかのベースの型 T の値を返す関数として lower, upper, median, norm があります。最初 3 つは文字通りの意味で、それぞれ区間の下端・上端・中点を返す関数です。norm は、区間が含む数の最大絶対値を返します。

訳註(?). 本家ドキュメントには「これは、実数に対する絶対値のように、数学でいう所のノルムになっている」という風なことが書かれていました。しかし先に注意したように、そもそも区間全体の為す空間がベクトル空間ではないので、区間全体の集合が norm によってノルム空間を為す訳ではありません。まぁそうは言ってもあながち間違いでもなくて、実際、ノルムの満たすべき性質（例えば三角不等式）はおよそ満たしています。私はこれをノルムと呼んでいるのを他で見たことがなくて、普通は magnitude などと呼ばれるものだと思います。ちなみに最小絶対値の方は mignitude と呼びますが、Boost の区間演算ライブラリでは特に実装されていなさそうです。

今日もここまで。あまり進まなかった… 次は区間特有の関数に触れたいと思います。ではでは。

今回のサンプルプログラムと出力例: https://gist.github.com/2941941
次: http://qiita.com/items/abcf31a7b1787604a4f4