kv-0.4.38 - kashiの日記

2016/12/10(土)kv-0.4.38

kvライブラリを0.4.38に更新しました。

今回は、区間演算にmidrad、max、minを追加しました。また、dka.hppの中にあるvdka (精度保証付きDKA法) で、解が重複する場合の処理をきちんと行なうようにしました。

midradとmax/minの追加は、どちらも同僚のS先生の要請によるものです。少し詳しく説明します。

midrad

midradは、区間 I に対して、

midrad(I, m, r);

とすると中心の近似mと、半径rを同時に計算します。これは、中心m、半径rの区間で I を包含する、すなわち、上端下端型の区間から中心半径型の区間への変換に使われることを意図しています。単純に

m = mid(I);
r = rad(I);

とすれば良さそうなものですが、中心mに丸め誤差が入るので、これだと I を包含しない可能性があるのです。例えば、ε=2^-52として、I=[1,1+3ε]とします。このとき、midとradを別々に計算すると、

m = 1+2ε
r = 1.5ε

となりますが、[m-r,m+r]は I を含まないことが分かります。midradを使うと、ちゃんと

m = 1+2ε
r = 2ε

になります。ミスが起こりやすいポイントなので気をつけて下さい。

max, min

max, minは区間同士のmax, minです。max(I1, I2)は、x1∈I1、x2∈I2としたときのmax(x1,x2)の取り得る範囲の区間を返します。minも同様です。計算は、

max([a,b],[c,d]) = [max(a,c),max(b,d)]
min([a,b],[c,d]) = min(a,c),max(b.d)]

のように行います。

この関数は、後述するようにVisual C++と少し相性が悪いこともあって実装していなかったのですが、実装しないと思わぬバグに遭遇する可能性があると同僚のS先生に指摘され、実装することにしました。例えば、max, minを実装していないversionのkvで、

#include <kv/interval.hpp>
#include <kv/rdouble.hpp>

typedef kv::interval<double> itv;

using namespace std;

int main()
{
    cout << max(itv(2., 4.), itv(3., 5.)) << endl;
    cout << min(itv(3., 5.), itv(2., 4.)) << endl;
}

のようなプログラムを動かすと、予想に反して

[2,4]
[3,5]

という結果が帰って来ます。これは、std::max, std::minが呼び出されてしまっており、恐らくmaxは

template <class T> 
T max(const T& a, const T& b)
{
    return (a < b) ? b : a;
}

のような実装になっていて、[2,4]<[3,5]は偽なので[2,4]が返されてしまっていると推測されます。よく考えれば分かるとは言え、これは事故を誘発しやすいので、ちゃんと区間用のmaxとminを実装しました。kv-0.4.38で上のプログラムを動かすと、ちゃんと

[3,5]
[2,4]

となります。

なお、Visual C++はwindows.hというヘッダファイルを持っており、これをincludeするとマクロでmaxとminを定義してしまうという大変極悪な仕様になっています。これがincludeされていると、max,minを定義している箇所がマクロ置換されてしまいinterval.hppがコンパイル出来ません。windows.hとkvをどうしても同時に使いたければ、

#define NOMINMAX
#include <windows.h>

のようにwindows.hをincludeする前にNOMINMAXを定義して下さい。