按位整数立方根算法答案

【问题标题】：Bitwise integer cube root algorithm按位整数立方根算法
【发布时间】：2011-07-10 13:17:48
【问题描述】：

以下是计算整数平方根的简单方法：

int isqrt(int num)
{
    int root=0;
    int b = 0x8000;
    int a=0, c=0;

    while (b) {
        c = a|b;

        if (c*c <= num)
            a |= b;   

        b >>= 1;
    }       
}

巧妙地（感谢Wikipedia），可以这样优化：

int sqrt(short num)
{
    int op = num;
    int res = 0;
    int one = 1 << 30;

    while (one > op)
        one >>= 2;

    while (one != 0) {
        if (op >= res + one) {
            op -= res + one;
            res = (res >> 1) + one;
        }
        else
          res >>= 1;
        one >>= 2;
    }
    return res;
}

我的问题：可以为整数立方根编写类似优化的算法吗？（这是在一个不喜欢做乘法的小型微控制器上运行的）

【问题讨论】：

我不知道，但是您可以通过添加 7、19、37、61 等来生成连续整数的三次方，您可以通过添加 12、18、24、30 来获得这些数字， 36 等。它不是特别聪明或快速，但考虑到 2^32 的整数立方根仍然只有 1625，它不应该进行那么多迭代（所有迭代都包括几个加法和比较，没有 mults）。编辑：所以事实证明有一种方法。很高兴知道！
是的，算法可以扩展到立方根，即使没有乘法。请参阅此代码：hackersdelight.org/HDcode/icbrt.c.txt 并考虑购买他的代码来自的书 Hackers Delight。如果您每年必须经常解决此类问题，那么您绝对应该阅读它！

标签： algorithm math optimization

【解决方案1】：

根据this SO question 和标记的答案，从 Hacker's Delight 书中你可以找到这个实现：

int icbrt2(unsigned x) {
   int s;
   unsigned y, b, y2;

   y2 = 0;
   y = 0;
   for (s = 30; s >= 0; s = s - 3) {
      y2 = 4*y2;
      y = 2*y;
      b = (3*(y2 + y) + 1) << s;
      if (x >= b) {
         x = x - b;
         y2 = y2 + 2*y + 1;
         y = y + 1;
      }
   }
   return y;
}

【讨论】：

【解决方案2】：

正如其他人所提到的，这是 Hacker's Delight 代码的（极端）C# 优化版本。

供参考（在我的电脑上）：Math.Sqrt 大约需要 35 ns，cbrt

使用小数的乘法，但很容易用移位替换它们补充道。例如最大的乘法（最后一行）： "12 * z" ==> "(z

很难判断代码的大小是否适合小型微控制器。

第一步：二分查找，“if”语句，大值（ >= 1u

第二步：跳入展开的循环，即“标签”。

private static uint cbrt(uint x)
{
    uint y = 2, z = 4, b = 0;
    if (x < 1u << 24)
        if (x < 1u << 12)
            if (x < 1u << 06)
                if (x < 1u << 03)
                    return x == 0u ? 0u : 1u;
                else
                    return x < 27u ? 2u : 3u;
            else
                if (x < 1u << 09) goto L8; else goto L7;
        else
            if (x < 1u << 18)
                if (x < 1u << 15) goto L6; else goto L5;
            else
                if (x < 1u << 21) goto L4; else goto L3;
    else
        if (x >= 1u << 30) goto L0;
        else
            if (x < 1u << 27) goto L2; else goto L1;

L0: x -= 1u << 30; if (x >= 19u << 27)
    { x -= 19u << 27; z = 9; y = 3; } goto M0;
L1: x -= 1u << 27; if (x >= 19u << 24)
    { x -= 19u << 24; z = 9; y = 3; } goto M1;
L2: x -= 1u << 24; if (x >= 19u << 21)
    { x -= 19u << 21; z = 9; y = 3; } goto M2;
L3: x -= 1u << 21; if (x >= 19u << 18)
    { x -= 19u << 18; z = 9; y = 3; } goto M3;
L4: x -= 1u << 18; if (x >= 19u << 15)
    { x -= 19u << 15; z = 9; y = 3; } goto M4;
L5: x -= 1u << 15; if (x >= 19u << 12)
    { x -= 19u << 12; z = 9; y = 3; } goto M5;
L6: x -= 1u << 12; if (x >= 19u << 09)
    { x -= 19u << 09; z = 9; y = 3; } goto M6;
L7: x -= 1u << 09; if (x >= 19u << 06)
    { x -= 19u << 06; z = 9; y = 3; } goto M7;
L8: x -= 1u << 06; if (x >= 19u << 03)
    { x -= 19u << 03; z = 9; y = 3; } goto M8;

M0: y *= 2; z *= 4; b = 3 * y + 3 * z + 1 << 24; 
    if (x >= b) { x -= b; z += 2 * y + 1; y += 1; }
M1: y *= 2; z *= 4; b = 3 * y + 3 * z + 1 << 21; 
    if (x >= b) { x -= b; z += 2 * y + 1; y += 1; }
M2: y *= 2; z *= 4; b = 3 * y + 3 * z + 1 << 18;
    if (x >= b) { x -= b; z += 2 * y + 1; y += 1; }
M3: y *= 2; z *= 4; b = 3 * y + 3 * z + 1 << 15;
    if (x >= b) { x -= b; z += 2 * y + 1; y += 1; }
M4: y *= 2; z *= 4; b = 3 * y + 3 * z + 1 << 12; 
    if (x >= b) { x -= b; z += 2 * y + 1; y += 1; }
M5: y *= 2; z *= 4; b = 3 * y + 3 * z + 1 << 09; 
    if (x >= b) { x -= b; z += 2 * y + 1; y += 1; }
M6: y *= 2; z *= 4; b = 3 * y + 3 * z + 1 << 06; 
    if (x >= b) { x -= b; z += 2 * y + 1; y += 1; }
M7: y *= 2; z *= 4; b = 3 * y + 3 * z + 1 << 03; 
    if (x >= b) { x -= b; z += 2 * y + 1; y += 1; }
M8: y *= 2; return x <= 3 * y + 12 * z ? y : y + 1;
}

【讨论】：