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としています。 - 数学関数は全て自力で実装しています。実装方法についてはここを参照して下さい。
2015/11/27(金)区間演算の実装について(1)
区間演算は原理は簡単ですが、実際に実装となると案外ハマリポイントが多く、難しいものです。簡単な区間演算の実装例を通して、少し解説してみようと思います。
使いやすい区間演算を実装するには、区間型を定義してその区間型を組み込み型と同等に扱える、いわゆる演算子多重定義(operator overloading)のできるプログラミング言語が望ましいです。自分はC++を愛用していますので、ここでもC++で実装します。両端点の型はdouble、方向付き丸めはC99のfenv.hとfesetroundを使います。
最も単純で最低限の機能のみの実装を示します (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 std::ostream& operator<<(std::ostream& s, const interval& x) {
s << '[';
s << x.inf;
s << ',';
s << x.sup;
s << ']';
return s;
}
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;
}
const double& lower() const {
return inf;
}
const double& upper() const {
return sup;
}
double& lower() {
return inf;
}
double& upper() {
return sup;
}
};
#endif // INTERVAL_HPP
使える演算は加減乗除と平方根で、coutでの表示にも対応してます。lowerとupperは下限と上限へのアクセサ。乗除算は手抜きせずに一応ちゃんと符号による場合分けを行っています。C++の初心者でも入門書を読みながら実装できるレベルでしょう。次にこのライブラリの使用例を示します (test-interval.cc)。
#include "interval.hpp"
int main()
{
interval x;
interval y = 1.;
interval z(1.);
x = 1.;
y = 10.;
z = x / y;
std::cout << z << "\n";
std::cout.precision(17);
std::cout << z << "\n";
// copy
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";
// operation with constant
std::cout << x + 1 << "\n";
std::cout << x + 1. << "\n";
// compound assignment operator
z += x;
z += 1.;
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";
}
実行すると、[0.1,0.1] [0.099999999999999992,0.10000000000000001] [4,6] [-3,-1] [3,8] [0.25,0.66666666666666674] [2,3] [2,3] 3 4 [3.5,4]のようになります。この2つのファイルを一応アップロードしておきます。
interval-supersimple.tar.gz
これを使って、大きく誤差が出る計算をしてみます。2次方程式ax2+bx+c=0でa=1, b=1015, c=1014として、2つの解のうちの大きい方を、
- 普通に解の公式を使う
- 解の公式の分子を有理化する
まず普通にdoubleで。
#include <iostream>
#include <cmath>
int main()
{
double a, b, c, x, y;
std::cout.precision(17);
a = 1.;
b = 1e15;
c = 1e14;
x = (-b + sqrt(b * b - 4. * a * c)) / (2. * a);
y = 2 * c / (-b - sqrt(b * b - 4. * a * c));
std::cout << x << "\n";
std::cout << y << "\n";
}
のようなプログラムを実行すると、-0.125 -0.10000000000000002のようになります。この計算は後者がほぼ正しく、前者は大きな誤差が入っています。
これを区間演算で計算してみます。プログラムは、
#include <iostream>
#include <cmath>
#include "interval.hpp"
int main()
{
interval a, b, c, x, y;
std::cout.precision(17);
a = 1.;
b = 1e15;
c = 1e14;
x = (-b + sqrt(b * b - 4. * a * c)) / (2. * a);
y = 2 * c / (-b - sqrt(b * b - 4. * a * c));
std::cout << x << "\n";
std::cout << y << "\n";
}
のような感じ。演算子多重定義のおかげで最小限の変更で区間演算が行えます。計算結果は、[-0.1875,-0.0625] [-0.10000000000000003,-0.099999999999999992]となり、前者の区間幅が極端に広いことで、計算結果が怪しいと気づくことが出来ます。
さて、学習用としてはこれでいいのですが、実際に使おうとするとこのレベルの実装ではいろいろと問題点があります。少し記事が長くなったので、問題点は次の記事で。