sqrt，完美正方形和浮点错误
在大多数语言的sqrt函数中（尽管这里我最感兴趣的是C和Haskell），是否有任何保证完美正方形的平方根将被准确返回？ 例如，如果我执行sqrt(81.0) == 9.0 ，那是否安全或者sqrt是否有可能返回8.999999998或9.00000003？
如果无法保证数值精度，那么检查数字是否为完美正方形的首选方法是什么？ 取平方根，拿地板和天花板，确保它们回到原始数字？
谢谢！
在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); }
我需要依赖于dd的空asm volatile块，否则你的编译器可能很聪明并且“优化”了dd的计算。
我使用了fenv.h的一些奇怪的函数，即feclearexcept和fetestexcept 。 查看他们的man页可能是个好主意。
您可能能够工作的另一个策略是计算平方根，检查它是否在尾数的低26位中设置了位，并且如果有，则进行抱怨。 我尝试以下方法。
我需要检查d是否为零，否则它可以为-0.0返回true 。
编辑 ：Eric Postpischil建议用尾数进行黑客攻击可能会更好。 鉴于上面的issquare在另一个流行的编译器中不起作用，我倾向于同意。 我认为以下代码有效：
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位的效果。 首先，当这些位清楚时，我们会准确地回过头来。
根据本文讨论了certificateIEEE浮点平方根的正确性：
IEEE-754二进制浮点运算标准[1]要求分频或平方根运算的结果计算为无限精度，然后舍入到指定精度的两个最接近的浮点数之一围绕着无限精确的结果
由于可以在浮点中精确表示的完美正方形是整数，并且其平方根是可以精确表示的整数，因此完美正方形的平方根应始终完全正确。
当然，不能保证您的代码将使用符合IEEE的浮点库执行。
@tmyklebu完美地回答了这个问题。 作为补充，让我们看看一个可能效率较低的替代品，用于在没有asm指令的情况下测试完美的分数平方。
假设我们有一个符合IEEE 754标准的sqrt，它可以正确地舍入结果。
假设已经处理了exception值（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 ，几乎不完美的正方形任何分数......
所以我们可以结合两个测试：
#include 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)或使用sqared_significand_fits_in_precision中的脏位黑客...
而不是做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 */