20

我正在寻找 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;
}
4

6 回答 6

12

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

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

于 2010-12-02T05:09:24.590 回答
4

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

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

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

我使用的初始估计是基于对值 ( x )中有效位数的“四舍五入” 。给定 ( x ) 中的 ( b ) 有效位，我们可以说：. 我在没有证据的情况下声明（尽管它应该相对容易证明）：2^(b - 1) <= x < 2^b2^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)); */
}

A crbt call probably isn't all that useful - unlike the sqrt call which can be efficiently implemented on modern hardware. That said, I've seen gains for sets of values under 2^53 (exactly represented in IEEE-754 doubles), which surprised me.

The only downside is the division by: (r * r) - this can be slow, as the latency of integer division continues to fall behind other advances in ALUs. The division by a constant: (3) is handled by reciprocal methods on any modern optimising compiler.

It's interesting that Intel's 'Icelake' microarchitecture will significantly improve integer division - an operation that seems to have been neglected for a long time. I simply won't trust the 'Hacker's Delight' answer until I can find a sound theoretical basis for it. And then I have to work out which variant is the 'correct' answer.

于 2019-06-24T13:48:25.570 回答
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;

或一些这样的。

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

于 2010-12-02T18:31:47.863 回答
3

如果pow太昂贵，您可以使用 count-leading-zeros 指令来获得结果的近似值，然后使用查找表，然后使用一些 Newton 步骤来完成它。

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];

给定band top8，您可以使用查找表（在我的代码中，8K 条目）找到cuberoot(n). 使用一些牛顿步骤（见comingstorm's answer）来完成它。

于 2010-12-02T21:38:42.210 回答
1 
// 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;
}
于 2013-09-24T15:11:21.817 回答
0

我会研究如何手工完成，然后将其转化为计算机算法，以 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迭代之间的值以及使用移位而不是乘法来优化 的计算。

于 2010-12-02T06:25:59.560 回答