整数立方根
我正在寻找64位（无符号）立方根的快速代码。 （我正在使用C并使用gcc进行编译，但我想大多数所需的工作都是语言和编译器无关的。）我将通过ulong表示一个64位的unisgned整数。
给定输入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 = 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; }
“Hacker’s Delight”这本书有针对这个和许多其他问题的算法。 代码在这里在线。 编辑 ：该代码无法正常使用64位整数，书中有关如何为64位修复它的说明有些令人困惑。 这里有一个正确的64位实现（包括测试用例）。
我怀疑你的squareroot函数是否“正确” – 对于参数而言应该是ulong a ，而不是n :)（但是使用cbrt而不是sqrt方法相同，尽管不是所有的C数学库都有立方根函数）。
您可以尝试Newton的步骤来修复舍入错误：
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;
或者其他一些。
（我承认我很懒;正确的方法是仔细检查以确定哪些（如果有的话）检查和增量实际上是必要的......）
如果pow太昂贵，你可以使用count-leading-zero指令来得到结果的近似值，然后使用查找表，然后使用一些牛顿步骤来完成它。
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)的良好近似值。 使用一些牛顿步骤（参见暴风雨的答案）来完成它。
我将研究如何手动完成 ，然后将其转换为计算机算法，在基础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的值，并使用移位而不是乘法。
// 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; }