どの浮動小数点比較がより正確で、なぜですか？

Question

私は、平方根を計算するためのニュートン法のさまざまな実装を試しています。 1つの重要な決定は、いつアルゴリズムを終了するかということです。

明らかにそれがxの値が大きいことは十分との平方根を表すために可能ではないかもしれないので、y*yとyはxの平方根の現在の推定値であるxの差の絶対値を使用するようにしないだろう精度。

だから、私は相対的な基準を使うことになっています。純粋に私はこのようなものを使用したでしょう：

static int sqrt_good_enough(float x, float y) { 
    return fabsf(y*y - x)/x < EPS; 
} 

これはうまくいくようです。しかし、最近、私はカーニハンとPlaugerのプログラミングスタイルの要素を読み始めており、それらは、第1章では、同じアルゴリズムのためのFortranプログラムを与え、その終了基準、Cに翻訳され、次のようになります。

static int sqrt_good_enough(float x, float y) { 
    return fabsf(x/y - y) < EPS * y; 
} 

両方数学的には等価である ですが、一方の形式を他方よりも優先する理由はありますか？

Answer 1

実際にfabsf（y * y - x）/（y * y）< EPSを最初に書くのでない限り、2つは数学的には全く同じではありません。

しかし、ここでの表現は、ニュートンの反復でyを計算するための式に一致させることです。たとえば、yの式がy =（y + x/y）/ 2の場合、KernighanとPlaugerのスタイルを使用する必要があります。 y =（y * y + x）/（2 * y）の場合は、（y * y - x）/（y * y）< EPSを使用してください。

一般に、終了基準は、abs（y（n + 1）-y（n））が十分に小さい（すなわち、y（n + 1）* EPSよりも小さい）ことでなければならない。これが2つの式が一致する理由です。それらが正確に一致しない場合、スケーリングが異なるために、y（n）の差が浮動小数点誤差よりも小さい間に、残差が十分に小さくないと終端テストが判断する可能性があります。 y（n）が変化しなくなり、終了基準が決して満たされないため、結果は無限ループになります。例えば

次のMATLABコードは、あなたの最初の例として、ソルバーまったく同じニュートンであるが、それは永遠に実行されます。

x = 6.800000000000002 
yprev = 0 
y = 2 
while abs(y*y - x) > eps*abs(y*y) 
    yprev = y; 
    y = 0.5*(y + x/y); 
end 

それのC/C++のバージョンが同じ問題を抱えています。

Answer 2

これらはまだ同等ではありません。下の方程式は数学的にはfabsf(y*y - x)/(y*y) < EPSに相当します。あなたが見ている問題は、y*yがオーバーフローした場合（おそらくxがFLT_MAXとyが不都合に選択されているため）、終了が発生しない可能性があります。次のインタラクションでは、倍精度を使用します。

>>> import math 
>>> x = (2.0 - 2.0 ** -52) * 2.0 ** 1023 
>>> y = x/math.sqrt(x) 
>>> y * y - x 
inf 
>>> y == 0.5 * (y + x/y) 
True 

EDIT：コメント（今削除されました）が指摘されているので、繰り返しと終了テストの間で操作を共有することもできます。

EDIT2：どちらもサブ正規の問題xを持っている可能性があります。 professionalsは、両極端の複雑さを避けるためにxを正規化します。