获取 sqrt(n) 整数部分的最快方法？答案

【问题标题】：Fastest way to get the integer part of sqrt(n)?获取 sqrt(n) 整数部分的最快方法？
【发布时间】：2011-06-23 06:02:54
【问题描述】：

我们知道如果n 不是一个完美的正方形，那么sqrt(n) 就不是一个整数。由于我只需要整数部分，我觉得调用sqrt(n) 不会那么快，因为计算小数部分也需要时间。

所以我的问题是，

我们能否只得到sqrt(n)的整数部分而不计算sqrt(n)的实际值？该算法应该比sqrt(n)（定义在<math.h>或<cmath>）更快？

如果可能，您也可以在asm 块中编写代码。

【问题讨论】：

大多数 CPU 在硬件中执行 sqrt，因此仅计算整数部分不太可能加快速度。
这是一个更确定性算法的有趣链接：embedded-systems.com/98/9802fe2.htm
@Gabe: sqrt() C 库中的 sqrt() 不太可能在所有机器上直接实现为硬件 sqrt 指令，因为硬件可能无法处理 IEEE 754 要求的所有极端情况。如果你不关心，你可以使用内联 asm 或 gcc 的-ffast-math 直接获取硬件。
也许这个link可以帮助你。
您是否对您的应用程序进行了概要分析？您确定需要提高 sqrt(n) 速度吗？

标签： c++ c algorithm math performance

【解决方案1】：

虽然我怀疑您可以通过搜索“快速整数平方根”找到很多选项，但这里有一些可能很好用的潜在新想法（每个想法独立，或者您可以将它们组合起来）：

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

【讨论】：

如果您想在二分搜索的每一步都对“索引”进行平方，请成为我的客人。这将是slooooooow。这就是为什么我建议预先计算它们。请注意，我说的是static const。计算它是免费的，因为它发生在您的程序编译之前。而且即使你支持全部 32 位整数，你的表也只有 256kb。
我在高性能环境中使用了策略 1，效果很好。我进一步提高了搜索的性能，因为我知道下一个要采用的 sqrt 可能与前一个 sqrt 很接近（上下文是图形的），它产生了惊人的性能差异。
@R.. ：我认为 (1) 不会比 sqrt 快；对 999999 个整数的列表进行二进制搜索很可能比 sqrt 慢！
@Nawaz：鉴于您显然足够关心这个问题，那么在谴责它之前对其进行基准测试怎么样。很大程度上取决于您的确切硬件......
@Nawaz 和 R：实际上我现在实施得更好了。它不是无枝的，但它把其他所有东西都吹到了水里：gist.github.com/3481607

【解决方案2】：

我会尝试Fast Inverse Square Root 技巧。

这是一种在没有任何分支的情况下获得 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

正如预期的那样，Fast 计算的性能比 Int 计算好得多。

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

【讨论】：

这是浮点数，但 nawaz 只需要整数值。
@Saeed：整数可以简单地转换为浮点数（来回），我很好奇这种方法的适用性。它当然是我能想到的唯一没有预计算表的无分支方法。在那之后......我想我们可以进行基准测试:)？
是的，它可以，但我认为方法（如我引用的文章）更快（因为它们只关心整数部分）但是应该以这种方式进行基准测试。
@Saeed：完成，正如预期的那样，快速逆向技巧表现更好，我猜无分支是有回报的
@Saeed：通常你已经得到了所有的包含。我想知道-ffast-math 会为sqrt 付出什么。

【解决方案3】：

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

【讨论】：

【解决方案4】：

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

在使用正确设置 (-march=native -m64 -O3) 进行分析后，上述速度快了很多。

好吧，有点老问题了，但“最快”的答案还没有给出。最快的（我认为）是二元平方根算法，在this Embedded.com article 中有详细说明。

基本上是这样的：

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 编译选项，时间分别为 2100 毫秒和 3100 毫秒。请注意，这甚至没有使用一行汇编程序，因此它可能仍然更快。

上述代码适用于 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

【讨论】：

试试我的无分支变体，看看它是否更快。
@R.：不，它慢了大约 3 倍。
您的编译器是否使用cmov 作为您的版本？
你也可以试试这个带有（实现定义的）签名右移的简化版本：for (s=squares, i=128; i; i=i>>1) s += s[i]-x-1>>31 & i; return s-squares;
@R.：不，它不使用cmov。此外，手动展开循环实际上快了 20% 左右。这是两个版本的 asm 输出（注意我制作的是 32 位版本）：gist.github.com/3481749 完整的 32 位版本可以在这里下载：gist.github.com/3481770

【解决方案5】：

要进行整数 sqrt，您可以使用牛顿方法的这种特殊化：

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)) 中。它也根本不使用浮点数，它也适用于任意精度的整数。

【讨论】：

我在 iPhone 硬件上试过那个，但由于 'b = N / a' 操作，它似乎很慢。
@AndrewTomazos — 不幸的是，您的函数无法返回 N ∈ { 0, 3, 8, 48, 63, 120, 143, ... } 的正确答案。
这仅在整数除法比 FP sqrt 更有效时才有用。或者，如果您需要处理太大而无法精确浮点表示的整数。现代 x86 硬件 (Intel Haswell) 可以在大约 25 个时钟周期延迟内转换为浮点数并返回，并执行双精度 FP sqrt。单个 32 位整数除法具有 22 到 29 个周期的延迟。从吞吐量的角度来看，整数牛顿迭代也表现不佳。

【解决方案6】：

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

如果：

数字范围有限
内存消耗并不重要
应用程序启动时间并不重要

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

所有这些条件都非常适合我的程序。特别是，int[10000000] 数组将消耗 40MB。

您对此有何看法？

【讨论】：

数组太大，缓存未命中可能比 sqrt 计算代价更高。

【解决方案7】：

在许多情况下，甚至不需要精确的整数 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) 相差不大，因此在近似精度足够时很有用。

【讨论】：

【解决方案8】：

在我的带有 gcc 和 -ffast-math 的计算机上，将 32 位整数转换为浮点数并使用 sqrtf 每 10^9 次操作需要 1.2 秒（没有 -ffast-math 需要 3.54 秒）。

以下算法每 10^9 使用 0.87 秒，但会牺牲一些准确性：尽管 RMS 误差仅为 0.79，但误差可能高达 -7 或 +1：

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 更快，至少使用 -ffast-math。

【讨论】：

两级表在这里也可能表现良好。

【解决方案9】：

如果你不介意一个近似值，我拼凑的这个整数 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;
}

它使用这篇Wikipedia 文章中描述的算法。在我的机器上，它几乎是 sqrt 的两倍 :)

【讨论】：

从技术上讲，这打破了严格的别名规则。在最近的 gcc (4.9) 下似乎不会引起问题，但符合要求的做法是union { float f; int32_t x } v; v.f = (float) x; v.x -= ... return (int)((float)v.x);。
学习一些东西对我很有帮助，但是你们是如何衡量速度的呢？我在 www 上找到的所有 sqrt(int) 函数都比 std::sqrt 慢。根据我的测量。顺便说一句，你的函数慢了 27%，准确率低了 3%。比标准

【解决方案10】：

这很短，以至于它 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);
}

英特尔内部版本（在 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)（如 cmets 中的 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);
}

【讨论】：

在将值复制到其中之前，您真的需要将 xmm0 归零吗？此外，内联可以禁用对周围代码的优化。使用来自here 的内置函数（即__builtin_ia32_cvtsi2ss、__builtin_ia32_sqrtss__builtin_ia32_cvtss2si）怎么样？在使用内联 asm 时，less is more.
英特尔内部版本不会向下舍入 - 它会四舍五入到最接近的 int。要解决此问题，请将t 添加到其最后一行：return _mm_cvtt_ss2si(xmm0);。这些比我的机器上的 sqrt() 快 5-6 倍）但是由于浮点数的舍入错误，当 num >= 16,785,407 时开始出现错误答案。要在 GCC Intrinsic 版本中解决此问题，请将第一行更改为 __v2df xmm0 = {0, 0}; 并将每个 ss 替换为 sd（警告：将速度减半）。由于某种原因，我在英特尔的Intrinsics Guide 中没有看到_mm_cvt_si2sd()。
您确定cvtsi2ss %1 确实安全吗？假设%1 是edi，那么rdi 的上半部分可能有一些位。也许您在自己的测试中不会注意到它，但原则上它可能会发生。
@J.Schultke 很好，我将初始化i

【解决方案11】：

以下解决方案计算整数部分，即 floor(sqrt(x)) 精确，没有舍入误差。

其他方法的问题

使用float 或double 既不便携也不精确
@orlp 的isqrt 给出了像isqrt(100) = 15 这样的疯狂结果
基于大型查找表的方法超过 32 位是不实用的
使用快速逆sqrt非常不精确，你最好使用sqrtf
牛顿的方法需要昂贵的整数除法和良好的初始猜测

我的方法

我的基于bit-guessing approach proposed on Wikipedia。不幸的是，维基百科上提供的伪代码有一些错误，所以我不得不做一些调整：

// 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 次迭代。

compile output 只有大约 20 条指令。

性能

我通过统一生成输入，将我的方法与float-bases 的方法进行了基准测试。请注意，在现实世界中，大多数输入会比std::numeric_limits<...>::max() 更接近于零。

对于uint32_t，这比使用std::sqrt(float) 的性能要差于25x
对于uint64_t，这比使用std::sqrt(double) 的性能差于30x

准确度

与使用浮点数学的方法不同，此方法始终非常准确。

使用 sqrtf 可能会在 [2²⁸, 2³²) 范围内提供不正确的舍入。例如sqrtf(0xffffffff) = 65536，当平方根实际上是65535.99999。
在 [2⁶⁰, 2⁶⁴) 范围内，双精度无法始终如一地工作。例如sqrt(0x3fff...) = 2147483648，当平方根实际上是2147483647.999999。

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

结论

正如我所说，这是唯一能正确处理所有输入、避免整数除法且不需要查找表的解决方案。总之，如果您需要一种独立于精度并且不需要庞大的查找表的方法，那么这是您唯一的选择。它可能在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^floor(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 作为初始猜测快大约五倍。

【讨论】：

【解决方案12】：

或者只是做一个二分查找，不能写一个更简单的版本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;
}

【讨论】：

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

其他方法的问题

我的方法

性能

准确度

结论

使用牛顿法的替代方法

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