整数立方根

Question

我正在寻找 64 位（无符号）立方根的快速代码。 （我使用 C 并使用 gcc 进行编译，但我认为所需的大部分工作将与语言和编译器无关。）我将用 ulong 表示一个 64 位无符号整数。

给定一个输入 n，我要求（整数）返回值 r 是这样的

r * r * r <= n && n < (r + 1) * (r + 1) * (r + 1)

也就是说，我想要 n 的立方根，向下取整。基本代码如

return (ulong)pow(n, 1.0/3);

不正确，因为向范围末尾舍入。像

这样的简单代码

ulong
cuberoot(ulong n)
{
    ulong ret = pow(n + 0.5, 1.0/3);
    if (n < 100000000000001ULL)
        return ret;
    if (n >= 18446724184312856125ULL)
        return 2642245ULL;
    if (ret * ret * ret > n) {
        ret--;
        while (ret * ret * ret > n)
            ret--;
        return ret;
    }
    while ((ret + 1) * (ret + 1) * (ret + 1) <= n)
        ret++;
    return ret;
}

给出正确的结果，但比需要的慢。

此代码用于数学库，它将从各种函数中多次调用。速度很重要，但您不能指望温暖的缓存（所以像 2,642,245 项二进制搜索这样的建议就出来了）。

为了比较，这里是正确计算整数平方根的代码。

ulong squareroot(ulong a) {
    ulong x = (ulong)sqrt((double)a);
    if (x > 0xFFFFFFFF || x*x > a)
        x--;
    return x;
}

Answer 1

“Hacker's Delight”一书有解决这个问题和许多其他问题的算法。代码在线here。 编辑：该代码不适用于 64 位整数，并且书中关于如何将其修复为 64 位的说明有些令人困惑。正确的 64 位实现（包括测试用例）在线here。

我怀疑您的 squareroot 函数是否“正确”工作 - 参数应该是 ulong a，而不是 n :)（但同样的方法可以使用 cbrt 而不是 sqrt，尽管并非所有 C 数学库都有立方根函数）。

Answer 2

我已经在Modern Computer Arithmetic (Brent and Zimmerman) 中调整了1.5.2（kth 根）中提出的算法。对于(k == 3) 的情况，并且给出了对初始猜测的“相对”准确的高估 - 该算法似乎优于上面的“Hacker's Delight”代码。

不仅如此，MCA 作为文本还提供了理论背景以及正确性和终止标准的证明。

假设我们可以产生一个“相对”好的初始高估，我还没有找到超过 (7) 次迭代的案例。 （这是否与具有 2^6 位的 64 位值有效相关？）无论哪种方式，它都是对 HacDel 代码中 (21) 次迭代的改进 - 尽管具有线性 O(b) 收敛性有一个明显更快的循环体。

我使用的初始估计是基于对值中有效位数 (x) 的“四舍五入”。给定 (x) 中的 (b) 个有效位，我们可以说：2^(b - 1) <= x < 2^b。我在没有证据的情况下声明（虽然它应该相对容易证明）：2^ceil(b / 3) > x^(1/3)

---

static inline uint32_t u64_cbrt (uint64_t x)
{
    uint64_t r0 = 1, r1;

    /* IEEE-754 cbrt *may* not be exact. */

    if (x == 0) /* cbrt(0) : */
        return (0);

    int b = (64) - __builtin_clzll(x);
    r0 <<= (b + 2) / 3; /* ceil(b / 3) */

    do /* quadratic convergence: */
    {
        r1 = r0;
        r0 = (2 * r1 + x / (r1 * r1)) / 3;
    }
    while (r0 < r1);

    return ((uint32_t) r1); /* floor(cbrt(x)); */
}

---

crbt 调用可能不是那么有用 - 不像 sqrt 调用可以在现代硬件上有效实现。也就是说，我看到2^53 下的值集（精确地表示在 IEEE-754 双精度值中）得到了提升，这让我感到惊讶。

唯一的缺点是除法：(r * r) - 这可能会很慢，因为整数除法的延迟继续落后于 ALU 的其他进步。除以常量：(3) 由任何现代优化编译器上的倒数方法处理。

有趣的是，英特尔的“Icelake”微架构将显着改进整数除法——这一操作似乎被忽视了很长时间。在找到合理的理论基础之前，我根本不会相信“黑客的喜悦”的答案。然后我必须找出哪个变体是“正确”的答案。

Answer 3

您可以尝试使用牛顿法来修正舍入错误：

ulong r = (ulong)pow(n, 1.0/3);
if(r==0) return r; /* avoid divide by 0 later on */
ulong r3 = r*r*r;
ulong slope = 3*r*r;

ulong r1 = r+1;
ulong r13 = r1*r1*r1;

/* making sure to handle unsigned arithmetic correctly */
if(n >= r13) r+= (n - r3)/slope;
if(n < r3)   r-= (r3 - n)/slope;

一个牛顿步骤应该就足够了，但你可能会遇到一个错误（或者可能更多？）错误。您可以使用最终检查和增量步骤来检查/修复这些问题，就像在您的 OQ 中一样：

while(r*r*r > n) --r;
while((r+1)*(r+1)*(r+1) <= n) ++r;

或类似的。

（我承认我很懒惰；正确的做法是仔细检查以确定哪些（如果有）检查和增量的东西实际上是必要的......）

Answer 4

如果pow 太昂贵，您可以使用前导零指令来获得结果的近似值，然后使用查找表，然后使用一些牛顿步骤来完成它。

int k = __builtin_clz(n); // counts # of leading zeros (often a single assembly insn)
int b = 64 - k;           // # of bits in n
int top8 = n >> (b - 8);  // top 8 bits of n (top bit is always 1)
int approx = table[b][top8 & 0x7f];

鉴于b 和top8，您可以使用查找表（在我的代码中，8K 条目）找到cuberoot(n) 的良好近似值。使用一些牛顿步骤（参见comingstorm的答案）来完成它。

Answer 5

// On my pc: Math.Sqrt 35 ns, cbrt64 <70ns, cbrt32 <25 ns, (cbrt12 < 10ns)

// cbrt64(ulong x) is a C# version of:
// http://www.hackersdelight.org/hdcodetxt/acbrt.c.txt     (acbrt1)

// cbrt32(uint x) is a C# version of:
// http://www.hackersdelight.org/hdcodetxt/icbrt.c.txt     (icbrt1)

// Union in C#:
// http://www.hanselman.com/blog/UnionsOrAnEquivalentInCSairamasTipOfTheDay.aspx

using System.Runtime.InteropServices;  
[StructLayout(LayoutKind.Explicit)]  
public struct fu_32   // float <==> uint
{
[FieldOffset(0)]
public float f;
[FieldOffset(0)]
public uint u;
}

private static uint cbrt64(ulong x)
{
    if (x >= 18446724184312856125) return 2642245;
    float fx = (float)x;
    fu_32 fu32 = new fu_32();
    fu32.f = fx;
    uint uy = fu32.u / 4;
    uy += uy / 4;
    uy += uy / 16;
    uy += uy / 256;
    uy += 0x2a5137a0;
    fu32.u = uy;
    float fy = fu32.f;
    fy = 0.33333333f * (fx / (fy * fy) + 2.0f * fy);
    int y0 = (int)                                      
        (0.33333333f * (fx / (fy * fy) + 2.0f * fy));    
    uint y1 = (uint)y0;                                 

    ulong y2, y3;
    if (y1 >= 2642245)
    {
        y1 = 2642245;
        y2 = 6981458640025;
        y3 = 18446724184312856125;
    }
    else
    {
        y2 = (ulong)y1 * y1;
        y3 = y2 * y1;
    }
    if (y3 > x)
    {
        y1 -= 1;
        y2 -= 2 * y1 + 1;
        y3 -= 3 * y2 + 3 * y1 + 1;
        while (y3 > x)
        {
            y1 -= 1;
            y2 -= 2 * y1 + 1;
            y3 -= 3 * y2 + 3 * y1 + 1;
        }
        return y1;
    }
    do
    {
        y3 += 3 * y2 + 3 * y1 + 1;
        y2 += 2 * y1 + 1;
        y1 += 1;
    }
    while (y3 <= x);
    return y1 - 1;
}

private static uint cbrt32(uint x)
{
    uint y = 0, z = 0, b = 0;
    int s = x < 1u << 24 ? x < 1u << 12 ? x < 1u << 06 ? x < 1u << 03 ? 00 : 03 :
                                                         x < 1u << 09 ? 06 : 09 :
                                          x < 1u << 18 ? x < 1u << 15 ? 12 : 15 :
                                                         x < 1u << 21 ? 18 : 21 :
                           x >= 1u << 30 ? 30 : x < 1u << 27 ? 24 : 27;
    do
    {
        y *= 2;
        z *= 4;
        b = 3 * y + 3 * z + 1 << s;
        if (x >= b)
        {
            x -= b;
            z += 2 * y + 1;
            y += 1;
        }
        s -= 3;
    }
    while (s >= 0);
    return y;
}

private static uint cbrt12(uint x) // x < ~255
{
    uint y = 0, a = 0, b = 1, c = 0;
    while (a < x)
    {
        y++;
        b += c;
        a += b;
        c += 6;
    }
    if (a != x) y--;
    return y;
} 

Answer 6

我会research how to do it by hand，然后将其转换为计算机算法，以 2 为底而不是 10 为底。

我们最终得到一个类似于（伪代码）的算法：

Find the largest n such that (1 << 3n) < input.
result = 1 << n.
For i in (n-1)..0:
    if ((result | 1 << i)**3) < input:
        result |= 1 << i.

我们可以优化(result | 1 << i)**3的计算，观察按位或等价于加法，重构为result**3 + 3 * i * result ** 2 + 3 * i ** 2 * result + i ** 3，在迭代之间缓存result**3和result**2的值，使用移位代替乘法.