14

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

y*yとの間の絶対差を使用することは明らかに適切xではありません。yの平方根の現在の推定値は です。これはx、 の値が大きいと、xその平方根を十分な精度で表すことができない可能性があるためです。

したがって、相対的な基準を使用することになっています。単純に、私は次のようなものを使用していたでしょう：

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

そして、これは非常にうまく機能しているようです。しかし最近、カーニハンとプラウガーのThe Elements of Programming Styleを読み始め、第 1 章で同じアルゴリズムの Fortran プログラムを提供しています。その終了基準を C に翻訳すると、次のようになります。

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

どちらも数学的には同等ですが、どちらか一方を優先する理由はありますか?

4

2 に答える 2

5

それらはまだ同等ではありません。一番下のものは、数学的には と同等fabsf(y*y - x) / (y*y) < EPSです。あなたの問題は、y*yオーバーフローした場合 (おそらくxisFLT_MAXとyis が不運に選択されたため)、終了が発生しない可能性があることです。次の相互作用は double を使用します。

>>> 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

編集：コメント（現在は削除されています）が指摘したように、反復と終了テストの間で操作を共有することも良いことです。

EDIT2: どちらもおそらく subnormal に問題がありますx。専門家xは、両極端の合併症を避けるために正常化します。

于 2012-03-02T13:25:19.227 に答える
3

最初のものに対して 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 コードは、最初の例とまったく同じ Newton ソルバーですが、永久に実行されます。

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++ バージョンにも同じ問題があります。

于 2012-03-02T14:13:46.310 に答える