获得 sqrt (n)整数部分的最快方法？

小开

虽然我怀疑你可以通过搜索“快速整数平方根”找到大量的选项，这里有一些潜在的新想法可能会很好地工作(每个独立的，或者也许你可以组合它们) :

创建一个 static const数组，其中包含您希望支持的域中的所有完美正方形，并对其执行快速无分支二进制搜索。数组中的结果索引是平方根。
将数字转换为浮点数，并将其分解为尾数和指数。将指数减半，然后将尾数乘以某个神奇的因子(你的工作就是找到它)。这应该能够给出一个非常接近的近似值。包括一个最后的步骤来调整它，如果它不是准确的(或使用它作为一个起点为二进制搜索以上)。

小开

我会试试平方根倒数速算法的把戏。

这是一种在没有任何分支的情况下获得非常好的 1/sqrt(n)近似值的方法，基于一些不可移植的位操作(特别是在32位和64位平台之间)。

一旦你得到了它，你只需要反转结果，取出整数部分。

当然，可能还有更快的技巧，因为这个有点绕圈子。

让我们开始吧！

首先是一个小帮手:

// benchmark.h #include <sys/time.h> template <typename Func> double benchmark(Func f, size_t iterations) { f(); timeval a, b; gettimeofday(&a, 0); for (; iterations --> 0;) { f(); } gettimeofday(&b, 0); return (b.tv_sec * (unsigned int)1e6 + b.tv_usec) - (a.tv_sec * (unsigned int)1e6 + a.tv_usec); }

然后是主体:

#include <iostream> #include <cmath> #include "benchmark.h" class Sqrt { public: Sqrt(int n): _number(n) {} int operator()() const { double d = _number; return static_cast<int>(std::sqrt(d) + 0.5); } private: int _number; }; // http://www.codecodex.com/wiki/Calculate_an_integer_square_root class IntSqrt { public: IntSqrt(int n): _number(n) {} int operator()() const { int remainder = _number; if (remainder < 0) { return 0; } int place = 1 <<(sizeof(int)*8 -2); while (place > remainder) { place /= 4; } int root = 0; while (place) { if (remainder >= root + place) { remainder -= root + place; root += place*2; } root /= 2; place /= 4; } return root; } private: int _number; }; // http://en.wikipedia.org/wiki/Fast_inverse_square_root class FastSqrt { public: FastSqrt(int n): _number(n) {} int operator()() const { float number = _number; float x2 = number * 0.5F; float y = number; long i = *(long*)&y; //i = (long)0x5fe6ec85e7de30da - (i >> 1); i = 0x5f3759df - (i >> 1); y = *(float*)&i; y = y * (1.5F - (x2*y*y)); y = y * (1.5F - (x2*y*y)); // let's be precise return static_cast<int>(1/y + 0.5f); } private: int _number; }; int main(int argc, char* argv[]) { if (argc != 3) { std::cerr << "Usage: %prog integer iterations\n"; return 1; } int n = atoi(argv[1]); int it = atoi(argv[2]); assert(Sqrt(n)() == IntSqrt(n)() && Sqrt(n)() == FastSqrt(n)() && "Different Roots!"); std::cout << "sqrt(" << n << ") = " << Sqrt(n)() << "\n"; double time = benchmark(Sqrt(n), it); double intTime = benchmark(IntSqrt(n), it); double fastTime = benchmark(FastSqrt(n), it); std::cout << "Number iterations: " << it << "\n" "Sqrt computation : " << time << "\n" "Int computation : " << intTime << "\n" "Fast computation : " << fastTime << "\n"; return 0; }

结果是:

sqrt(82) = 9 Number iterations: 4096 Sqrt computation : 56 Int computation : 217 Fast computation : 119 // Note had to tweak the program here as Int here returns -1 :/ sqrt(2147483647) = 46341 // real answer sqrt(2 147 483 647) = 46 340.95 Number iterations: 4096 Sqrt computation : 57 Int computation : 313 Fast computation : 119

正如预期的那样，快点计算比内景计算执行得好得多。

哦，顺便说一句，sqrt更快:)

小开

如果您需要计算平方根的性能，我想您会计算很多平方根。那为什么不缓存答案呢？我不知道在你的例子中 N 的范围，也不知道你是否会计算同一个整数的平方根的许多倍，但是如果是，那么你可以在每次调用你的方法时缓存结果(在一个数组中如果不是太大的话将是最有效的)。

小开

编辑: 这个答案是愚蠢的-使用 (int) sqrt(i)

在使用 适当的设置(-march=native -m64 -O3)分析之后，上面的很多速度更快。

好吧，这个问题有点老套，但是“最快”的答案还没有给出。最快的(我认为)是二进制平方根算法，在这篇 Embedded.com 的文章中有完整的解释。

基本上可以归结为:

unsigned short isqrt(unsigned long a) { unsigned long rem = 0; int root = 0; int i; for (i = 0; i < 16; i++) { root <<= 1; rem <<= 2; rem += a >> 30; a <<= 2; if (root < rem) { root++; rem -= root; root++; } } return (unsigned short) (root >> 1); }

在我的机器(Q6600，Ubuntu 10.10)上，我通过取数字1-100000000的平方根进行分析。使用 iqsrt(i)耗时2750毫秒，使用 (unsigned short) sqrt((float) i)耗时3600毫秒。这是使用 g++ -O3完成的。使用 -ffast-math编译选项的时间分别为2100ms 和3100ms。请注意，这是没有使用甚至一个单一的汇编程序行，所以它可能仍然快得多。

上面的代码同时适用于 C 和 C + + ，并且对 Java 的语法也有很小的改变。

在有限的范围内，二进制搜索效果更好。在我的机器上，这把上面的版本从水中吹出4倍。遗憾的是，它的射程非常有限:

#include <stdint.h> const uint16_t squares[] = { 0, 1, 4, 9, 16, 25, 36, 49, 64, 81, 100, 121, 144, 169, 196, 225, 256, 289, 324, 361, 400, 441, 484, 529, 576, 625, 676, 729, 784, 841, 900, 961, 1024, 1089, 1156, 1225, 1296, 1369, 1444, 1521, 1600, 1681, 1764, 1849, 1936, 2025, 2116, 2209, 2304, 2401, 2500, 2601, 2704, 2809, 2916, 3025, 3136, 3249, 3364, 3481, 3600, 3721, 3844, 3969, 4096, 4225, 4356, 4489, 4624, 4761, 4900, 5041, 5184, 5329, 5476, 5625, 5776, 5929, 6084, 6241, 6400, 6561, 6724, 6889, 7056, 7225, 7396, 7569, 7744, 7921, 8100, 8281, 8464, 8649, 8836, 9025, 9216, 9409, 9604, 9801, 10000, 10201, 10404, 10609, 10816, 11025, 11236, 11449, 11664, 11881, 12100, 12321, 12544, 12769, 12996, 13225, 13456, 13689, 13924, 14161, 14400, 14641, 14884, 15129, 15376, 15625, 15876, 16129, 16384, 16641, 16900, 17161, 17424, 17689, 17956, 18225, 18496, 18769, 19044, 19321, 19600, 19881, 20164, 20449, 20736, 21025, 21316, 21609, 21904, 22201, 22500, 22801, 23104, 23409, 23716, 24025, 24336, 24649, 24964, 25281, 25600, 25921, 26244, 26569, 26896, 27225, 27556, 27889, 28224, 28561, 28900, 29241, 29584, 29929, 30276, 30625, 30976, 31329, 31684, 32041, 32400, 32761, 33124, 33489, 33856, 34225, 34596, 34969, 35344, 35721, 36100, 36481, 36864, 37249, 37636, 38025, 38416, 38809, 39204, 39601, 40000, 40401, 40804, 41209, 41616, 42025, 42436, 42849, 43264, 43681, 44100, 44521, 44944, 45369, 45796, 46225, 46656, 47089, 47524, 47961, 48400, 48841, 49284, 49729, 50176, 50625, 51076, 51529, 51984, 52441, 52900, 53361, 53824, 54289, 54756, 55225, 55696, 56169, 56644, 57121, 57600, 58081, 58564, 59049, 59536, 60025, 60516, 61009, 61504, 62001, 62500, 63001, 63504, 64009, 64516, 65025 }; inline int isqrt(uint16_t x) { const uint16_t *p = squares; if (p[128] <= x) p += 128; if (p[ 64] <= x) p += 64; if (p[ 32] <= x) p += 32; if (p[ 16] <= x) p += 16; if (p[ 8] <= x) p += 8; if (p[ 4] <= x) p += 4; if (p[ 2] <= x) p += 2; if (p[ 1] <= x) p += 1; return p - squares; }

一个32位版本可以在这里下载: https://gist.github.com/3481770

小开

要执行整数 sqrt，可以使用 Newton 方法的这种特殊化:

Def isqrt(N): a = 1 b = N while |a-b| > 1 b = N / a a = (a + b) / 2 return a

基本上对于任何 x，sqrt 都在(x... N/x)的范围内，所以我们只需要对新猜测的每个循环中的间隔进行平分。有点像二进制搜索，但收敛得更快。

这个收敛在 O (loglog (N))中非常快。它也根本不使用浮点数，而且对于任意精度的整数也能很好地工作。

小开

为什么没有人建议最快的方法？

如果:

数字的范围是有限的

内存消耗并不重要

应用程序启动时间并不重要

然后用 sqrt(x)(不需要使用函数 sqrt())创建填充 int[MAX_X](在启动时)。

所有这些条件都很适合我的计划。特别是，int[10000000]阵列将消耗 40MB。

你有什么想法？

小开

在许多情况下，甚至不需要精确的整数 sqrt 值，只需要对它进行很好的近似即可。(例如，它经常发生在 DSP 优化，当32位信号应该压缩到16位，或16位到8位，而不失去很多精度在零附近)。

我发现了一个有用的公式:

k = ceil(MSB(n)/2); - MSB(n) is the most significant bit of "n"

sqrt(n) ~= 2^(k-2)+(2^(k-1))*n/(2^(2*k))); - all multiplications and divisions here are very DSP-friendly, as they are only 2^k.

这个方程产生光滑曲线(n，sqrt (n)) ，它的值与实数 sqrt (n)相差不大，因此在近似精度足够的情况下可以使用。

小开

在我的计算机上，使用 gcc 和-fast-數学，将一个32位整数转换为 float 并使用 sqrtf 需要1.2 s/10 ^ 9次运算(如果不使用-fast-數学，则需要3.54 s)。

下面的算法使用了0.87 s/10 ^ 9，但牺牲了一些精度: 误差可能高达 -7或 + 1，尽管 RMS 误差只有0.79:

uint16_t SQRTTAB[65536]; inline uint16_t approxsqrt(uint32_t x) { const uint32_t m1 = 0xff000000; const uint32_t m2 = 0x00ff0000; if (x&m1) { return SQRTTAB[x>>16]; } else if (x&m2) { return SQRTTAB[x>>8]>>4; } else { return SQRTTAB[x]>>8; } }

表格的构造使用:

void maketable() { for (int x=0; x<65536; x++) { double v = x/65535.0; v = sqrt(v); int y = int(v*65535.0+0.999); SQRTTAB[x] = y; } }

我发现，使用 if 语句进一步细化二分法确实提高了准确性，但它也减缓了 sqrtf 的速度，至少在-fast-數学中是这样。

小开

如果你不介意一个近似值，我拼凑的这个整数 sqrt 函数怎么样。

int sqrti(int x) { union { float f; int x; } v; // convert to float v.f = (float)x; // fast aprox sqrt // assumes float is in IEEE 754 single precision format // assumes int is 32 bits // b = exponent bias // m = number of mantissa bits v.x -= 1 << 23; // subtract 2^m v.x >>= 1; // divide by 2 v.x += 1 << 29; // add ((b + 1) / 2) * 2^m // convert to int return (int)v.f; }

它使用本维基百科文章中描述的算法。在我的机器上，它的速度几乎是 sqrt 的两倍:)

小开

这是如此之短，它99% 内联:

static inline int sqrtn(int num) { int i = 0; __asm__ ( "pxor %%xmm0, %%xmm0\n\t" // clean xmm0 for cvtsi2ss "cvtsi2ss %1, %%xmm0\n\t" // convert num to float, put it to xmm0 "sqrtss %%xmm0, %%xmm0\n\t" // square root xmm0 "cvttss2si %%xmm0, %0" // float to int :"=r"(i):"r"(num):"%xmm0"); // i: result, num: input, xmm0: scratch register return i; }

为什么要清理 xmm0? cvtsi2ss的文档

目标操作数是一个 XMM 寄存器。结果存储在目标操作数的低双字中，上面三个双字保持不变。

GCC 内部版本(仅在 GCC 上运行) :

#include <xmmintrin.h> int sqrtn2(int num) { register __v4sf xmm0 = {0, 0, 0, 0}; xmm0 = __builtin_ia32_cvtsi2ss(xmm0, num); xmm0 = __builtin_ia32_sqrtss(xmm0); return __builtin_ia32_cvttss2si(xmm0); }

Intel 内部版本(在 GCC，Clang，ICC 上测试) :

#include <xmmintrin.h> int sqrtn2(int num) { register __m128 xmm0 = _mm_setzero_ps(); xmm0 = _mm_cvt_si2ss(xmm0, num); xmm0 = _mm_sqrt_ss(xmm0); return _mm_cvtt_ss2si(xmm0); }

它们都需要 SSE 1(甚至不需要 SSE 2)。

注意 : 这正是 GCC 使用 -Ofast计算 (int) sqrt((float) num)的方法。如果你想要更大的 i有更高的精度，那么我们可以计算 (int) sqrt((double) num)(正如 Gumby The Green 在评论中指出的那样) :

static inline int sqrtn(int num) { int i = 0; __asm__ ( "pxor %%xmm0, %%xmm0\n\t" "cvtsi2sd %1, %%xmm0\n\t" "sqrtsd %%xmm0, %%xmm0\n\t" "cvttsd2si %%xmm0, %0" :"=r"(i):"r"(num):"%xmm0"); return i; }

或者

#include <xmmintrin.h> int sqrtn2(int num) { register __v2df xmm0 = {0, 0}; xmm0 = __builtin_ia32_cvtsi2sd(xmm0, num); xmm0 = __builtin_ia32_sqrtsd(xmm0); return __builtin_ia32_cvttsd2si(xmm0); }

小开

下面的解决方案计算整数部分，准确地表示 floor(sqrt(x))，没有舍入错误。

其他方法的问题

使用 float或 double既不便携也不够精确

@ orlp 的 isqrt得到了和 isqrt(100) = 15一样疯狂的结果

基于大型查找表的方法在32位以上不实用

使用 快速反向 sqrt 快速反向 sqrt是非常不精确的，你最好使用 sqrtf

牛顿的方法需要昂贵的整数除法和良好的初始猜测

我的方法

我的是基于维基百科提出的比特猜测方法的。不幸的是，维基百科上提供的伪代码有一些错误，所以我不得不做一些调整:

// C++20 also provides std::bit_width in its <bit> header unsigned char bit_width(unsigned long long x) { return x == 0 ? 1 : 64 - __builtin_clzll(x); } template <typename Int, std::enable_if_t<std::is_unsigned<Int, int = 0>> Int sqrt(const Int n) { unsigned char shift = bit_width(n); shift += shift & 1; // round up to next multiple of 2 Int result = 0; do { shift -= 2; result <<= 1; // make space for the next guessed bit result |= 1; // guess that the next bit is 1 result ^= result * result > (n >> shift); // revert if guess too high } while (shift != 0); return result; }

bit_width可以在恒定的时间内计算，循环最多迭代 ceil(bit_width / 2)次。因此，即使对于64位整数，这最多也只是基本算术和按位运算的32次迭代。

编译输出只有大约20个指令。

表演

我已经基准测试我的方法对 float基地的通过生成输入一致。请注意，在现实世界中，大多数输入将比 std::numeric_limits<...>::max()更接近于零。

对于 uint32_t，它的性能比使用 std::sqrt(float)差

对于 uint64_t，它的性能比使用 std::sqrt(double)差

准确性

与使用浮点数学的方法不同，这种方法总是非常精确。

使用 sqrtf可以在[2²⁸,2³²)范围内提供不正确的舍入。例如，sqrtf(0xffffffff) = 65536，当平方根实际上是 65535.99999时。

在[2⁶⁰,2⁶⁴)范围内，双精度不能始终如一地工作。例如，当平方根实际上是 2147483647.999999时，就是 sqrt(0x3fff...) = 2147483648。

唯一覆盖所有64位整数的是 x86扩展精度 long double，因为它可以容纳整个64位整数。

结论

正如我所说，这是唯一正确处理所有输入、避免整数除法和不需要查找表的解决方案。总之，如果您需要一个与精度无关且不需要巨大查找表的方法，那么这是您唯一的选择。在 constexpr环境中，它可能特别有用，因为在 constexpr环境中，性能并不重要，而且获得100% 准确的结果可能更为重要。

采用牛顿法的替代方法

牛顿的方法可以相当快，当开始一个很好的猜测。对于我们的猜测，我们将四舍五入到下一个2的幂，并计算平方根在恒定的时间。对于任意数字2^X，我们可以使用2^X/2得到平方根。

template <typename Int, std::enable_if_t<std::is_unsigned_v<Int>, int> = 0> Int sqrt_guess(const Int n) { Int log2floor = bit_width(n) - 1; // sqrt(x) is equivalent to pow(2, x / 2 = x >> 1) // pow(2, x) is equivalent to 1 << x return 1 << (log2floor >> 1); }

注意，这并不是2^X/2，因为我们在右移时失去了一些精度，而是2^楼面(x/2)。还要注意的是，sqrt_guess(0) = 1实际上是在第一次迭代中避免被零除所必需的:

template <typename Int, std::enable_if_t<std::is_unsigned_v<Int>, int> = 0> Int sqrt_newton(const Int n) { Int a = sqrt_guess(n); Int b = n; // compute unsigned difference while (std::max(a, b) - std::min(a, b) > 1) { b = n / a; a = (a + b) / 2; } // a is now either floor(sqrt(n)) or ceil(sqrt(n)) // we decrement in the latter case // this is overflow-safe as long as we start with a lower bound guess return a - (a * a > n); }

这种替代方法的性能大致相当于第一个方案，但通常比第一个方案快几个百分点。然而，它在很大程度上依赖于有效的硬件划分，结果可能差异很大。

sqrt_guess的使用产生了巨大的差异。它大约比使用 1作为初始猜测快5倍。

小开

或者只是做一个二进制搜索，不能写一个简单的版本 imo:

uint16_t sqrti(uint32_t num) { uint16_t ret = 0; for(int32_t i = 15; i >= 0; i--) { uint16_t temp = ret | (1 << i); if(temp * temp <= num) { ret = temp; } } return ret; }

获得 sqrt (n)整数部分的最快方法？

编辑: 这个答案是愚蠢的-使用 (int) sqrt(i)

其他方法的问题

我的方法

表演

准确性

结论

采用牛顿法的替代方法

编辑: 这个答案是愚蠢的-使用 `(int) sqrt(i)`