2015/11/27(金)区間演算の実装について(2)
まず、入出力の問題です。区間型の内部で使われているdoubleは2進数で表現されているので、区間型に数値を代入するとき10進数→2進数変換が、区間型の数値を表示するとき2進数→10進数変換が行われます。精度保証付き数値計算を行うためには、この変換が適切に行われる必要があります。例えば、
#include "interval.hpp"
int main()
{
interval x;
std::cout.precision(17);
x = 0.1;
std::cout << x << "\n";
}
のようなプログラムを実行すると、[0.10000000000000001,0.10000000000000001]となってしまい、なぜかxは0.1を含んでいません。これは精度保証付き数値計算のビギナーが陥りやすい罠で、ソーステキストに書いた「0.1」はコンパイラが2進数に変換しますが、0.1は2進数だと無限小数になるのでどこかで打ち切る必要があり、正確に0.1になりません。実際、2進数53桁で0.1に最も近い数は、
1.1001100110011001100110011001100110011001100110011010 x 2-4になります。これは10進数だと、
0.1000000000000000055511151231257827021181583404541015625という数です。頑張って工夫して、例えば
fesetround(FE_DOWNWARD);
x.lower() = 0.1;
fesetround(FE_UPWARD);
x.upper() = 0.1;
と書いたとしても、この0.1を2進数に変換するのはコンパイル時なので、どちらの0.1も上で示した近似値になってしまい、意図通りの0.1を含む区間にはなりません。(あまり大きくない)整数は2進数で正しく表現できるので、
x = 1;
x /= 10;
のように書けば一応0.1を含む区間をxに入れることは出来ますが、毎回このように書くのは大変です。根本的に解決するには、ソーステキストの「0.1」をコンパイラに変換させないこと、すなわちソーステキストのレベルでは文字列で「"0.1"」のように保持し、それの2進数への変換はこちらでやってあげることになります。出力でも同じ問題があって、例えば
#include "interval.hpp"
int main()
{
interval x, y, z;
x = 1;
y = 10;
z = x / y;
std::cout << z << "\n";
std::cout.precision(17);
std::cout << z << "\n";
}
のようなプログラムを実行すると[0.1,0.1] [0.099999999999999992,0.10000000000000001]のようになり、前者はまるで幅の無い区間のように見えてしまいます。これは、表示をcout(というか標準ライブラリ)に任せているせいで、本当は10進文字列への変換を自力で行い、下端は下向き丸めで、上端は上向き丸めで行わなければなりません。
2つ目に、最適化の問題があります。以下に、
g++ (Ubuntu 4.8.4-2ubuntu1~14.04) 4.8.4で実行した結果を示します。他のシステムやコンパイラでは少し違うかも知れません。
まず、
#include "interval.hpp"
int main()
{
interval x, y;
std::cout.precision(17);
x = 1.;
y = 10.;
std::cout << x / y << "\n";
}
のようなプログラムを、オプション無しでコンパイルして実行すると、[0.099999999999999992,0.10000000000000001]となりますが、これを-O3をつけて最適化して実行すると、
[0.10000000000000001,0.10000000000000001]のように誤った結果が得られてしまいます。驚いたことに、最終行を2回実行するだけの
#include "interval.hpp"
int main()
{
interval x, y;
std::cout.precision(17);
x = 1.;
y = 10.;
std::cout << x / y << "\n";
std::cout << x / y << "\n";
}
だと、最適化を行っても行わなくても[0.099999999999999992,0.10000000000000001] [0.099999999999999992,0.10000000000000001]のように正しい結果が得られました。これは推測なのですが、これはコンパイラがプログラムを解析して、実行結果が分かりきっている部分を予めコンパイル時に計算しておく、いわゆる定数最適化のせいだと考えられます。コンパイル時には丸めの方向は変更されませんので、結果がおかしくなってしまいます。プログラムが単純だと除算の結果が常に1/10だと気づくが、ある程度以上複雑になるとそれに気づかないのでしょうか。最適化の影響を抑えるのはかなり難しい問題なのですが、この場合はその計算に関連する変数にvolatile修飾子を付けることがよく行われます。volatile修飾子は、その変数がいつどんなタイミングで変更されるか分からない特殊な変数であることをコンパイラに伝えるもので、コンパイラはその変数に関係する部分の最適化を抑制します。
第3の問題は、加減乗除と平方根以外の、例えばexpやsinなどの関数の問題です。sinやcosなどの非単調な関数の区間演算は元々難しいのですが、単調増加関数であるexpの区間演算を行うとき、sqrtと同様に
friend interval exp(const interval& x) {
interval r;
fesetround(FE_DOWNWARD);
r.inf = exp(x.inf);
fesetround(FE_UPWARD);
r.sup = exp(x.sup);
fesetround(FE_TONEAREST);
return r;
}
と書いて良いのか? 結論から言うと、このように書いてはいけません。IEEE754の方向付き丸めは、加減乗除と平方根のみに有効で、それ以外の演算に関してはどんな結果になるか保証されていません。そもそも加減乗除と平方根以外の演算では計算結果の精度の保証が何もありませんので、丸めの方向の指定どころの話ではありません。実際調べた結果がここにあります。すなわち、加減乗除と平方根以外の演算に関して区間演算を行うには、その関数を自力で実装する必要があります。以上により、区間演算をきちんと実装するにはかなり手間がかかることが分かります。以下に、これら3つの問題点を解決した実装例を挙げておきます。
([2015/12/23]より、kvライブラリの一部に含めました。)これは、上で挙げた3つ問題点を全て解決したものです。
- 文字列とdoubleの相互変換に関しては、dtostring, stringtodという関数で丸め方向の指定が出来る相互変換を実装しました。これにより、例えば x = "0.1"; と書けばxは0.1を含む区間になりますし、表示は必ず外側に丸めて行われます。
- 最適化対策は、例えば上向き丸めの加算は
static double add_up(const double& x, const double& y) { volatile double r, x1 = x, y1 = y; fesetround(FE_UPWARD); r = x1 + y1; fesetround(FE_TONEAREST); return r; }のように加算に関係する全ての変数をvolatileとしています。 - 数学関数は全て自力で実装しています。実装方法についてはここを参照して下さい。