2016/03/16(水)kv-0.4.29
細かいバグフィックスが主で、大きな変更や機能追加はありません。
小野 寛生様のご指摘により、conv-double.hpp, conv-dd.hppで配列の外側にアクセスする可能性のある危険な書き方を直しました。ご指摘ありがとうございました。また、Visual C++では丸めの向きを変更するのに_controlfpを用いていますが、_controlfp_sが推奨されている点もご指摘いただき、そのように直しました。また、ついでにVisual C++ 2013でmpfrを用いたものを除く全てのサンプルファイルがコンパイル可能なことを確認しました。
Visual C++での丸めの向きの変更に関しては、非標準である_controlfpではなくVisual C++ 2013で追加された標準規格のfesetroundを使いたいとずっと思っていました。しかし、最近上向き丸めと下向き丸めが逆になっているというとんでもないバグを抱えていることが分かり、当分これは使えないと判断しました。
fesetroundの丸めモードが逆になる
困ったものです。interval-simpleというkvの区間演算部分だけを切り出した簡単なライブラリを公開していますが、こちらは単純化のためfesetroundを用いた実装をしていたのですが、丸めの向きが逆では致命的過ぎるのでifdefで分けてVisual C++のときは_controlfp_sを使うように直しました。
(2016/4/4に追記)
Visual Studio 2015 update2 が出たので試してみましたが、同じ結果でした。x64では相変わらず上向き丸め下向き丸めが逆になります。2015/12/24(木)日本応用数理学会 三部会連携「応用数理セミナー」
「精度保証付き数値計算プログラムの実装について」(ap-seminar.pdf)
中身の大半はkvライブラリの紹介です。これで少しでもユーザが増えるといいなあ。
2015/12/23(水)kv-0.4.28
今回はライブラリ本体の変更はほとんどなく、highderiv.hpp(高階微分の計算)、test-rounding.cc(丸めモードの変更が正常に動いているか簡易チェック)を追加したくらいです。webのドキュメントは、psaやaffineにいくらか加筆しています。明日、応用数理セミナーの講演でこのライブラリの話をするのでその前に区切りとして全ての変更を放出しておきたかったという理由もあります。
test-rounding.ccの追加は、前回日記で触れたcygwinのsqrt問題があったので、とりあえず軽くチェックしようと思ったためです。いろいろな環境でテストしてみたところ、概ね以下のような感じです。
- Linux, Macは64bit OSで64bitコンパイルしていれば問題なし
- x86の32bitモードはFPUが使われてしまうためIEEE754に従っておらず、お勧め出来ない。
- Windowsでは、cygwin, msys2で最適化をしないと、sqrtに丸めの向きが変わらない問題が発生。VC++の64bitなら問題なし。
2015/12/16(水)cygwinのsqrtの丸めモード変更
問題は、cygwinのsqrtで丸めの向きが最適化をかけないときに変わらないというものです。環境は64bitのwindows7で、cygwinも64bitの最新のものです。gccは4.9.3でした。
#include <stdio.h>
#include <math.h>
#include <fenv.h>
int main()
{
volatile double x, y, z;
x = 2.;
fesetround(FE_DOWNWARD);
y = sqrt(x);
fesetround(FE_UPWARD);
z = sqrt(x);
fesetround(FE_TONEAREST);
if (y == z) {
printf("sqrt error\n");
}
}
このようなプログラムで、cc cygwin-sqrt.cのように最適化オプションをつけないでコンパイルしたとき、丸めの向きが上向きでも下向きでも結果が同じになってしまいます。
cc -O3 cygwin-sqrt.cのように最適化をかければ問題ありません。いろいろ試すと、-O0ではアウト、-O1, -O2, -O3はセーフでした。
volatileを使わない場合最適化をかけることで結果がおかしくなるのはさほど珍しくありませんが、このように最適化をかけない状態で結果が異常になるのはとても珍しいです。
cc -Sでアセンブリソースを表示させてみると、最適化無しでは単に外部のsqrt関数を
call sqrt
と呼んでいるのに対して、-O3だと
sqrtsd %xmm0, %xmm1
のようにSSE2で直接計算するコードが生成されているので、どうやら外部(libm?)のsqrtに問題がありそうです。以前にVisual C++で似たようなことがあったので、それとも関係あるのでしょうか。どなたか詳しい方に教えて欲しいものです。
これで丸めの向きが変わらないのはIEEE754的にまずいので、そんな腐った環境は無視するというのも一つの考え方でしょうが、こういう事例が多くあるようならsqrtに関しては丸めモード変更に期待しないで何とかすることを検討するべきかも。kvライブラリでいえば-DKV_NOHWROUNDを付ければ問題なくなるので、それをデフォルトにするとか。
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としています。 - 数学関数は全て自力で実装しています。実装方法についてはここを参照して下さい。