寺田さん,
平鍋です.
> さて,今回は浮動小数点演算/幾何演算のテストについて質問に参りました.
来ましたね,この話題.
> 私は (3次元)CAD 屋です.従って,普段書いているプログラムはウェブサービ
> ス等ではなく,浮動小数点演算が頻出する幾何計算プログラムです.この浮動小
> 数点演算というヤツは整数演算と比べるとなかなかの曲者のように思います.浮
> 動小数点演算をテストする場合には,整数演算とは違うコツがいるのではないか
> と思い,質問に来ました.
> ...
> 外積 outer_product() のテスト:
> Vector a = ... ; // 乱数を使ってランダムなベクトルを生成
> Vector b = ... ; // 乱数を使ってランダムなベクトルを生成
> Vector c = outer_product( a, b ) ; // 外積を計算
> CPPUNIT_ASSERT( inner_product( a, c ) == 0.0 ) ; // ダメ!
> CPPUNIT_ASSERT( inner_product( b, c ) == 0.0 ) ; // ダメ!
> このテストは,outer_product() が正しく書かれていても,殆ど失敗します.何
> 故なら,浮動小数点演算は誤差を含むため,a と c の内積がピッタリ 0 になる
> ことは少ないからです.
> 従って,アサートを次のように書き換えることになります.
> CPPUNIT_ASSERT( fabs(inner_product( a, c )) < 1.0e-4 ) ;
> CPPUNIT_ASSERT( fabs(inner_product( b, c )) < 1.0e-4 ) ;
> この手のテストの成否は,外積演算の入力 a, b の値に依存します.私の作った
> テストは a, b を乱数で生成していますので,テストの実行のたびに,このテス
> トは成功したり失敗したりします.
> なんとも後味の悪いテストですが,どうにかならないものでしょうか.
私もこの問題については,おそらく5年以上悩んでいます(その後,
5年以上忘れていましたが).そして,解が見つけられていません.
より問題を単純化して,
double mul(double x, double y) { return x*y; }
double div(double x, double y) { return x/y; }
double a = ...;
double b = ...;
double c = mul(a, b);
assert(a == div(c, b));
assert(b == div(c, a));
でさえ,assert が失敗します.これを,どれくらいのepsilon
で押さえるか.
double esp = ...;
assert(fabs(mul(a, div(c, b))) < eps);
assert(fabs(mul(b, div(c, a))) < eps);
eps a, b に依存します.これをやるには,IEEE 754 の仕様から,
* と / の誤差マップ( eps = f(a, b) となる f )を作る必要があ
ります.先に寺田さんが提示された問題に拡張すると,やはり,
nner_product と outer_product の誤差マップを
inner_product/outer_product の「仕様として」提示しなければ正
確なテストはできない.ということになります.これを,全幾何関
数について作成するのは,(論理的に可能ですが)至難と思います.
# こういうアプローチ(IEEE 754 仕様からボトムアップで幾何計算
# 関数の誤差を仕様として定義していく)についての文献,ご存じ
# ないですか?
ところで,XP 的に考えて,「どうせ全ケースのテストはできない
から,アルゴリズムが正しいことだけでもテストしよう」という姿
勢に立てば,上記の誤差マップを安全サイドに単純化することが可
能だと思います.これなら,実用可能かも.たとえば,a, b を
1e-6 〜 1e+6 に限って,eps = min(a, b, 1e-6) でテストする,
とか.
CAD でよく問題になる「近傍」の扱いも同じですよね.結局,ある
近傍を「等しい」と見なす,という考え方が,ブロークンコンセプ
トなんだ,とも思います.つまり,近傍の近傍は近傍だ.という論
法を三段論法で重ねると,有限範囲のすべての点は,近傍になって
しまう,ということ.
また便乗質問になってしまいますが,そういえば,以前,こういう
double と epsilon を使った形状処理ではなく,ディスクリートな
同一判定と多倍長数を使った形状処理の論文(本?)を読んだことが
ありますが,失念しました.寺田さんご存知ないでしょうか? これ
を使った形状処理アプリケーションは,世の中にないのでしょうか
ね?
以上