sqrt
函数(这里主要关注C和Haskell),它是否保证返回一个完全平方数的平方根是精确的?例如,如果我执行sqrt(81.0) == 9.0
,这样安全吗?或者sqrt
会返回8.999999998或9.00000003?如果不能保证数值精度,那么如何检查一个数字是完全平方数呢?取平方根,获取最小整数和最大整数,确保它们的平方等于原始数字吗?
谢谢!
sqrt
函数(这里主要关注C和Haskell),它是否保证返回一个完全平方数的平方根是精确的?例如,如果我执行sqrt(81.0) == 9.0
,这样安全吗?或者sqrt
会返回8.999999998或9.00000003?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
是否为零,因为否则它可能会返回 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;
}
}
从a
中加减67108864.0
将擦除尾数的低26位。只有当这些位在一开始就被清除时,我们才能精确地得到a
。
1/d
会引发异常,可以通过使用signbit(d)
来避免。当然,C实现必须支持asm
扩展。 - Eric Postpischilasm
是不是不必要的?在 double dd = sqrt(d);
之后有一个序列点,所以在调用 fetestexcept
之前,该操作的所有副作用应该已经完成了。 - Eric Postpischilasm
块不是用于排序的;我试图告诉编译器 dd
并没有死亡,最好不要将其优化掉。无论是否执行了 #pragma STDC FENV_ACCESS ON
。无论是 gcc
还是 clang
,似乎都将 sqrt
视为没有副作用的纯函数。(我的 clang
版本甚至不在意依赖于 dd
的 asm
块,并完全省略了它的计算。我猜,如果您想要一个可移植的实现,则确实需要进行位操作。叹气。) - tmyklebusqrt(x)
的低26位清除并不能保证平方根是精确的。 - Eric Postpischil@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
,几乎不是任何分数的完全平方数...
因此,我们可以结合两个测试:
#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 的 dirty bit hacks...
ilogb
。 - tmyklebu不要使用 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 */
f
是一个精确的数字,那么特定的操作是否会产生精确的结果?” - Gabe