我正在查看另一个问题(此处),有人正在寻找一种方法来获取 x86 程序集中的 64 位整数的平方根。
事实证明这很简单。解决方案是转换为浮点数,计算 sqrt,然后再转换回来。
我需要在 C 中做一些非常相似的事情,但是当我研究等价物时,我有点卡住了。我只能找到一个接受双打的 sqrt 函数。双精度数不具备存储大型 64 位整数而不引入显着舍入误差的精度。
是否有一个我可以使用的具有long double
sqrt 函数的通用数学库?
我正在查看另一个问题(此处),有人正在寻找一种方法来获取 x86 程序集中的 64 位整数的平方根。
事实证明这很简单。解决方案是转换为浮点数,计算 sqrt,然后再转换回来。
我需要在 C 中做一些非常相似的事情,但是当我研究等价物时,我有点卡住了。我只能找到一个接受双打的 sqrt 函数。双精度数不具备存储大型 64 位整数而不引入显着舍入误差的精度。
是否有一个我可以使用的具有long double
sqrt 函数的通用数学库?
没有必要long double
;平方根可以用double
(如果它是 IEEE-754 64 位二进制)来计算。将 64 位整数转换为的舍入误差double
与此问题几乎无关。
舍入误差最多为 2 53的一部分。这会导致 2 54中最多一个部分的平方根出现误差。由于将数学结果四舍五入到格式,sqrt
它本身的舍入误差小于 2 53的一部分。double
这些错误的总和很小;64 位整数(四舍五入为 53 位)的最大可能平方根为 2 32 ,因此 2 54中三个部分的误差小于 0.00000072。
对于一个uint64_t x
,考虑sqrt(x)
。我们知道这个值在 的精确平方根的 .00000072 范围内x
,但我们不知道它的方向。如果我们将其调整为sqrt(x) - 0x1p-20
,那么我们知道我们的值小于但非常接近 的平方根x
。
然后此代码计算 的平方根x
,截断为整数,前提是操作符合 IEEE 754:
uint64_t y = sqrt(x) - 0x1p-20;
if (2*y < x - y*y)
++y;
(2*y < x - y*y
等效于(y+1)*(y+1) <= x
除了它避免包装 64 位整数 if y+1
is 2 32。)
Function sqrtl()
,取 a long double
,是 C99 的一部分。
请注意,您的编译平台不必实现long double
为 80 位扩展精度。它只需要double
像double
. GCC 和 Clang 确实在 Intel 处理器上编译long double
为 80 位扩展精度。
如果您只想计算整数的 sqrt,使用分治法应该在最多 32 次迭代中找到结果:
uint64_t mysqrt (uint64_t a)
{
uint64_t min=0;
//uint64_t max=1<<32;
uint64_t max=((uint64_t) 1) << 32; //chux' bugfix
while(1)
{
if (max <= 1 + min)
return min;
uint64_t sqt = min + (max - min)/2;
uint64_t sq = sqt*sqt;
if (sq == a)
return sqt;
if (sq > a)
max = sqt;
else
min = sqt;
}
调试留给读者作为练习。
是的,标准库有sqrtl()
(从 C99 开始)。
在这里,我们收集了几个观察结果以得出解决方案:
x
是无符号整数类型,则 tnenx >> 1 == x / 2
和x << 1 == x * 2
. sqrt(x)
数学上等价于exp(log(x)/2.0)
。 IntExp2( IntLog2(x) / 2) "==" IntSqrtDn(x)
数和以 2 为底的指数,我们可以获得一个公平的估计: "="
IntExp2( IntLog2(x) / 2 + 1) "==" IntSqrtUp(x)
,我们将获得整数平方根的“以上”近似值。 x
。 http://en.wikipedia.org/wiki/Integer_square_root
最终的算法需要一些数学证明以确保它始终正常工作,但我现在不会这样做......我将向您展示最终的程序,而不是:
#include <stdio.h> /* For printf()... */
#include <stdint.h> /* For uintmax_t... */
#include <math.h> /* For sqrt() .... */
int IntLog2(uintmax_t n) {
if (n == 0) return -1; /* Error */
int L;
for (L = 0; n >>= 1; L++)
;
return L; /* It takes < 64 steps for long long */
}
uintmax_t IntExp2(int n) {
if (n < 0)
return 0; /* Error */
uintmax_t E;
for (E = 1; n-- > 0; E <<= 1)
;
return E; /* It takes < 64 steps for long long */
}
uintmax_t IntSqrtDn(uintmax_t n) { return IntExp2(IntLog2(n) / 2); }
uintmax_t IntSqrtUp(uintmax_t n) { return IntExp2(IntLog2(n) / 2 + 1); }
int main(void) {
uintmax_t N = 947612934; /* Try here your number! */
uintmax_t sqrtn = IntSqrtDn(N), /* 1st approx. to sqrt(N) by below */
sqrtn0 = IntSqrtUp(N); /* 1st approx. to sqrt(N) by above */
/* The following means while( abs(sqrt-sqrt0) > 1) { stuff... } */
/* However, we take care of subtractions on unsigned arithmetic, just in case... */
while ( (sqrtn > sqrtn0 + 1) || (sqrtn0 > sqrtn+1) )
sqrtn0 = sqrtn, sqrtn = (sqrtn0 + N/sqrtn0) / 2; /* Newton iteration */
printf("N==%llu, sqrt(N)==%g, IntSqrtDn(N)==%llu, IntSqrtUp(N)==%llu, sqrtn==%llu, sqrtn*sqrtn==%llu\n\n",
N, sqrt(N), IntSqrtDn(N), IntSqrtUp(N), sqrtn, sqrtn*sqrtn);
return 0;
}
最后存储的值sqrtn
是 的整数平方根N
。
程序的最后一行只显示了所有的值,以防万一。
所以,你可以尝试不同的值,N
看看会发生什么。
如果我们在 while 循环中添加一个计数器,我们将看到不会超过几次迭代。
备注:abs(sqrtn-sqrtn0)<=1
在整数设置下工作时,需要验证是否始终满足条件。如果不是,我们将不得不修复算法。
备注 2:在初始化语句中,观察sqrtn0 == sqrtn * 2 == sqrtn << 1
. 这避免了我们的一些计算。
// sqrt_i64 returns the integer square root of v.
int64_t sqrt_i64(int64_t v) {
uint64_t q = 0, b = 1, r = v;
for( b <<= 62; b > 0 && b > r; b >>= 2);
while( b > 0 ) {
uint64_t t = q + b;
q >>= 1;
if( r >= t ) {
r -= t;
q += b;
}
b >>= 2;
}
return q;
}
可以通过使用clz
机器代码指令来优化 for 循环。