6

我正在查看另一个问题（此处），有人正在寻找一种方法来获取 x86 程序集中的 64 位整数的平方根。

事实证明这很简单。解决方案是转换为浮点数，计算 sqrt，然后再转换回来。

我需要在 C 中做一些非常相似的事情，但是当我研究等价物时，我有点卡住了。我只能找到一个接受双打的 sqrt 函数。双精度数不具备存储大型 64 位整数而不引入显着舍入误差的精度。

是否有一个我可以使用的具有long doublesqrt 函数的通用数学库？

4

6 回答 6

16

没有必要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+1is 2 32。)

于 2013-08-29T02:06:43.507 回答
5

Function sqrtl()，取 a long double，是 C99 的一部分。

请注意，您的编译平台不必实现long double为 80 位扩展精度。它只需要double像double. GCC 和 Clang 确实在 Intel 处理器上编译long double为 80 位扩展精度。

于 2013-08-28T22:46:33.297 回答
2

如果您只想计算整数的 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;
    }

调试留给读者作为练习。

于 2013-08-29T01:09:38.763 回答
2

是的，标准库有sqrtl()（从 C99 开始）。

于 2013-08-28T22:45:52.983 回答
2

在这里，我们收集了几个观察结果以得出解决方案：

  1. 在标准 C >= 1999 中，保证非净整数具有位表示，正如任何以 2 为基数的数字所期望的那样。
    ----> 因此，我们可以相信这种类型的数字的位操作。
  2. 如果x是无符号整数类型，则 tnenx >> 1 == x / 2和x << 1 == x * 2.
    （！）但是：很可能位运算应该比它们的算术对应物更快地完成。
  3. sqrt(x)数学上等价于exp(log(x)/2.0)。
  4. 如果我们考虑整数的截断对IntExp2( IntLog2(x) / 2) "==" IntSqrtDn(x)数和以 2 为底的指数，我们可以获得一个公平的估计： "="
  5. 如果我们写IntExp2( IntLog2(x) / 2 + 1) "==" IntSqrtUp(x)，我们将获得整数平方根的“以上”近似值。
  6. 在 (4.) 和 (5.) 中获得的近似值有点粗略（它们将 sqrt(x) 的真实值包含在 2 的两个连续幂之间），但它们可能是任何搜索算法的一个很好的起点为广场的屋顶x。
  7. 如果我们对真实解有一个很好的初步近似，那么求平方根 的牛顿算法可能对整数很有效。

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. 这避免了我们的一些计算。

于 2013-08-30T00:17:12.660 回答
0 
// 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 循环。

于 2021-05-15T14:41:35.320 回答