sqrt,完美正方形和浮点错误

Question

在sqrt大多数语言的功能中(虽然这里我最感兴趣的是C和Haskell),是否有任何保证完美正方形的平方根将被准确返回？例如,如果我这样做sqrt(81.0) == 9.0,那是安全还是有机会sqrt返回8.999999998或9.00000003？

如果无法保证数值精度,那么检查数字是否为完美正方形的首选方法是什么？取平方根,拿地板和天花板,确保它们回到原始数字？

谢谢!

Answer 1

在IEEE 754浮点中,如果双精度值x是非负表示数y的平方(即y*y == x且y*y的计算不涉及任何舍入,上溢或下溢),然后sqrt(x)将返回y.

这都是因为要求sqrt通过IEEE 754标准正确舍入.也就是说,对于任何 x ,sqrt(x)将是与x的实际平方根最接近的两倍.那个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.

编辑: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在这些位清楚的时候准确回来.

Answer 2

根据本文讨论了证明IEEE浮点平方根的正确性:

IEEE-754二进制浮点运算标准[1]要求分频或平方根运算的结果计算为无限精度,然后舍入到指定精度的两个最接近的浮点数之一围绕着无限精确的结果

由于可以在浮点中精确表示的完美平方是整数,并且其平方根是可以精确表示的整数,因此完美平方的平方根应始终完全正确.

当然,不能保证您的代码将使用符合IEEE的浮点库执行.