开方，完美的广场和浮点错误 [英] sqrt, perfect squares and floating point errors

问题描述

在大多数语言中的开方功能（尽管在这里我使用C和Haskell最感兴趣），是否有任何保证一个完美的正方形的平方根会精确地返回？例如，如果我这样做的sqrt（81.0）== 9.0 ，是安全的还是有一个机会，开方将返回8.999999998或9.00000003？

In the sqrt function of most languages (though here I'm mostly interested in C and Haskell), are there any guarantees that the square root of a perfect square will be returned exactly? For example, if I do sqrt(81.0) == 9.0, is that safe or is there a chance that sqrt will return 8.999999998 or 9.00000003?

如果数值precision不能保证，这将是preferred的方法来检查，一些是完全平方？就拿平方根，得到了地板和天花板，并确保它们方回原来的号码？

If numerical precision is not guaranteed, what would be the preferred way to check that a number is a perfect square? Take the square root, get the floor and the ceiling and make sure they square back to the original number?

感谢您！

推荐答案

在IEEE 754浮点，如果双precision值x是一个非负重presentable数y的平方（即ÿ * Y == x和y * y的计算不涉及任何舍入，溢出或下溢），那么SQRT（x）的将返回年。

In IEEE 754 floating-point, if the double-precision value x is the square of a nonnegative representable number y (i.e. y*y == x and the computation of y*y does not involve any rounding, overflow, or underflow), then sqrt(x) will return y.

这一切都是因为需要开方由IEEE 754标准要正确圆润。也就是说，开方（x）的，为的任何的X，将是最接近双x的实际平方根。这开方的作品完美广场是这一事实的简单推论。

This is all because sqrt is required to be correctly-rounded by the IEEE 754 standard. That is, sqrt(x), for any x, will be the closest double to the actual square root of x. That sqrt works for perfect squares is a simple corollary of this fact.

如果您要检查一个浮点数是否是一个完美的广场，这里是最简单的code我能想到的：

If you want to check whether a floating-point number is a perfect square, here's the simplest code I can think of:

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挥发性块，它依赖于 DD ，否则你的编译器可能是聪明，优化走了DD的计算。

I need the empty asm volatile block that depends on dd because otherwise your compiler might be clever and "optimise" away the calculation of dd.

我用了几个奇怪的职能从 fenv.h ，即 feclearexcept 和 fetestexcept 。这可能是一个好主意，看看他们的人页。

I used a couple of weird functions from fenv.h, namely feclearexcept and fetestexcept. It's probably a good idea to look at their man pages.

这你也许可以使工作的另一个策略是计算平方根，检查是否已尾数的低26位设置位，如果确实如此抱怨。我试试下面这种方法。

Another strategy that you might be able to make work is to compute the square root, check whether it has set bits in the low 26 bits of the mantissa, and complain if it does. I try this approach below.

和我需要检查 D 是否为零，否则它可以返回真正为 -0.0 。

And I needed to check whether d is zero because otherwise it can return true for -0.0.

修改：埃里克Postpischil建议，与尾数黑客周围可能会更好。鉴于上述 issquare 没有在另一个受欢迎的编译工作，铛，我倾向于同意。我认为以下code ++工程：

EDIT: Eric Postpischil suggested that hacking around with the mantissa might be better. Given that the above issquare doesn't work in another popular compiler, clang, I tend to agree. I think the following code works:

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 回到什么时候这些位是摆在首位清晰。

Adding and subtracting 67108864.0 from a has the effect of wiping the low 26 bits of the mantissa. We will get a back exactly when those bits were clear in the first place.

这篇关于开方，完美的广场和浮点错误的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！