外野がいきなりツッコミ入れるのもどうかと思うので、ここで。
この問題は、 m68k FPE のデバッグをしている時に、libm の exp() のテストがおかしいんじゃないかという疑問を持ったことに端を発する。
2012/05 時点で、m68k FPE の exp(x) が stub 実装で、常に x を返すという痛々しい事実はまず忘れておいて、そもそもこのテストコードの意味を考えてみる。
テストとは、
「信頼されたもの」を使って、「信頼出来るかわからないもの」が信頼出来るかどうか調べる行為である。
今回問題になってるテストコードは
exp(x + y) - exp(x) * exp(y) < eps
さてここで、「信頼されていないもの=調べたいもの」には .N をつけてみよう。
exp.N(x + y) - exp.N(x) * exp.N(y) < eps
左辺に信頼されていないものを集めて、右辺に信頼出来るものを集める。が、もうこれ以上移項できない。
そうするとこの式を見て分かるように、このテストコードでは、「左辺の式が eps 以下である」ことしかテストできていない。exp 自体が信頼できるかどうかは、このコードでは分からない。
仮に exp 自体がすでに信用出来るものである、とすると、このテストコードがテストすることになるのは、数学領域の「exp(x+y) = exp(x)*exp(y)」という式が数学的に正しいかどうか、であって libm の exp ではない。
このテストコードを、まがりなりにも libm のテストとして成立するように書きなおすとしたら、
#define Ex (計算済みのexp(x)を定数にしたもの)
#define Ey (計算済みのexp(y)を定数にしたもの)
exp(x + y) - Ex * Ey < eps
とすれば、
exp.N(x + y) < eps + Ex * Ey
となって、「exp(x + y) が信頼出来るかどうか」のテストコードとして成立する。
もちろんこうなると、信頼されている Ex * Ey は事前に計算してよいので、
#define Exy (計算済みの exp(x) * exp(y) を定数にしたもの)
exp.N(x + y) - Exy < eps
としてよい。
また、x+y も信頼されているので事前に計算してよい。z=x+y とおけば、
exp.N(z) - Ez < eps
となる。すでにこれは exp[2f]?_product のテストではない。掛け算どこにも出てこない。
以上、libm のテスト exp[2f]?_product は意味が無い、という考察でした。
もしかして、数学の正しさをテストしたかったんならごめん、ライブラリのテストコードじゃなくて他所でやって・・・
