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

Question

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

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

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

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

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

Answer 1

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

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

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

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

---

11年后的你好。

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

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

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

Answer 2

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

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

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

Answer 3

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

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

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

Answer 4

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

Answer 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 倍。

Answer 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;
}