浮動小数点数に対するテスト

どうやってテストするべきなのかよく判らないもののひとつに浮動小数点数がある。どこをどうやってテストしていいかがそもそも全然判らないものよりはマシかもしれないけれど。
複雑な数値計算で精度がすごく重要とか場合でなければ、多少の誤差はあまり気にしなくても良いのだろうか。たとえば

$\sin^2 \theta = \frac{1}{2}(1 - \cos (2\theta))$

という等式がある。この式でθが0に近い場合、cos(2θ)≒1なので右辺を計算すると桁落ちが起こる(一方、左辺はほぼ正確な値と考えて良いはず)。

gosh> (define (lhs theta)
        (expt (sin theta) 2))
lhs
gosh> (define (rhs theta)
        (/ (- 1.0 (cos (* 2.0 theta))) 2.0))
rhs
gosh> (lhs 0.001)
9.99999666666711e-7
gosh> (rhs 0.001)
9.999996666842925e-7

(本題に関係ないけど、aのb乗をSchemeでは(expt a b)と書くのは何か違和感がある)。
あまり違いがないように見えるので相対誤差を調べる。xとyの相対誤差を $\frac{|x-y|}{{\rm max}(|x|, |y|)}$ で測ることにすると

gosh> (define (relative-error x y)
        (/ (abs (- x y)) (max (abs x) (abs y))))
relative-error
gosh> (relative-error (lhs 0.001) (rhs 0.001))
1.758144522985574e-11

となり、やはり誤差は小さいように見える。でも計算機イプシロン=2^-52との比をとると

gosh> (define m-eps (expt 2.0 -52))
m-eps
gosh> m-eps
2.220446049250313e-16
gosh> (/ (relative-error (lhs 0.001) (rhs 0.001)) m-eps)
79179.7901858131

となり、今度はかなり大きな誤差に思えてくる(計算機イプシロンの意味は「計算機イプシロンのこと」で書いたので省略)。
次に、右辺の $\frac{1}{2}(1 - \cos (2\theta))$ を間違えて $\frac{1}{2}(1 - \cos (\theta))$ にした場合を考える。

gosh> (define (wrong-rhs theta)
        (/ (- 1.0 (cos theta)) 2.0))
wrong-rhs
gosh> (wrong-rhs 0.001)
2.4999997916275163e-7
gosh> (relative-error (lhs 0.001) (wrong-rhs 0.001))
0.7499999375039053
gosh> (/ (relative-error (lhs 0.001) (wrong-rhs 0.001)) m-eps)
3.377699439070483e15

そもそも式が間違っているので、さっきと比べて誤差がはるかに大きい。
別の例。方程式

x²-1000.001x + 1 = 0

は、(x - 1000)(x - 0.001)=0と因数分解できるので、解は1000と0.001。
でも解の公式 $\frac{-b \pm \sqrt{b^2 -4ac}}{2a}$ を使って解を求めると桁落ちが起こる。

gosh> (define (solve-quadratic-equation a b c)
        (let ((root-D (sqrt (- (expt b 2) (* 4.0 a c)))))
          (list (/ (+ (- b) root-D) (* 2.0 a))
                (/ (- (- b) root-D) (* 2.0 a)))))
solve-quadratic-equation
gosh> (define x1 (list-ref (solve-quadratic-equation 1.0 -1000.001 1.0) 1))
x1
gosh> x1
9.999999999763531e-4
gosh> (relative-error 0.001 x1)
2.3646883062777846e-11
gosh> (/ (relative-error 0.001 x1) m-eps)
106496.09375

解く方程式を

x²-100000.00001x + 1 = (x - 100000)(x - 0.00001) = 0

にすると、桁落ちによる誤差がさらに大きくなる。

gosh> (define x2 (list-ref (solve-quadratic-equation 1.0 -100000.00001 1.0) 1))
x2
gosh> x2
1.0000003385357559e-5
gosh> (relative-error 0.00001 x2)
3.385356411845042e-7
gosh> (/ (relative-error 0.00001 x2) m-eps)
1.524628987490165e9

相対誤差の許容値を変えてテストしてみると、次のような感じになる。

gosh> (define (make-float-eq . opt)
        (let ((tolerance (get-optional opt 1e-15)))
          (lambda (x y) (<= (relative-error x y) tolerance))))
make-float-eq
gosh> (use gauche.test)
#<undef>
gosh> (test* "" 0.001 x1 (make-float-eq))
test , expects 0.001 ==> ERROR: GOT 9.999999999763531e-4
#<undef>
gosh> (test* "" 0.001 x1 (make-float-eq 1e-10))
test , expects 0.001 ==> ok
#<undef>
gosh> (test* "" 0.00001 x2 (make-float-eq))
test , expects 1.0e-5 ==> ERROR: GOT 1.0000003385357559e-5
#<undef>
gosh> (test* "" 0.00001 x2 (make-float-eq 1e-10))
test , expects 1.0e-5 ==> ERROR: GOT 1.0000003385357559e-5
#<undef>
gosh> (test* "" 0.00001 x2 (make-float-eq 1e-5))
test , expects 1.0e-5 ==> ok
#<undef>

調べてみると、Ruby1.9のMiniTest::Assertionsのassert_in_epsilonが、相対誤差で比較していて(ただし相対誤差の定義がここで使ったのとは少し違い、二つの数のmaxではなくminを取っている)、許容値のデフォルト値は0.001になっていた。許容値は意外と大きくても良いのだろうか。
しかしassert_in_epsilonにバグがあった。絶対値を取っていないので負数のときおかしくなる。

irb(main):001:0> require "minitest/unit"
=> true
irb(main):002:0> include MiniTest::Assertions
=> Object
irb(main):003:0> assert_in_epsilon(-100, -100)
MiniTest::Assertion: Expected -100 - -100 (0) to be < -0.1.