是否可以推出更快的 sqrt 版本

Posted 2023-02-16

【中文标题】是否可以推出更快的 sqrt 版本

【英文标题】：Is it possible to roll a significantly faster version of sqrt

【发布时间】：2010-04-14 13:29:13

【问题描述】：

在我正在分析的应用程序中，我发现在某些情况下，此函数能够占用总执行时间的 10% 以上。

多年来，我看到了有关使用偷偷摸摸的浮点技巧更快地实现 sqrt 的讨论，但我不知道这些东西在现代 CPU 上是否已经过时。

正在使用 MSVC++ 2008 编译器，供参考...虽然我认为 sqrt 不会增加太多开销。

有关modf 函数的类似讨论，另请参见此处。

编辑：作为参考，this 是一种广泛使用的方法，但它实际上更快吗？现在 SQRT 到底有多少个周期？

【问题讨论】：

你能发布一些代码吗？优化 sqrt 的最好方法是去掉它，或者至少减少对它的调用次数，这可能是可能的。

代码是来自第 3 方的长而复杂的软体物理建模。没有几个执行 sqrt 的内部循环可以使用 length^2 代替 length :)

单精度还是双精度？您需要什么精度？

不要使用“快速反平方根”。如果您愿意接受一个近似值，硬件rsqrtss（近似倒数平方根）要快得多。

请参阅***.com/questions/31555260/… 以获得（近似）rsqrtps + 牛顿迭代的体面版本，对于单精度 float，给出 +/-2ulp。请参阅***.com/questions/32002277/… 了解-mrecip 编译器优化，它应该控制rsqrt、but doesn't seem to actually do so (only rcp for 1/x) 的自动使用。

【参考方案1】：

是的，即使没有诡计也是可能的：

牺牲精度换取速度：sqrt 算法是迭代的，用更少的迭代重新实现。

查找表：要么只是迭代的起点，要么与插值相结合，让您一直到达那里。

缓存：您是否总是对相同的有限值集进行排序？如果是这样，缓存可以很好地工作。我发现这在图形应用程序中很有用，在这些应用程序中，要为许多相同大小的形状计算相同的东西，因此可以有效地缓存结果。

---

11年后的你好。

考虑到这仍然会偶尔获得投票，我想我会添加一个关于性能的注释，现在它甚至比那时更受内存访问的限制。在优化这样的事情时，您绝对必须使用现实的基准（理想情况下，您的整个应用程序） - 您的应用程序的内存访问模式将对查找表和缓存等解决方案产生巨大影响，并且只需比较优化版本的“周期”会让您误入歧途：将程序时间分配给单个指令也非常困难，您的分析工具可能会在此处误导您。

在相关说明中，如果适合您的用例，请考虑使用 simd/矢量化指令来计算平方根，例如 _mm512_sqrt_ps 或类似指令。

看看英特尔optimisation reference manual 的第 15.12.3 节，它描述了近似方法，带有矢量化指令，这可能也可以很好地转化为其他架构。

【讨论】：

我总是很难相信手动执行即使是少量的迭代也会比内置的 SQRT 指令更快......但我猜 SQRT 并不神奇，它仍然如此里面的迭代。

您有任何衡量标准...您看到了多少改进？

里程数因使用情况而异 :) 您确实必须分析自己的使用场景才能查看有效的方法。关于 fsqrt 指令，您可能仍然可以使用它，但是通过不检查错误条件来加快速度：这取决于您的编译器正在生成的汇编程序。

使用 quake sqrt 算法。它的魔力不在于迭代，而在于初始值的选择。但是 O3 级别的 gcc 实现了它，所以我的建议是使用它。

@rytis 著名的地震算法适用于平方根的倒数 (1/sqrt(x))，而不适用于 sqrt(x)

【参考方案2】：

这里有一个很棒的比较表： http://assemblyrequired.crashworks.org/timing-square-root/

长话短说，SSE2 的 ssqrts 比 FPU fsqrt 快大约 2 倍，近似 + 迭代大约是 4 倍（总共 8 倍）。

另外，如果您尝试使用单精度 sqrt，请确保您确实得到了这样的结果。我听说至少有一个编译器会将浮点参数转换为双精度，调用双精度 sqrt，然后再转换回浮点数。

【讨论】：

那个链接现在失效了。它仍然可以在 archive.org 上找到：web.archive.org/web/20210208132927/http://…

【参考方案3】：

您很可能通过更改您的算法而不是更改它们的实现来获得更多的速度改进：尽量少调用sqrt()而不是加快调用速度。 （如果您认为这是不可能的 - 您提到的 sqrt() 的改进就是：用于计算平方根的 算法 的改进。）

由于经常使用它，您的标准库对sqrt() 的实现很可能在一般情况下几乎是最佳的。除非您有一个受限域（例如，如果您需要较低的精度），否则算法可以采取一些捷径，否则不太可能有人提出更快的实现。

请注意，由于该函数占用了您 10% 的执行时间，即使您设法提出一个只占用 std::sqrt() 75% 时间的实现，这仍然只会使您的执行时间减少2.5%。对于大多数应用程序，用户甚至不会注意到这一点，除非他们使用手表进行测量。

【讨论】：

这方面的一个例子可能是：如果你想找到离你最近的东西，你可以比较距离的平方而不是取每个的平方，因为距离*距离 > 距离。或者您可以运行一个预处理步骤，预先计算所有内容的距离对。 （显然我在想象某种 2-D 或 3-D 数据结构）。

+1 意识到对使用频率较低的一段代码的重大改进最终会在大图中得到几乎可以忽略不计的改进。

为什么人们似乎认为原始发帖人是智障？只需回答他们的问题，而不是告诉他们不应该尝试做他们正在做的事情。也许他们有充分的理由去做他们正在做的事情。 10% 的代码时间对于一个功能来说是一大块时间，如果操作简单，则值得优化。我不敢相信这种无益的回应得到了如此多的支持。

实际上，我发现这个回复经过深思熟虑，很好地解释了各种可能性，并权衡了改进实现与选择可能更好的不同方法的利弊。 +1 提供真正有用且组合良好的答案。

“你没有问正确的问题”实际上应该是评论而不是答案，因为它没有回答所提出的问题。没错，但问题很具体。

【参考方案4】：

您需要sqrt 的准确度如何？您可以很快得到合理的近似值：请参阅 Quake3 出色的 inverse square root 函数以获得灵感（请注意，代码是 GPL 的，因此您可能不想直接集成它）。

【讨论】：

它有多快？如果 1/sqrt 不行，又需要 sqrt，额外除法是不是还是比普通版快？

@John：他们有。毕竟，这就是他们创建该功能的原因。但这并不意味着它会对您的情况有所帮助。

@John：我无法在您的系统上进行测试，并且在引用函数中进行浮点修改的系统变化太重要了，无法忽略。

所谓的“快速反平方根”在现代硬件上不是“快速”。在过去 10 年设计的几乎所有处理器上，都有更快的替代方案。在许多情况下，硬件平方根指令会更快。许多有更快的硬件逆平方根估计（SSE 上的rsqrtss，ARMv7 上的rsqrte，等等）。

@jemfinch ... 它是 VS2008，在典型的 PC 上运行。不只是我的电脑。

【参考方案5】：

不知道你是否解决了这个问题，但我以前读过它，似乎最快的做法是将sqrt 函数替换为内联汇编版本；

您可以看到大量备选方案的描述here。

最好的是这个魔法sn-p：

double inline __declspec (naked) __fastcall sqrt(double n)

    _asm fld qword ptr [esp+4]
    _asm fsqrt
    _asm ret 8
 

它比具有相同精度的标准 sqrt 调用快约 4.7 倍。

【讨论】：

您能否建议在 GCC 中的实施？对于浮点类型应该如何修改？谢谢！

那为什么编译器不这样做呢？

我真的不知道。

什么鬼？这不会自动矢量化，并且在内联时可能会中断（esp+4 可能不是n，并且在内联汇编中有一个ret）。它强制您的输入通过存储/重新加载。这是在 Haswell 上 fsqrt 的 10-23 个周期之上的另外约 5 个延迟周期。 sqrtpd（双）是 16c 延迟，而sqrtps（单）是 11c 延迟。您的编译器必须非常糟糕才能赢得胜利（也许您在没有优化的情况下编译？）另外，ret 8?这是一个 3uop 指令。

@will 测试不准确，因为例如他将sqrt(int) 与std::sqrt(double) 进行比较，这是错误的，内部 std::sqrt 算法速度对于所有参数类型都不相同。顺便说一句，我重新实现了他的测试，只有 2 或 3 个变体更快，但与标准变体相比，准确性大大降低..

【参考方案6】：

这是一个快速的方法，查找表只有 8KB。错误约为结果的 0.5%。您可以轻松地放大表格，从而减少错误。运行速度比常规 sqrt() 快约 5 倍

// LUT for fast sqrt of floats. Table will be consist of 2 parts, half for sqrt(X) and half for sqrt(2X).
const int nBitsForSQRTprecision = 11;                       // Use only 11 most sagnificant bits from the 23 of float. We can use 15 bits instead. It will produce less error but take more place in a memory. 
const int nUnusedBits   = 23 - nBitsForSQRTprecision;       // Amount of bits we will disregard
const int tableSize     = (1 << (nBitsForSQRTprecision+1)); // 2^nBits*2 because we have 2 halves of the table.
static short sqrtTab[tableSize]; 
static unsigned char is_sqrttab_initialized = FALSE;        // Once initialized will be true

// Table of precalculated sqrt() for future fast calculation. Approximates the exact with an error of about 0.5%
// Note: To access the bits of a float in C quickly we must misuse pointers.
// More info in: http://en.wikipedia.org/wiki/Single_precision
void build_fsqrt_table(void)
    unsigned short i;
    float f;
    UINT32 *fi = (UINT32*)&f;

    if (is_sqrttab_initialized)
        return;

    const int halfTableSize = (tableSize>>1);
    for (i=0; i < halfTableSize; i++)
         *fi = 0;
         *fi = (i << nUnusedBits) | (127 << 23); // Build a float with the bit pattern i as mantissa, and an exponent of 0, stored as 127

         // Take the square root then strip the first 'nBitsForSQRTprecision' bits of the mantissa into the table
         f = sqrtf(f);
         sqrtTab[i] = (short)((*fi & 0x7fffff) >> nUnusedBits);

         // Repeat the process, this time with an exponent of 1, stored as 128
         *fi = 0;
         *fi = (i << nUnusedBits) | (128 << 23);
         f = sqrtf(f);
         sqrtTab[i+halfTableSize] = (short)((*fi & 0x7fffff) >> nUnusedBits);
    
    is_sqrttab_initialized = TRUE;


// Calculation of a square root. Divide the exponent of float by 2 and sqrt() its mantissa using the precalculated table.
float fast_float_sqrt(float n)
    if (n <= 0.f) 
        return 0.f;                           // On 0 or negative return 0.
    UINT32 *num = (UINT32*)&n;
    short e;                                  // Exponent
    e = (*num >> 23) - 127;                   // In 'float' the exponent is stored with 127 added.
    *num &= 0x7fffff;                         // leave only the mantissa 

    // If the exponent is odd so we have to look it up in the second half of the lookup table, so we set the high bit.
    const int halfTableSize = (tableSize>>1);
    const int secondHalphTableIdBit = halfTableSize << nUnusedBits;
    if (e & 0x01) 
        *num |= secondHalphTableIdBit;  
    e >>= 1;                                  // Divide the exponent by two (note that in C the shift operators are sign preserving for signed operands

    // Do the table lookup, based on the quaternary mantissa, then reconstruct the result back into a float
    *num = ((sqrtTab[*num >> nUnusedBits]) << nUnusedBits) | ((e + 127) << 23);
    return n;

【讨论】：

我预计缓存未命中会比sqrt 更糟糕地用于实际代码，而不是在微基准测试中测试 sqrt 吞吐量的代码。 x86 在硬件中具有非常快的sqrt，甚至更快的rsqrt（近似倒数 sqrt），您可以使用 newton-raphson 步骤对其进行优化。