9

ほとんどの言語のsqrt関数では (ここでは主に C と Haskell に興味があります)、完全平方の平方根が正確に返されるという保証はありますか? たとえば、私がそうする場合、それsqrt(81.0) == 9.0は安全sqrtですか、それとも 8.999999998 または 9.00000003 を返す可能性はありますか?

数値の精度が保証されていない場合、数値が完全な二乗であることを確認するための推奨される方法は何ですか? 平方根を取り、床と天井を取得し、それらが元の数に戻っていることを確認しますか?

ありがとうございました！

4

4 に答える 4

12

IEEE 754 浮動小数点では、倍精度値 x が非負の表現可能な数値 y の 2 乗である場合 (つまり、y*y == x であり、y*y の計算に丸め、オーバーフロー、またはアンダーフローが含まれない場合)、 sqrt(x) は y を返します。

これはすべて、IEEE 754 標準で sqrt を正しく丸める必要があるためです。つまり、任意のx の sqrt(x) は、x の実際の平方根に最も近い double になります。sqrt が完全な正方形で機能することは、この事実の単純な結果です。

浮動小数点数が完全な二乗であるかどうかを確認したい場合、私が考えることができる最も単純なコードは次のとおりです。

int issquare(double d) {
  if (signbit(d)) return false;
  feclearexcept(FE_INEXACT);
  double dd = sqrt(d);
  asm volatile("" : "+x"(dd));
  return !fetestexcept(FE_INEXACT);
}

asm volatileに依存する空のブロックが必要ですdd。そうしないと、コンパイラが賢く、 の計算を「最適化」する可能性があるためですdd。

からいくつかの奇妙な関数fenv.h、つまりfeclearexceptとを使用しましたfetestexcept。彼らのmanページを見るのはおそらく良い考えです。

うまくいく可能性のある別の戦略は、平方根を計算し、仮数の下位 26 ビットにビットが設定されているかどうかを確認し、設定されている場合は文句を言うことです。以下でこのアプローチを試します。

そして、dがゼロかどうかを確認する必要がありtrueました-0.0。

EDIT : Eric Postpischil は、仮数部をハックする方が良いかもしれないと提案しました。issquare上記が別の一般的なコンパイラで機能しないことを考えるとclang、私は同意する傾向があります。次のコードが機能すると思います。

int _issquare2(double d) {
  if (signbit(d)) return 0;
  int foo;
  double s = sqrt(d);
  double a = frexp(s, &foo);
  frexp(d, &foo);
  if (foo & 1) {
    return (a + 33554432.0) - 33554432.0 == a && s*s == d;
  } else {
    return (a + 67108864.0) - 67108864.0 == a;
  }
}

加算および減算67108864.0はa、仮数の下位 26 ビットを消去する効果があります。aそれらのビットが最初にクリアされた正確な時刻に戻ります。

于 2013-01-04T07:28:31.097 に答える
6

この論文によると、IEEE 浮動小数点平方根の正確性の証明について説明しています。

The IEEE-754 Standard for Binary Floating-Point Arithmetic [1] requires that the result of a divide or square root operation be calculated as if in infinite precision, and then rounded to one of the two nearest floating-point numbers of the specified precision that surround the infinitely precise result

Since a perfect square that can be represented exactly in floating-point is an integer and its square root is an integer that can be precisely represented, the square root of a perfect square should always be exactly correct.

Of course, there's no guarantee that your code will execute with a conforming IEEE floating-point library.

于 2013-01-04T07:32:49.240 に答える
1

@tmyklebu は質問に完全に答えました。補足として、asm ディレクティブを使用せずに分数の完全な 2 乗をテストするための、おそらく効率の悪い代替手段を見てみましょう。

結果を正しく丸める IEEE 754 準拠の sqrt があるとします。
例外的な値 (Inf/Nan) とゼロ (+/-) が既に処理されているとします。が奇数である場所に分解し
ましょう。 そして、n ビットにまたがる場所: .sqrt(x)I*2^mI
I1+2^(n-1) <= I < 2^n

n > 1+floor(p/2)whereが浮動小数点精度の場合p(例: 倍精度で p=53 および n>27)
Then 2^(2n-2) < I^2 < 2^2n. 奇数である
ため、奇数でもあり、したがって > p ビットにまたがります。 したがって、この精度で表現可能な浮動小数点の正確な平方根ではありません。II^2
I

しかし、 が与えられた場合、それは完全な正方形であるI^2<2^pと言えますか? 答えは明らかにノーです。テイラー展開はx

sqrt(I^2+e)=I*(1+e/2I - e^2/4I^2 + O(e^3/I^3))

したがって、平方根までは ... に正しく丸められe=ulp(I^2)ます(最も近い偶数に丸める、切り捨てる、またはフロア モード)。sqrt(ulp(I^2))rsqrt(I^2+e)=I

したがって、それを主張する必要がありsqrt(x)*sqrt(x) == xます。
But above test is not sufficient, for example, assuming IEEE 754 double precision, sqrt(1.0e200)*sqrt(1.0e200)=1.0e200, where 1.0e200 is exactly 99999999999999996973312221251036165947450327545502362648241750950346848435554075534196338404706251868027512415973882408182135734368278484639385041047239877871023591066789981811181813306167128854888448 whose first prime factor is 2^613, hardly a perfect square of any fraction...

したがって、両方のテストを組み合わせることができます。

#include <float.h>
bool is_perfect_square(double x) {
    return sqrt(x)*sqrt(x) == x
        && squared_significand_fits_in_precision(sqrt(x));
}
bool squared_significand_fits_in_precision(double x) {
    double scaled=scalb( x , DBL_MANT_DIG/2-ilogb(x));
    return scaled == floor(scaled)
        && (scalb(scaled,-1)==floor(scalb(scaled,-1)) /* scaled is even */
            || scaled < scalb( sqrt((double) FLT_RADIX) , DBL_MANT_DIG/2 + 1));
}

編集: 整数の場合に制限したい場合は、それを確認するfloor(sqrt(x))==sqrt(x)か、squared_significand_fits_in_precision でダーティ ビット ハックを使用することもできます...

于 2013-01-04T23:14:59.040 に答える
0

する代わりにsqrt(81.0) == 9.0、試してみてください9.0*9.0 == 81.0。これは、平方が浮動小数点の大きさの制限内にある限り常に機能します。

編集：「浮動小数点の大きさ」が何を意味するのか、おそらく不明確でした。私が言いたいのは、IEEE double の場合は 2**53 未満で、精度を失うことなく保持できる整数値の範囲内に数値を維持することです。また、平方根が整数であることを確認するための別の操作があることも期待していました。

double root = floor(sqrt(x) + 0.5);  /* rounded result to nearest integer */
if (root*root == x && x < 9007199254740992.0)
    /* it's a perfect square */
于 2013-01-04T05:50:54.847 に答える