2015/11/27(金)区間演算の実装について(1)
区間演算は原理は簡単ですが、実際に実装となると案外ハマリポイントが多く、難しいものです。簡単な区間演算の実装例を通して、少し解説してみようと思います。
使いやすい区間演算を実装するには、区間型を定義してその区間型を組み込み型と同等に扱える、いわゆる演算子多重定義(operator overloading)のできるプログラミング言語が望ましいです。自分はC++を愛用していますので、ここでもC++で実装します。両端点の型はdouble、方向付き丸めはC99のfenv.hとfesetroundを使います。
最も単純で最低限の機能のみの実装を示します (interval.hpp)。
#ifndef INTERVAL_HPP
#define INTERVAL_HPP
#include <iostream>
#include <cmath>
#include <stdexcept>
#include <fenv.h>
class interval {
double inf;
double sup;
public:
interval() {
inf = 0.;
sup = 0.;
}
interval(const double& x) {
inf = x;
sup = x;
}
interval(const double& x, const double& y) {
inf = x;
sup = y;
}
friend interval operator+(const interval& x, const interval& y) {
interval r;
fesetround(FE_DOWNWARD);
r.inf = x.inf + y.inf;
fesetround(FE_UPWARD);
r.sup = x.sup + y.sup;
fesetround(FE_TONEAREST);
return r;
}
friend interval operator+(const interval& x, const double& y) {
interval r;
fesetround(FE_DOWNWARD);
r.inf = x.inf + y;
fesetround(FE_UPWARD);
r.sup = x.sup + y;
fesetround(FE_TONEAREST);
return r;
}
friend interval operator+(const double& x, const interval& y) {
interval r;
fesetround(FE_DOWNWARD);
r.inf = x + y.inf;
fesetround(FE_UPWARD);
r.sup = x + y.sup;
fesetround(FE_TONEAREST);
return r;
}
friend interval& operator+=(interval& x, const interval& y) {
x = x + y;
return x;
}
friend interval& operator+=(interval& x, const double& y) {
fesetround(FE_DOWNWARD);
x.inf = x.inf + y;
fesetround(FE_UPWARD);
x.sup = x.sup + y;
fesetround(FE_TONEAREST);
return x;
}
friend interval operator-(const interval& x, const interval& y) {
interval r;
fesetround(FE_DOWNWARD);
r.inf = x.inf - y.sup;
fesetround(FE_UPWARD);
r.sup = x.sup - y.inf;
fesetround(FE_TONEAREST);
return r;
}
friend interval operator-(const interval& x, const double& y) {
interval r;
fesetround(FE_DOWNWARD);
r.inf = x.inf - y;
fesetround(FE_UPWARD);
r.sup = x.sup - y;
fesetround(FE_TONEAREST);
return r;
}
friend interval operator-(const double& x, const interval& y) {
interval r;
fesetround(FE_DOWNWARD);
r.inf = x - y.sup;
fesetround(FE_UPWARD);
r.sup = x - y.inf;
fesetround(FE_TONEAREST);
return r;
}
friend interval& operator-=(interval& x, const interval& y) {
x = x - y;
return x;
}
friend interval& operator-=(interval& x, const double& y) {
fesetround(FE_DOWNWARD);
x.inf = x.inf - y;
fesetround(FE_UPWARD);
x.sup = x.sup - y;
fesetround(FE_TONEAREST);
return x;
}
friend interval operator-(const interval& x) {
interval r;
r.sup = - x.inf;
r.inf = - x.sup;
return r;
}
friend interval operator*(const interval& x, const interval& y) {
interval r;
double tmp;
if (x.inf >= 0.) {
if (y.inf >= 0.) {
fesetround(FE_DOWNWARD);
r.inf = x.inf * y.inf;
fesetround(FE_UPWARD);
r.sup = x.sup * y.sup;
} else if (y.sup <= 0.) {
fesetround(FE_DOWNWARD);
r.inf = x.sup * y.inf;
fesetround(FE_UPWARD);
r.sup = x.inf * y.sup;
} else {
fesetround(FE_DOWNWARD);
r.inf = x.sup * y.inf;
fesetround(FE_UPWARD);
r.sup = x.sup * y.sup;
}
} else if (x.sup <= 0.) {
if (y.inf >= 0.) {
fesetround(FE_DOWNWARD);
r.inf = x.inf * y.sup;
fesetround(FE_UPWARD);
r.sup = x.sup * y.inf;
} else if (y.sup <= 0.) {
fesetround(FE_DOWNWARD);
r.inf = x.sup * y.sup;
fesetround(FE_UPWARD);
r.sup = x.inf * y.inf;
} else {
fesetround(FE_DOWNWARD);
r.inf = x.inf * y.sup;
fesetround(FE_UPWARD);
r.sup = x.inf * y.inf;
}
} else {
if (y.inf >= 0.) {
fesetround(FE_DOWNWARD);
r.inf = x.inf * y.sup;
fesetround(FE_UPWARD);
r.sup = x.sup * y.sup;
} else if (y.sup <= 0.) {
fesetround(FE_DOWNWARD);
r.inf = x.sup * y.inf;
fesetround(FE_UPWARD);
r.sup = x.inf * y.inf;
} else {
fesetround(FE_DOWNWARD);
r.inf = x.inf * y.sup;
tmp = x.sup * y.inf;
if (tmp < r.inf) r.inf = tmp;
fesetround(FE_UPWARD);
r.sup = x.inf * y.inf;
tmp = x.sup * y.sup;
if (tmp > r.sup) r.sup = tmp;
}
}
fesetround(FE_TONEAREST);
return r;
}
friend interval operator*(const interval& x, const double& y) {
interval r;
if (y >= 0.) {
fesetround(FE_DOWNWARD);
r.inf = x.inf * y;
fesetround(FE_UPWARD);
r.sup = x.sup * y;
} else {
fesetround(FE_DOWNWARD);
r.inf = x.sup * y;
fesetround(FE_UPWARD);
r.sup = x.inf * y;
}
fesetround(FE_TONEAREST);
return r;
}
friend interval operator*(const double& x, const interval& y) {
interval r;
if (x >= 0.) {
fesetround(FE_DOWNWARD);
r.inf = x * y.inf;
fesetround(FE_UPWARD);
r.sup = x * y.sup;
} else {
fesetround(FE_DOWNWARD);
r.inf = x * y.sup;
fesetround(FE_UPWARD);
r.sup = x * y.inf;
}
fesetround(FE_TONEAREST);
return r;
}
friend interval& operator*=(interval& x, const interval& y) {
x = x * y;
return x;
}
friend interval& operator*=(interval& x, const double& y) {
x = x * y;
return x;
}
friend interval operator/(const interval& x, const interval& y) {
interval r;
if (y.inf > 0.) {
if (x.inf >= 0.) {
fesetround(FE_DOWNWARD);
r.inf = x.inf / y.sup;
fesetround(FE_UPWARD);
r.sup = x.sup / y.inf;
} else if (x.sup <= 0.) {
fesetround(FE_DOWNWARD);
r.inf = x.inf / y.inf;
fesetround(FE_UPWARD);
r.sup = x.sup / y.sup;
} else {
fesetround(FE_DOWNWARD);
r.inf = x.inf / y.inf;
fesetround(FE_UPWARD);
r.sup = x.sup / y.inf;
}
} else if (y.sup < 0.) {
if (x.inf >= 0.) {
fesetround(FE_DOWNWARD);
r.inf = x.sup / y.sup;
fesetround(FE_UPWARD);
r.sup = x.inf / y.inf;
} else if (x.sup <= 0.) {
fesetround(FE_DOWNWARD);
r.inf = x.sup / y.inf;
fesetround(FE_UPWARD);
r.sup = x.inf / y.sup;
} else {
fesetround(FE_DOWNWARD);
r.inf = x.sup / y.sup;
fesetround(FE_UPWARD);
r.sup = x.inf / y.sup;
}
} else {
fesetround(FE_TONEAREST);
throw std::domain_error("interval: division by 0");
}
fesetround(FE_TONEAREST);
return r;
}
friend interval operator/(const interval& x, const double& y) {
interval r;
if (y > 0.) {
fesetround(FE_DOWNWARD);
r.inf = x.inf / y;
fesetround(FE_UPWARD);
r.sup = x.sup / y;
} else if (y < 0.) {
fesetround(FE_DOWNWARD);
r.inf = x.sup / y;
fesetround(FE_UPWARD);
r.sup = x.inf / y;
} else {
fesetround(FE_TONEAREST);
throw std::domain_error("interval: division by 0");
}
fesetround(FE_TONEAREST);
return r;
}
friend interval operator/(const double& x, const interval& y) {
interval r;
if (y.inf > 0. || y.sup < 0.) {
if (x >= 0.) {
fesetround(FE_DOWNWARD);
r.inf = x / y.sup;
fesetround(FE_UPWARD);
r.sup = x / y.inf;
} else {
fesetround(FE_DOWNWARD);
r.inf = x / y.inf;
fesetround(FE_UPWARD);
r.sup = x / y.sup;
}
} else {
fesetround(FE_TONEAREST);
throw std::domain_error("interval: division by 0");
}
fesetround(FE_TONEAREST);
return r;
}
friend interval& operator/=(interval& x, const interval& y) {
x = x / y;
return x;
}
friend interval& operator/=(interval& x, const double& y) {
x = x / y;
return x;
}
friend std::ostream& operator<<(std::ostream& s, const interval& x) {
s << '[';
s << x.inf;
s << ',';
s << x.sup;
s << ']';
return s;
}
friend interval sqrt(const interval& x) {
interval r;
if (x.inf < 0.) {
throw std::domain_error("interval: sqrt of negative value");
}
fesetround(FE_DOWNWARD);
r.inf = sqrt(x.inf);
fesetround(FE_UPWARD);
r.sup = sqrt(x.sup);
fesetround(FE_TONEAREST);
return r;
}
const double& lower() const {
return inf;
}
const double& upper() const {
return sup;
}
double& lower() {
return inf;
}
double& upper() {
return sup;
}
};
#endif // INTERVAL_HPP
使える演算は加減乗除と平方根で、coutでの表示にも対応してます。lowerとupperは下限と上限へのアクセサ。乗除算は手抜きせずに一応ちゃんと符号による場合分けを行っています。C++の初心者でも入門書を読みながら実装できるレベルでしょう。次にこのライブラリの使用例を示します (test-interval.cc)。
#include "interval.hpp"
int main()
{
interval x;
interval y = 1.;
interval z(1.);
x = 1.;
y = 10.;
z = x / y;
std::cout << z << "\n";
std::cout.precision(17);
std::cout << z << "\n";
// copy
x = interval(1., 2.);
y = interval(3., 4.);
// basic four operations
std::cout << x + y << "\n";
std::cout << x - y << "\n";
std::cout << x * y << "\n";
std::cout << x / y << "\n";
// operation with constant
std::cout << x + 1 << "\n";
std::cout << x + 1. << "\n";
// compound assignment operator
z += x;
z += 1.;
z = interval(3., 4.);
// access to endpoints
std::cout << z.lower() << "\n";
std::cout << z.upper() << "\n";
z.lower() = 3.5;
std::cout << z << "\n";
}
実行すると、[0.1,0.1] [0.099999999999999992,0.10000000000000001] [4,6] [-3,-1] [3,8] [0.25,0.66666666666666674] [2,3] [2,3] 3 4 [3.5,4]のようになります。この2つのファイルを一応アップロードしておきます。
interval-supersimple.tar.gz
これを使って、大きく誤差が出る計算をしてみます。2次方程式ax2+bx+c=0でa=1, b=1015, c=1014として、2つの解のうちの大きい方を、
- 普通に解の公式を使う
- 解の公式の分子を有理化する
まず普通にdoubleで。
#include <iostream>
#include <cmath>
int main()
{
double a, b, c, x, y;
std::cout.precision(17);
a = 1.;
b = 1e15;
c = 1e14;
x = (-b + sqrt(b * b - 4. * a * c)) / (2. * a);
y = 2 * c / (-b - sqrt(b * b - 4. * a * c));
std::cout << x << "\n";
std::cout << y << "\n";
}
のようなプログラムを実行すると、-0.125 -0.10000000000000002のようになります。この計算は後者がほぼ正しく、前者は大きな誤差が入っています。
これを区間演算で計算してみます。プログラムは、
#include <iostream>
#include <cmath>
#include "interval.hpp"
int main()
{
interval a, b, c, x, y;
std::cout.precision(17);
a = 1.;
b = 1e15;
c = 1e14;
x = (-b + sqrt(b * b - 4. * a * c)) / (2. * a);
y = 2 * c / (-b - sqrt(b * b - 4. * a * c));
std::cout << x << "\n";
std::cout << y << "\n";
}
のような感じ。演算子多重定義のおかげで最小限の変更で区間演算が行えます。計算結果は、[-0.1875,-0.0625] [-0.10000000000000003,-0.099999999999999992]となり、前者の区間幅が極端に広いことで、計算結果が怪しいと気づくことが出来ます。
さて、学習用としてはこれでいいのですが、実際に使おうとするとこのレベルの実装ではいろいろと問題点があります。少し記事が長くなったので、問題点は次の記事で。