30

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

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

正在使用 MSVC++ 2008 编译器，以供参考……尽管我认为 sqrt 不会增加太多开销。

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

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

4

6 回答 6

34

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

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

  2. 查找表：要么仅用于迭代的起点，要么与插值相结合以使您一直到达那里。

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

11年后的你好。

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

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

  2. 看一下英特尔优化参考手册的第 15.12.3 节，它描述了近似方法，带有矢量化指令，这可能也可以很好地转化为其他架构。

于 2010-04-14T13:42:49.957 回答
15

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

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

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

于 2010-04-14T15:29:09.040 回答
8

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

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

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

于 2010-04-14T13:32:16.757 回答
3

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

于 2010-04-14T13:44:00.050 回答
2

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

你可以在这里看到大量替代品的描述。

最好的是这个魔术片段：

double inline __declspec (naked) __fastcall sqrt(double n)
{
    _asm fld qword ptr [esp+4]
    _asm fsqrt
    _asm ret 8
}

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

于 2013-11-23T16:37:20.873 回答
1

这是一个查找表只有 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;
}
于 2015-11-11T18:12:21.727 回答