sqrt、完美平方和浮点错误答案

【问题标题】：sqrt, perfect squares and floating point errorssqrt、完美平方和浮点错误
【发布时间】：2012-12-18 13:44:42
【问题描述】：

在大多数语言的sqrt 函数中（虽然我对C 和Haskell 最感兴趣），是否有任何保证可以准确返回完美平方的平方根？例如，如果我执行sqrt(81.0) == 9.0，那是否安全，或者sqrt 是否有可能返回 8.999999998 或 9.00000003？

如果不能保证数值精度，那么检查一个数字是否为完美平方的首选方法是什么？取平方根，得到地板和天花板，并确保它们平方回到原来的数字？

谢谢！

【问题讨论】：

是的，不保证数值精度。
What every computer scientist should know about floating point arithmetic
@MitchWheat：众所周知，给定不精确的数字，FP 操作将产生不精确的结果。我认为 OP 专门谈论可以精确表示的精确数字，即使在 FP 中也是如此。他的问题更像是“鉴于f 是一个精确的数字，特定的操作会产生精确的结果吗？”
@Mitch: 9.0 是一个整数，因此它可以精确地表示为 1.001*2^3 形式的二进制浮点数，并且其数量级小于所有 IEEE 格式的尾数，是以任何格式精确表示（0x41100000 作为单精度或 0x4022000000 作为双精度）。
@AlokSave 迷信！平方根是基本的 IEEE 754 之一，必须正确舍入。请参阅 tmyklebu 的回答。

标签： c haskell floating-point

【解决方案1】：

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

这都是因为 IEEE 754 标准要求 sqrt 正确舍入。也就是说，对于 any 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);
}

我需要依赖于dd 的空asm volatile 块，因为否则您的编译器可能会很聪明并“优化”掉dd 的计算。

我使用了来自fenv.h 的几个奇怪的函数，即feclearexcept 和fetestexcept。看看他们的man 页面可能是个好主意。

另一种可能可行的策略是计算平方根，检查它是否在尾数的低 26 位中设置了位，如果有，请抱怨。我在下面尝试这种方法。

我需要检查d 是否为零，否则它可以为-0.0 返回true。

编辑：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;
  }
}

从a 中加减67108864.0 具有擦除尾数的低26 位的效果。当这些位一开始就被清除时，我们将得到a。

【讨论】：

在可以说该代码实现“Is d a square?”之前，必须满足一些要求。您需要从 C 到 IEEE 754 或其他合适的浮点规范的绑定，这是 C 标准不需要的，许多 C 实现也不提供。您必须包含并使用适当的#pragma 将 FENV_ACCESS 设置为打开，或者知道它对于正在使用的 C 实现是打开的。 1/d 引发异常，可以通过使用 signbit(d) 来避免。当然，C 实现必须支持asm 扩展。
在某些处理器上使用 feclearexcept 和 fetestexcept 访问浮点环境可能非常慢。检查有效数（不是尾数）的低位的策略会更快，因为这可以通过普通算术来完成。上溢、下溢和次正规将不是问题，因为平方根远离浮点范围的极端。但是，我会担心使用的精度。我不清楚 C 附件 F 是否会阻止实现使用比标称类型更高的精度，这是在没有附件 F 的情况下允许的。
不应该asm 是不必要的吗？在double dd = sqrt(d); 之后有一个序列点，因此所有副作用都应该在调用fetestexcept 之前完成。
asm 块不用于排序；我试图告诉 copmiler dd 没有死，最好不要优化它。不管我是否完成了#pragma STDC FENV_ACCESS ON。 gcc 和 clang 似乎都将 sqrt 视为没有副作用的纯函数。（我的clang 副本甚至不在乎有一个asm 块取决于dd；它完全省略了它的计算。我想如果你想要一个可移植的实现，你确实需要做一些位摆弄。叹息.)
啊。在这种情况下，正确的实现不会优化 dd 的初始化，因为它包含 fetestexcept 访问的副作用。不正确的实现可能会，但asm 只是针对特定不正确实现的黑客攻击。

【解决方案2】：

根据this paper，讨论证明IEEE浮点平方根的正确性：

二进制浮点的 IEEE-754 标准算术 [1] 要求除法或平方的结果根运算被计算为无限精度，并且然后四舍五入到两个最近的浮点数之一包围的指定精度的数字无限精确的结果

由于可以用浮点数精确表示的完全平方是一个整数，而它的平方根是一个可以精确表示的整数，所以完全平方的平方根应该总是完全正确的。

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

【讨论】：

假设 sqrt 例程的实现类似于长除法（重复移位和减法），算法以零余数结束并产生精确结果。

【解决方案3】：

@tmyklebu 完美地回答了这个问题。作为补充，让我们看看在没有 asm 指令的情况下测试分数的完全平方的效率可能较低的替代方法。

假设我们有一个符合 IEEE 754 的 sqrt，它可以正确地对结果进行四舍五入。
假设已经处理了异常值 (Inf/Nan) 和零 (+/-)。
让我们将sqrt(x) 分解为I*2^m，其中I 是一个奇数。
其中I 跨越 n 位：1+2^(n-1) <= I < 2^n。

如果 n > 1+floor(p/2) 其中 p 是浮点精度（例如 p=53 和 n>27 在双精度）
然后2^(2n-2) < I^2 < 2^2n.
由于I 是奇数，I^2 也是奇数，因此跨越 > p 位。
因此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。
但上述试验是不充分的，例如，假设IEEE 754双精度，sqrt(1.0e200)*sqrt(1.0e200)=1.0e200，其中1.0e200正是99999999999999996973312221251036165947450327545502362648241750950346848435554075534196338404706251868027512415973882408182135734368278484639385041047239877871023591066789981811181813306167128854888448其第一素因数是2^613，任何分数... P>的几乎完美的方

所以我们可以结合两个测试：

#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 中使用脏位黑客...

【讨论】：

该死的。希望我以前知道ilogb。
不幸的是，虽然 ilogb 和 scalb 在良好的库中可用，但直到 C++11 才出现在 C 标准中......所以我作弊了。

【解决方案4】：

试试9.0*9.0 == 81.0，而不是sqrt(81.0) == 9.0。只要平方在浮点幅度的范围内，这将始终有效。

编辑：我可能不清楚我所说的“浮点幅度”是什么意思。我的意思是将数字保持在可以保持而没有精度损失的整数值范围内，对于 IEEE 双精度值小于 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 */

【讨论】：

当你用无穷大代替 81.0 和 1e300 代替 9.0 时它会失败。（至少对于双打。）或者不理会 81.0 并使用 -9.0 作为 9.0。
@EricPostpischil，我看到这个答案引起了一些混乱。希望我已经澄清了。
但是 scalb( FLT_RADIX , 2*DBL_MANT_DIG) 也是一个完美的正方形，不是吗？
确实不是，因为它是 2 的奇次幂。可能你的意思是 scalbn(1.0, 2*DBL_MANT_DIG)。