有没有可能开发出一个速度显著更快的sqrt函数?

34
在我分析的一个应用程序中,我发现在某些情况下,该函数能够占据总执行时间的10%以上。
多年来,我看到了有关使用狡猾的浮点技巧实现更快的sqrt函数的讨论,但我不知道这些东西在现代CPU上是否已经过时。
参考使用了MSVC++ 2008编译器...虽然我认为sqrt不会增加太多开销。
请参见这里,了解modf函数的类似讨论。
编辑:作为参考,this是一种广泛使用的方法,但它实际上是否更快?SQRT现在需要多少个周期?

4
你能发一些代码吗?优化sqrt的最佳方法是摆脱它,或者至少减少调用它的次数,这可能是可行的。 - IVlad
2
代码又长又复杂,使用了第三方的软体物理建模。不是一些内部循环做平方根,而是可以用长度的平方代替长度的情况 :) - Mr. Boy
3
单精度还是双精度?你需要什么样的准确度? - Paul R
5
不要使用“快速反平方根算法”。如果您可以接受近似值,硬件的rsqrtss(近似倒数平方根)会更快。 - Stephen Canon
2
参见 https://dev59.com/aVwZ5IYBdhLWcg3wWfOl 得到一个(近似)rsqrtps和牛顿迭代的不错版本,对于单精度float可达到+/-2ulp。参见 https://dev59.com/H1wY5IYBdhLWcg3w9r7a 了解 -mrecip 编译器优化,它应该控制自动使用 rsqrt(但实际上只有 1/x 的 rcp)。 - Peter Cordes
显示剩余3条评论
6个回答

39

是的,即使没有欺骗也是可能的:

  1. 为了速度而牺牲准确性:sqrt算法是迭代的,可以重新实现更少的迭代。

  2. 查找表:可以用于迭代的起点,或与插值结合使用以帮助你一直到达目标。

  3. 缓存:你总是对同一组有限的值进行平方根运算吗?如果是,缓存可以很好地工作。我在图形应用程序中发现这很有用,因为相同大小的许多形状正在计算相同的东西,因此结果可以被有用地缓存。


来自未来的11年问候。

考虑到这仍然偶尔会得到投票,我想添加一条关于性能的注释,现在比以前更受内存访问的限制。当优化像这样的东西时,您绝对必须使用逼真的基准测试(理想情况下是整个应用程序)-您应用程序的内存访问模式将对查找表和缓存等解决方案产生巨大影响,仅仅比较您优化版本的“循环”将使您大失所望:还很难将程序时间分配给单个指令,您的分析工具可能会误导您。

  1. 在相关说明上,请考虑使用simd /向量化指令计算平方根,例如_mm512_sqrt_ps或类似产品,如果它们适合您的用例。

  2. 请参阅英特尔的optimisation reference manual第15.12.3节,其中描述了近似方法,使用向量化指令,这些方法可能非常适合其他体系结构。


2
我总是很难相信,即使是手动执行少量迭代,也可能比内置的SQRT指令更快...但是我猜SQRT也不是魔法,它内部仍然进行迭代。 - Mr. Boy
你有任何指标可以衡量吗?你看到了多少改进? - Mr. Boy
4
里程随使用方式而异 :) 您确实需要分析自己的使用场景来确定哪种方法最适合。关于fsqrt指令,您可能仍然可以使用它,但是可以通过不检查错误条件来稍微加快速度:这取决于编译器生成的汇编代码。 - James
使用Quake sqrt算法。它的魔力不在于迭代,而在于选择初始值。然而,gcc在O3级别实现了它,所以我的建议是使用它。 - rytis
@rytis,著名的Quake算法是用于求倒数平方根(1/sqrt(x)),而不是sqrt(x) - James
当然,你是对的,不过 1/(1/sqrt(x)) 可以作为一种妥协。另外,当你提到 LUTs 时,我想知道你是否考虑了在 C++ 中实现 Gal 的精确表格?我很想了解更多信息。 - rytis

18

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

简而言之,SSE2的ssqrts比FPU fsqrt快约两倍,而一种近似+迭代方法比它更快约4倍(总共8倍)。

此外,如果你试图获取单精度平方根,请确保你实际上正在得到它。我听说至少有一个编译器会将浮点参数转换为双精度,调用双精度sqrt,然后再转换回浮点数。


1
那个链接现在已经失效了。但你仍然可以在archive.org上找到它:https://web.archive.org/web/20210208132927/http://assemblyrequired.crashworks.org/timing-square-root/ - Bart

7

改变你的算法比改变它们的实现可能会带来更多的速度提升:尝试减少对sqrt()的调用而不是让调用变得更快。 (如果你认为这不可能-你提到的sqrt()的改进只是用于计算平方根的算法的改进。)

由于sqrt()使用非常频繁,你标准库的实现在一般情况下可能已经接近最优了。除非你有一个受限的域(例如,如果你需要较少的精度),算法可以采取一些捷径,否则很难想出一个更快的实现。

请注意,由于该函数使用了10%的执行时间,即使你成功地想出了一个只需75%的时间就能完成std::sqrt()的实现,这仍然只会将你的执行时间降低 2.5%。对于大多数应用程序用户来说,他们甚至不会注意到这一点,除非他们使用手表来测量。


2
+1 对于这一认识的肯定:对于那些不经常使用的代码进行大幅度改进,最终在“大局”上带来的改善几乎可以忽略不计。 - Thomas Matthews
79
为什么人们似乎默认原作者很愚蠢?回答他们的问题,而不是告诉他们不应该做他们正在做的事情。也许他们有一个很好的理由去做这件事。10% 的代码时间对于一个函数来说是很大的一块时间,如果可以简单地进行优化,那么它就值得优化。我无法相信这个毫无帮助的回复竟然得到了这么多的赞同。 - Joe
3
实际上,我发现这个回复非常经过深思熟虑,解释了可能性,并权衡了改进实现与选择可能更好的不同方法之间的利弊。给予真正有帮助且组织得很好的答案一个+1。 - leeor_net
9
“你没有问对问题”其实应该作为评论而非答案,因为它并没有回答所问的问题。虽然说这个说法是正确的,但问题相当具体。 - Mr. Boy
@Nicolaus:是的,有时会发生这种情况,但我认为这并不常见。在我的经验中,一个函数使用10%的执行时间是很少见的事情,而显著减少这样一个函数的CPU时间对你的应用程序几乎没有什么好处。 - sbi
显示剩余3条评论

3
你需要多精确的平方根函数?你可以快速获得合理的近似值:参考Quake3的优秀的倒数平方根函数(请注意,代码是GPL许可的,因此您可能不想直接集成它)。

它有多快呢?如果1/sqrt不好用,你需要sqrt,那么额外的除法是否仍然比正常版本更快? - Mr. Boy
@John:他们有。毕竟这就是他们创建该函数的原因。但这并不意味着它会在你的情况下有太大帮助。 - sbi
@John:我无法在你的系统上进行测试,而且参考函数中所做的浮点数操作会导致系统变化太大,不能忽略这个重要因素。 - jemfinch
4
所谓的“快速反平方根”在现代硬件上并不“快速”。在过去10年中设计的几乎所有处理器上,都有更快的替代方法。在许多处理器上,硬件平方根指令会更快。许多处理器还具有更快的硬件反平方根估计(SSE上的rsqrtss,ARMv7上的rsqrte等)。 - Stephen Canon
@jemfinch... 是在典型的个人电脑上运行的VS2008。不仅仅是我的电脑。 - Mr. Boy

2

不知道你是否已经解决了这个问题,但我之前读过相关内容,似乎最快的方法是用内联汇编版本替换 sqrt 函数;

你可以在这里看到许多替代方案的介绍。

最好的方法是这个神奇的代码片段:

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

这个方法比标准的sqrt函数快大约4.7倍,精度相同。


12
为什么编译器没有这样做? - Mr. Boy
我真的不知道。 - will
5
什么鬼?那不会自动矢量化,而且当内联时可能会出错(esp+4可能不是 n,并且内联 asm 中有一个 ret)。 它强制要求您的输入通过存储/重新加载来进行。这又增加了大约5个周期的延迟,加上Haswell上10-23个周期的fsqrt延迟。sqrtpd(double)为16个周期的延迟,而sqrtps(single)为11个周期的延迟。 除非您的编译器真的很差(也许您没有进行优化编译?),否则您的编译器不应选择这种方式。此外,ret 8?那是一个3uop指令。 - Peter Cordes
我没有写这篇文章,所以目前还不确定(不查看一下就不知道)他是如何编译和运行测试的。这也让我感到惊讶。 - will
1
@will的测试结果不准确,因为他比较了sqrt(int)std::sqrt(double),这是错误的。实际上,内部的std::sqrt算法速度对于所有参数类型并不相同。顺便说一句,我重新实现了他的测试,只有2或3个变量略微更快,但与std变体相比,精度大大降低了。 - metablaster
显示剩余3条评论

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

4
我预计缓存未命中会使得在实际代码中使用该方法比使用sqrt更糟糕,尤其是针对测试sqrt吞吐量的微基准测试。x86架构具有非常快速的硬件sqrt,甚至还有更快的rsqrt(近似倒数平方根),你可以通过牛顿-拉弗逊步骤来进一步优化。 - Peter Cordes

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接