最終更新: 2022/1/6

区間演算ライブラリ

柏木 雅英

1. はじめに

2. 特徴

3. サポートする環境

4. ファイル構成

interval.hpp
(下請け: convert.hpp, constants.hpp)
rdouble.hpp
(下請け: conv-double.hpp, hwround.hpp, rdouble-hwround.hpp, rdouble-nohwround.hpp)

なお、doubleに対してinterval.hppだけincludeした場合は、区間演算らしき計算をするが 両端点の値を計算するときに丸めの変更は行われず、 全く精度保証されないことに注意。 有理数などの丸めの変更の必要の無い型に対してはこれでもいいが、doubleを 使う場合はrdouble.hppも必ず一緒にincludeすること。

5. 使い方

5.1 基本的な使い方 

#include <kv/interval.hpp>
のようにincludeすれば使える。ただし、前述したように、 丸めの向きの変更の必要のある型を内部型として使う場合は、 そのためのheader fileも必要。例えば内部がdoubleの場合は、 
#include <kv/interval.hpp>
#include <kv/rdouble.hpp>
のようにrdouble.hppも同時にincludeすること。

interval型はkv名前空間の中にある。またtemplate型であり、 例えば内部にdoubleを入れるならば、interval<double> という型名で 使用する。従って、例えば変数宣言は 

kv::interval<double> x;
のようになる。kv::を付けるのが面倒だと思うなら、using namespace kv; などとすれば省略出来る。または、 
typedef kv::interval<double> itv;
などとtypedefして、 
itv x;
のように使うのもお勧め。 以下簡単のためこの表記で説明する。

5.2 宣言と代入

変数の宣言は、 
    itv x;
のようにする。宣言と同時に初期化するには、 
    itv x(1.);
    itv x(1., 2.);
のような感じ。引数一つのコンストラクタがexplicitなので、 
    itv x = 1.;
と書くとエラーになることに注意。 
    itv x = (itv)1.;
    itv x = itv(1.);
    itv x = itv(1., 2.);
は大丈夫。

代入演算子はオーバーロードしてあって、 

    itv x;
    x = 1.;
    x = (itv)1.;
    x = itv(1.);
    x = itv(1., 2.);
はいずれもエラーにならない(一つ目の書き方が有効)。
注意   初期化や代入時の引数 (上の1.や2.) は、intervalの 内部に使っている型だけでなく、「内部に使っている型に変換出来る型」も 指定可能。例えば interval<dd>型のxに対してx = 1.とか出来る。 このあたりの詳細は ここを参照。
また、 
    x.assign(1., 2.);
のようにも書ける(boost.intervalとの互換性のため)。

また、文字列での初期化及び代入が可能で、 

    itv x("0.1");
    itv x("0.1", "0.2");
や、 
    itv x;
    x = "0.1";
    x = (itv)"0.1";
    x = itv("0.1");
    x = itv("0.1", "0.2");
のように書ける。この場合、"0.1"を指定すれは 下端は0.1以下の最大のdouble、 上端は0.1以上の最小のdoubleとなり、真の0.1を包含する区間になる。 ("0.1", "0.2")なら、 下端は0.1以下の最大のdouble、 上端は0.2以上の最小のdouble。

ただし、文字列の変換処理は遅いので、 ループの中など多数回実行される箇所で使用する場合は、予めループの外で 変数に代入するようにするなど工夫すること。

5.3 表示 

    std::cout << z << "\n";
のように書ける。ここでの2進10進変換も、正確に表現出来ない場合は 外側に丸められる。桁数はcout.precisionの設定に従うので、 
    std::cout.precision(17);
    std::cout << z << "\n";
とすれば17桁で表示される。

5.4 四則演算 

    itv x, y, z;
    x = 1.; y = 2.;

    z = x + y;
    z = x - y;
    z = - x;
    z = x * y;
    z = x / y;

    z = x + 1.;
    z = 2. + x;
    z = x + "0.1";
    
    x += y;
    x += 1.;
    x += "0.1";
など。定数との演算で使える定数(上の1.や2.)は、内部で使っている型だけでなく、 「内部で使っている型に変換出来る型」も指定可能。

5.5 数学関数 

    std::cout << sqrt(itv(2.5, 3.5)) << "\n";
    std::cout << exp(itv(2.5, 3.5)) << "\n";
    std::cout << expm1(itv(-0.25, 0.25)) << "\n";
    std::cout << log(itv(0.75, 1.25)) << "\n";
    std::cout << log(itv(0., 1.)) << "\n";
    std::cout << log1p(itv(-0.25, 0.25)) << "\n";
    std::cout << sin(itv(-0.25, 0.25)) << "\n";
    std::cout << cos(itv(-0.25, 0.25)) << "\n";
    std::cout << tan(itv(-0.25, 0.25)) << "\n";
    std::cout << atan(itv(-0.25, 0.25)) << "\n";
    std::cout << asin(itv(-0.25, 0.25)) << "\n";
    std::cout << acos(itv(-0.25, 0.25)) << "\n";
    std::cout << atan2(itv(1.), itv(1.)) << "\n";
    std::cout << sinh(itv(-0.25, 0.25)) << "\n";
    std::cout << cosh(itv(-0.25, 0.25)) << "\n";
    std::cout << tanh(itv(-0.25, 0.25)) << "\n";
    std::cout << asinh(itv(-0.25, 0.25)) << "\n";
    std::cout << acosh(itv(1.5, 2.)) << "\n";
    std::cout << atanh(itv(-0.25, 0.25)) << "\n";
    // integer power
    std::cout << pow(itv(2., 3.), -2) << "\n";
    // general power
    std::cout << pow(itv(2., 3.), itv(2., 3)) << "\n";
など。rdouble.hppをincludeするなど、きちんと内部に使った型に対する 丸めが考慮されていれば、これらの数学関数の計算結果は 全て精度保証される。 が、必ずしも最小の幅の結果を返すわけではなく、浮動小数点数数個分の 区間幅を持つ結果が返ることが多い。また、速度もあまり速くない。

数学関数の実装方法の詳細は、 「シンプルな区間数学関数の実装 (Version:2018/9/5)」 を見て欲しい。

5.6 等号、不等号 

    std::cout << (x < y) << "\n";
    std::cout << (x < 1.) << "\n";
    std::cout << (1. < x) << "\n";
    std::cout << (itv(1.) == 1.) << "\n";
    std::cout << (itv(1., 2.) != 3.) << "\n";
    std::cout << ("0.20001" > x) << "\n";
など。区間内の全ての点に対してその等号または不等号が成立するとき真になる。

5.7 端点へのアクセス

下端、上端の値を取り出したり、書き換えたりすることが出来る。 
    // access to endpoints
    std::cout << x.lower() << "\n";
    std::cout << x.upper() << "\n";
    x.lower() = 3.5;

5.8 division_part1, division_part2

除算 (区間x, yに対する x / y) は、yが0を含んでいたら0除算のエラーとなる 仕様になっている。しかし、端点として±∞を使ってよいならば、 例えば1/[-1,1]は[-∞,-1]∪[1,∞]を返す仕様にすることも考えられる。

division_part1, division_part2は、そのような仕様の除算を行うためのものである。 従って、端点として±∞が返される可能性がある。 区間が2つの集合に分かれる場合があるので、division_part1でbool値partedを 返し、もしこれが偽ならば区間は分かれていない、真ならば区間が分かれているので division_part2を呼ぶともう一つの区間が計算できる、という仕様。 (この仕様はboost.intervalと同じ) 

    // division_part1, division_part2
    // calculate X / (Y ∖ 0). the result may be divided into two part.
    bool parted;
    std::cout << division_part1(itv(1., 2.), itv(-3., 4.), parted) << "\n";
    // if the result has division_part2, parted is set to true.
    if (parted) std::cout << division_part2(itv(1., 2.), itv(-3., 4.)) << "\n";
計算方法の具体的な場合分けは、 端点に無限大を持つ区間演算まとめの除算の節の通り。 除数が真に[0,0]の場合に限って、0除算のエラーとなる。

5.9 その他の機能

その他、以下のような機能がある。 
    // static functions
    std::cout << itv::whole() << "\n";
    std::cout << itv::hull(1., 2.) << "\n";
    std::cout << itv::hull(1., z) << "\n";
    std::cout << itv::hull(x, y) << "\n";
wholeは実数全体(-∞, ∞)を返す。
hullはその引数を含む区間を返す。引数は、両方T、片方がTで片方がinterval<T>、両方interval<T> がある。 
    std::cout << width(z) << "\n";
    std::cout << rad(z) << "\n";
    std::cout << median(z) << "\n";
    std::cout << mid(z) << "\n";
widthは区間幅。「上端-下端」を上向き丸めで計算したもの。
radは区間の半径。「(上端-下端)/2」を上向き丸めで計算したもの。
medianとmidは区間の中点。「(上端+下端)/2」を特に丸め制御無しで計算したもの。戻り値はTで、真の中点とは ずれている可能性がある。 
    std::cout << norm(z) << "\n";
    std::cout << mag(z) << "\n";
    std::cout << abs(itv(2., 3.)) << "\n";
    std::cout << mig(itv(-2., 1.)) << "\n";
absは絶対値。区間Iのabsは、abs(I)={abs(x) | x ∈ I} のように、Iに含まれる点の絶対値の取る集合(区間) を表す。
magとnormは同じで、mag(I)=norm(I)=max{abs(x) | x ∈ I}のように、Iに含まれる点の絶対値の最大値。戻り値はT型。
migは、mig(I)=min{abs(x) | x ∈ I}のように、Iに含まれる点の絶対値の最小値。戻り値はT型。 
    std::cout << in(3.9, z) << "\n";
    std::cout << subset(x, y) << "\n";
    std::cout << proper_subset(y, z) << "\n";
    std::cout << overlap(x, y) << "\n";
    std::cout << intersect(y, z) << "\n";
in(x, I)は、x ∈ Iなら真。xはT型でIはinterval<T>型。
subset(I, J)は、I ⊂ Jなら真。I, Jはinterval<T>型。
proper_subset(I, J)は、 I ⊂ int(J) (Jの内点集合)なら IがJの真部分集合なら真。I, Jはinterval<T>型。
overlap(I, J)は I ∩ J ≠ ∅なら真。
intersect(I, J)は I ∩ J を返す。ただし、overlap(I, J)が真の場合しか使ってはいけない。(本ライブラリは空集合を表せない。) 
    // numeric constants
    std::cout << kv::constants<itv>::pi() << "\n";
    std::cout << kv::constants<itv>::e() << "\n";
    std::cout << kv::constants<itv>::ln2() << "\n";
予め用意された定数として、pi, e, loge2が使える。

5.10 midrad

midとradを同時に計算するmidrad関数がある。 
    midrad(I, m, r);
で、区間Iの中心(の近似)がmに、半径がrに得られる。 この関数は、上端下端型の区間を中心と半径で包含するときに用いるためのものである。 mはmid(I)と同じだが、rはmax(I.upper()-m, m-I.lower())で計算される。 単に 
    m = mid(I);
    r = rad(I);
とすると、mが真の中心からずれる可能性があるため、中心m、半径rの区間はIを包含しない可能性がある。

5.11 max, min

2つの区間の最大値、最小値を区間として得る関数。 
    I3 = max(I1, I2);
    I4 = min(I1, I2);
max([a,b],[c,d])=[max(a,c),max(b,d)]、 min([a,b],[c,d])=[min(a,c),min(b,d)]である。 この関数を使わずに例えばstd::maxが使われてしまうと、例えばmax(x,y)の実装が 
  (x > y) ? x : y;
のような感じだった場合、 
  max([2,4],[1,3]) = [1,3]
と予期せぬ結果が得られてしまう。

6. double以外の型を内部型として使う

まだちゃんと書いてない。多分、 が出来ればいいはず。

精度保証するためには、更にdoubleにとってのrdouble.hppに相当する、 丸めを扱うための機能が必要。これについては次節のカスタマイズを参照。

7. カスタマイズ

区間演算ライブラリのカスタマイズを見てください。

8. 内部型の異なる区間同士の変換

内部型の異なる区間同士の変換を見てください。

9. その他