区間演算をstep by stepで実装する

1. 区間演算をC言語で実装する。

上端と下端をdouble型で持つ区間演算を、C言語で、最も簡単な形で実装してみる。 C言語でよくあるライブラリの作法に従って、実装をinterval.c、そのプロトタイプ宣言をinterval.hに書いている。まずはinterval.cを示す。

interval_add, interval_sub, interval_mul, interval_divが加減乗除を行う関数、interval_sqrtが平方根、 interval_minusが符号反転、interval_init, interval_init2は初期化、 interval_printが表示。

乗除算はちゃんと場合分けせず全種類計算して最大・最小を取る手抜き実装になっている。また、除算と平方根でエラーが起きたときはその場でexitするようになっている。

interval.c

#include "interval.h"

interval interval_init(double x)
{
    interval z;
    z.inf = x;
    z.sup = x;
    return z;
}

interval interval_init2(double x, double y)
{
    interval z;
    z.inf = x;
    z.sup = y;
    return z;
}

interval interval_add(interval x, interval y)
{
    interval z;

    fesetround(FE_DOWNWARD);
    z.inf = x.inf + y.inf;
    fesetround(FE_UPWARD);
    z.sup = x.sup + y.sup;
    fesetround(FE_TONEAREST);

    return z;
}

interval interval_sub(interval x, interval y)
{
    interval z;

    fesetround(FE_DOWNWARD);
    z.inf = x.inf - y.sup;
    fesetround(FE_UPWARD);
    z.sup = x.sup - y.inf;
    fesetround(FE_TONEAREST);

    return z;
}

interval interval_minus(interval x)
{
    interval z;

    z.inf = - x.sup;
    z.sup = - x.inf;

    return z;
}

interval interval_mul(interval x, interval y)
{
    interval z;
    double tmp;

    fesetround(FE_DOWNWARD);
    z.inf = x.inf * y.inf;
    tmp = x.inf * y.sup;
    if (tmp < z.inf) z.inf = tmp;
    tmp = x.sup * y.inf;
    if (tmp < z.inf) z.inf = tmp;
    tmp = x.sup * y.sup;
    if (tmp < z.inf) z.inf = tmp;

    fesetround(FE_UPWARD);
    z.sup = x.inf * y.inf;
    tmp = x.inf * y.sup;
    if (tmp > z.sup) z.sup = tmp;
    tmp = x.sup * y.inf;
    if (tmp > z.sup) z.sup = tmp;
    tmp = x.sup * y.sup;
    if (tmp > z.sup) z.sup = tmp;

    fesetround(FE_TONEAREST);

    return z;
}

interval interval_div(interval x, interval y)
{
    interval z;
    double tmp;

    if (y.inf <= 0 && y.sup >= 0) {
        printf("zero divide!\n");
        exit(1);
    }

    fesetround(FE_DOWNWARD);
    z.inf = x.inf / y.inf;
    tmp = x.inf / y.sup;
    if (tmp < z.inf) z.inf = tmp;
    tmp = x.sup / y.inf;
    if (tmp < z.inf) z.inf = tmp;
    tmp = x.sup / y.sup;
    if (tmp < z.inf) z.inf = tmp;

    fesetround(FE_UPWARD);
    z.sup = x.inf / y.inf;
    tmp = x.inf / y.sup;
    if (tmp > z.sup) z.sup = tmp;
    tmp = x.sup / y.inf;
    if (tmp > z.sup) z.sup = tmp;
    tmp = x.sup / y.sup;
    if (tmp > z.sup) z.sup = tmp;

    fesetround(FE_TONEAREST);

    return z;
}

interval interval_sqrt(interval x)
{
    interval z;

    if (x.inf < 0) {
        printf("sqrt of negative value!\n");
        exit(1);
    }

    fesetround(FE_DOWNWARD);
    z.inf = sqrt(x.inf);

    fesetround(FE_UPWARD);
    z.sup = sqrt(x.sup);

    fesetround(FE_TONEAREST);

    return z;
}

void interval_print(interval x)
{
    printf("[%.17g,%.17g]\n", x.inf, x.sup);
}

以下のinterval.hは、肝心の構造体intervalの宣言と、上のinterval.cで定義された関数のプロトタイプ宣言を集めたもの。構造体の宣言ではtypedefを使って使用時にstructと書かなくていいようにしている。 #ifdef -- #define -- #endifはよくあるインクルードガードで、二度以上includeした場合に最初の一回しか読み込まれないための工夫。

interval.h

#ifndef INTERVAL_H
#define INTERVAL_H

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <fenv.h>

typedef struct {
    double inf;
    double sup;
} interval;

interval interval_init(double x);
interval interval_init2(double x, double y);
interval interval_add(interval x, interval y);
interval interval_sub(interval x, interval y);
interval interval_minus(interval x);
interval interval_mul(interval x, interval y);
interval interval_div(interval x, interval y);
interval interval_sqrt(interval x);
void interval_print(interval x);

#endif /* INTERVAL_H */

以下のtest-interval.cは、これらの関数の使用例。

test-interval.c

#include "interval.h"

int main()
{
    interval x, y, z;

    x = interval_init(1);
    y = interval_init(10);
    z = interval_div(x, y);

    interval_print(z);

    x = interval_init2(1, 2);
    y = interval_init2(3, 4);

    /* basic four operations */
    z = interval_add(x, y);
    interval_print(z);
    z = interval_sub(x, y);
    interval_print(z);
    z = interval_mul(x, y);
    interval_print(z);
    z = interval_div(x, y);
    interval_print(z);

    /* square root */
    z = interval_sqrt(x);
    interval_print(z);

    /* negate */
    z = interval_minus(x);
    interval_print(z);
}

このように、区間演算を使用するプログラムではinterval.hをincludeし、 compile時に本体interval.cと一緒に使う。以下は実行結果。

$ cc test-interval.c interval.c -lm
$ ./a.out 
[0.099999999999999992,0.10000000000000001]
[4,6]
[-3,-1]
[3,8]
[0.25,0.66666666666666674]
[1,1.4142135623730951]
[-2,-1]
$

2. 2次方程式で大きな誤差が出る例題を区間演算してみる

以下は、2次方程式 \[ ax^2 + bx + c = 0 \] で、\(a=1, b=10^{15}, c=10^{14}\)としたときの解のうち大きい方を、そのまま解の公式 \[ x = \frac{-b+\sqrt{b^2 - 4ac}}{2a} \] で計算した場合と、分子を有理化した \[ x = \frac{2c}{-b-\sqrt{b^2 - 4ac}} \] で大きく計算結果が異なってしまう例:

quad-double.c

#include <stdio.h>
#include <math.h>

int main()
{
    double a, b, c, r;

    a = 1.;
    b = 1e15;
    c = 1e14;

    r = (-b + sqrt(b * b - 4 * a * c)) / (2 * a);
    printf("%.17g\n", r);

    r = 2 * c / (-b - sqrt(b * b - 4 * a * c));
    printf("%.17g\n", r);
}

これを実行すると、

$ cc quad-double.c -lm
$ ./a.out
-0.125
-0.10000000000000002
$

のように、数学的には同じ値になるはずのものが大きく異なってしまう。

どちらの値がおかしいかは、区間演算を行うと分かる。次のプログラムは、上のプログラムの計算を区間演算に置き換えたものである。

quad-interval.c

#include "interval.h"

int main()
{
    interval a, b, c, d, r, t1, t2;

    a = interval_init(1);
    b = interval_init(1e15);
    c = interval_init(1e14);

    /* r = (-b + sqrt(b * b - 4 * a * c)) / (2. * a); */

    t1 = interval_mul(b, b);
    t2 = interval_mul(interval_mul(interval_init(4), a), c);
    d = interval_sub(t1, t2);

    d = interval_sqrt(d);

    t1 = interval_sub(d, b);
    t2 = interval_mul(interval_init(2), a);

    r = interval_div(t1, t2);

    interval_print(r);

    /* r = 2 * c / (-b - sqrt(b * b - 4 * a * c)); */

    t1 = interval_mul(b, b);
    t2 = interval_mul(interval_mul(interval_init(4), a), c);
    d = interval_sub(t1, t2);
    d = interval_sqrt(d);

    t1 = interval_mul(interval_init(2), c);
    t2 = interval_minus(interval_add(b, d));

    r = interval_div(t1, t2);

    interval_print(r);
}

これを実行すると、

$ cc quad-interval.c interval.c -lm
$ ./a.out 
[-0.1875,-0.0625]
[-0.10000000000000003,-0.099999999999999992]
$

となる。これを見ると、前者は極端に区間幅が広くなっていて計算結果が信頼できず、後者は区間幅が狭く信頼できる計算であったことが分かる。

ただ、このプログラムを見ると、C言語では自分で定義した構造体intervalに対して加減乗除の演算子が使えないためそれらを関数呼び出しに置き換えており、プログラムを大きく書き換えることを強いられてしまっている。 C言語の拡張であるC++には演算子多重定義の機能があり、自分で定義した型に対して演算子を用いたときの動作を定義することができる。

3. 区間演算のライブラリをC++に書き換える

上で作ったinterval.c/interval.hによるC言語版の区間演算ライブラリをC++で書き換えて、演算子多重定義が使えるようにしてみる。書き換えたポイントは、次の通り。

interval.hpp単一で機能する、header onlyのスタイルにした。
オブジェクト指向言語なので、構造体intervalの定義の中に、その構造体に関係する関数群を一緒に書くスタイル。
interval_init, interval_init2は、コンストラクタと呼ばれる、 interval型が生成されるときに自動的に呼ばれる関数として実装する。コンストラクタは、クラス名(構造体名)と同じ名前で戻り値がない、特殊な関数の形で書く。
(演算子多重定義) operator+という特別な名前の関数を定義すると、 +演算子を使用している箇所がその関数の呼び出しに置き換えられる。
関数の前についているfriendという修飾子は、メンバ関数ではない一般の関数ではあるが、interval型への自由なアクセスが許される、という意味。
平方根は、C言語では同名関数の定義は許されないため、 (C言語が標準で持っているsqrtという名前を避けて) interval_sqrtという名前にしていたが、C++では引数の型が違えば同名の関数の定義が許されているため、同じsqrtという名前が使える。
C++での「std::cout << x;」のような画面出力は、演算子「<<」への多重定義を用いて、C言語のように%d, %fなどを使い分ける必要なしにさまざまな型の表示を実現している。つまり、同様に「<<」演算子の定義をintervalに対して行ってあげれば、 std::coutにinterval型を流し込めるようになる。
除算と平方根での引数のエラーは、C言語では単にプログラムを終了させてしまっていたが、C++版ではC++にしかないエラーハンドリングの機能を使って、throwで例外を発生させている。throwで発生させた例外は、try-catch構文を使うと、単にプログラムが終了する以外の他のことをさせることもできるようになる。

以上のように書き換えたものを示す。

interval.hpp

#ifndef INTERVAL_HPP
#define INTERVAL_HPP

#include <iostream>
#include <cmath>
#include <stdexcept>
#include <fenv.h>

struct interval {

    double inf;
    double sup;

    interval() {
        inf = 0;
        sup = 0;
    }

    interval(double x) {
        inf = x;
        sup = x;
    }

    interval(double x, double y) {
        inf = x;
        sup = y;
    }

    friend interval operator+(interval x, 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-(interval x, 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-(interval x) {
        interval r;

        r.sup = - x.inf;
        r.inf = - x.sup;

        return r;
    }

    friend interval operator*(interval x, interval y) {
        interval r;
        double tmp;

        fesetround(FE_DOWNWARD);
        r.inf = x.inf * y.inf;
        tmp = x.inf * y.sup;
        if (tmp < r.inf) r.inf = tmp;
        tmp = x.sup * y.inf;
        if (tmp < r.inf) r.inf = tmp;
        tmp = x.sup * y.sup;
        if (tmp < r.inf) r.inf = tmp;

        fesetround(FE_UPWARD);
        r.sup = x.inf * y.inf;
        tmp = x.inf * y.sup;
        if (tmp > r.sup) r.sup = tmp;
        tmp = x.sup * y.inf;
        if (tmp > r.sup) r.sup = tmp;
        tmp = x.sup * y.sup;
        if (tmp > r.sup) r.sup = tmp;

        fesetround(FE_TONEAREST);

        return r;
    }

    friend interval operator/(interval x, interval y) {
        interval r;
        double tmp;

        if (y.inf <= 0. && y.sup >= 0.) {
            throw std::domain_error("interval: division by 0");
        }

        fesetround(FE_DOWNWARD);
        r.inf = x.inf / y.inf;
        tmp = x.inf / y.sup;
        if (tmp < r.inf) r.inf = tmp;
        tmp = x.sup / y.inf;
        if (tmp < r.inf) r.inf = tmp;
        tmp = x.sup / y.sup;
        if (tmp < r.inf) r.inf = tmp;

        fesetround(FE_UPWARD);
        r.sup = x.inf / y.inf;
        tmp = x.inf / y.sup;
        if (tmp > r.sup) r.sup = tmp;
        tmp = x.sup / y.inf;
        if (tmp > r.sup) r.sup = tmp;
        tmp = x.sup / y.sup;
        if (tmp > r.sup) r.sup = tmp;

        fesetround(FE_TONEAREST);

        return r;
    }

    friend interval sqrt(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;
    }

    friend std::ostream& operator<<(std::ostream& s, interval x) {
        s << '[' << x.inf << ',' << x.sup << ']';
        return s;
    }

};

#endif // INTERVAL_HPP

使用例は以下の通り。最後の行の「x+1」がうまくコンパイル出来ているのが少し不思議ではあるが、自動的に引数が1つのコンストラクタを呼び出して、「operator+(x, interval(1))」のように処理されている。

test-interval.cc

#include "interval.hpp"

int main()
{
    std::cout.precision(17);

    // constructors
    interval x;
    interval y(1);
    interval z(2, 3);

    x = 1;
    y = 10;
    std::cout << x / y << "\n";

    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";

    // square root
    std::cout << sqrt(x) << "\n";

    // negate
    std::cout << -x << "\n";

    // operation with constant
    std::cout << x + 1 << "\n";
}

以下は実行結果。

$ c++ test-interval.cc 
$ ./a.out 
[0.099999999999999992,0.10000000000000001]
[4,6]
[-3,-1]
[3,8]
[0.25,0.66666666666666674]
[1,1.4142135623730951]
[-2,-1]
[2,3]
$

このC++のライブラリ使うと、2次方程式の例題の区間演算化は、ほぼ変数宣言の型doubleをintervalに書き換えるだけで済むようになる。

quad-interval.cc

#include <iostream>
#include <cmath>
#include "interval.hpp"

int main()
{
    interval a, b, c, r;

    std::cout.precision(17);

    a = 1.;
    b = 1e15;
    c = 1e14;

    r = (-b + sqrt(b * b - 4 * a * c)) / (2 * a);
    std::cout << r << "\n";

    r = 2 * c / (-b - sqrt(b * b - 4 * a * c));
    std::cout << r << "\n";
}

この実行結果は以下の通り。

$ c++ quad-interval.cc 
$ ./a.out 
[-0.1875,-0.0625]
[-0.10000000000000003,-0.099999999999999992]
$

4. C++のコードをもうちょっとまともに

上で書いたC++のコードをもうちょっと改善し、もう少しC++らしく直してみる。

interval xに対する「x+1」のようなコードに対してコンストラクタによる自動型変換に任せず、「operator+(interval, x, double y)」のような片方がdoubleである場合専用の関数を別に作る。
乗除算を全種類計算するような雑なことをせず引数の正負で場合分けを行う。
+=, -=, *=, /=のような演算子を使えるようにする。
inf, supはprivateにして、lower(), upper()というアクセサ(accessor)経由でアクセスするようにする。
関数に対するdouble型の引数の受け渡しを「const参照」にして、データのコピーを減らす。

以上の方針で書き換えたものを示す。

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 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;
    }

    friend std::ostream& operator<<(std::ostream& s, const interval& x) {
        s << '[' << x.inf << ',' << x.sup << ']';
        return s;
    }

    const double& lower() const {
        return inf;
    }
    const double& upper() const {
        return sup;
    }
    double& lower() {
        return inf;
    }
    double& upper() {
        return sup;
    }
};

#endif // INTERVAL_HPP

使用例は以下の通り。

test-interval.cc

#include "interval.hpp"

int main()
{
    std::cout.precision(17);

    // constructors
    interval x;
    interval y(1);
    interval z(2, 3);

    x = 1;
    y = 10;
    std::cout << x / y << "\n";

    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";

    // square root
    std::cout << sqrt(x) << "\n";

    // negate
    std::cout << -x << "\n";

    // operation with constant
    std::cout << x + 1 << "\n";

    // compound assignment operator
    z += x;
    z += 1;
    std::cout << z << "\n";

    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";
}

実行結果は以下の通り。

$ c++ test-interval.cc 
$ ./a.out 
[0.099999999999999992,0.10000000000000001]
[4,6]
[-3,-1]
[3,8]
[0.25,0.66666666666666674]
[1,1.4142135623730951]
[-2,-1]
[2,3]
[4,6]
3
4
[3.5,4]
$

Implementing Interval Arithmetic step by step