9

在sqrt大多数语言的函数中（尽管我对 C 和 Haskell 最感兴趣），是否有任何保证可以准确返回完美平方的平方根？例如，如果我这样做sqrt(81.0) == 9.0，那是安全的还是有机会sqrt返回 8.999999998 或 9.00000003？

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

谢谢！

4

4 回答 4

12

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

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

加法和减法具有擦除尾数低 26 位的效果67108864.0。a我们会a在这些位一开始就清楚的时候回来。

于 2013-01-04T07:28:31.097 回答
6

根据这篇论文，讨论了证明 IEEE 浮点平方根的正确性：

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

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

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

于 2013-01-04T07:32:49.240 回答
1

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

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

If n > 1+floor(p/2)where pis floating point precision (eg p=53 and n>27 in double precision)
Then 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，几乎没有任何部分的完全平方...

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

#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 双精度值小于 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 回答