## 2015/11/27(金)区間演算の実装について(1)

```#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
```

```#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.);
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つの解のうちの大きい方を、
• 普通に解の公式を使う
• 解の公式の分子を有理化する
の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]
```
となり、前者の区間幅が極端に広いことで、計算結果が怪しいと気づくことが出来ます。

さて、学習用としてはこれでいいのですが、実際に使おうとするとこのレベルの実装ではいろいろと問題点があります。少し記事が長くなったので、問題点は次の記事で。
OK キャンセル 確認 その他