为什么pow(int,int)这么慢?

Fan*_*ang 48 c++ performance cmath pow

我一直在做一些项目Euler练习,以提高我对C++的知识.

我写了以下函数:

int a = 0,b = 0,c = 0;

for (a = 1; a <= SUMTOTAL; a++)
{
    for (b = a+1; b <= SUMTOTAL-a; b++)
    {
        c = SUMTOTAL-(a+b);

        if (c == sqrt(pow(a,2)+pow(b,2)) && b < c)
        {
            std::cout << "a: " << a << " b: " << b << " c: "<< c << std::endl;
            std::cout << a * b * c << std::endl;
        }
    }
}
Run Code Online (Sandbox Code Playgroud)

这计算在17毫秒.

但是,如果我改变了这条线

if (c == sqrt(pow(a,2)+pow(b,2)) && b < c)
Run Code Online (Sandbox Code Playgroud)

if (c == sqrt((a*a)+(b*b)) && b < c)
Run Code Online (Sandbox Code Playgroud)

计算在2毫秒内完成.是否有一些明显的实现细节pow(int, int)我缺少哪个使得第一个表达式计算速度慢得多?

Rin*_*g Ø 68

pow() 使用真实的浮点数并在引擎盖下使用公式

pow(x,y) = e^(y log(x))
Run Code Online (Sandbox Code Playgroud)

计算x^y.将int被转换为double之前调用pow.(log是自然对数,基于e)

x^2pow()因此使用比慢x*x.

根据相关评论进行编辑

  • 使用pow整数指数甚至可能会产生不正确的结果(PaulMcKenzie)
  • 除了使用双重类型的数学函数外,pow还有函数调用(而x*x不是)(jtbandes)
  • 事实上,许多现代编译器都会使用常量整数参数来优化pow,但这不应该依赖.

  • 我的观点是,C++标准中的任何内容都不能保证`pow`能给出精确的结果,即使是整数指数也是如此.使用查找表或其他一些保证没有舍入错误的方法让生活变得愉快和计算能量,无论使用哪个库. (6认同)
  • 它不仅速度慢,而且可以获得不精确的答案,即使对于整数基数和指数也是如此. (5认同)
  • @YanZhou - 它不会总是给出确切的结果,否则[这永远不会被问到](http://stackoverflow.com/questions/25678481/why-does-pown-2-return-24-when-n -5与 - 我的编译器和-OS?noredirect = 1&LQ = 1) (3认同)
  • @PaulMcKenzie正如我所说的,信誉良好的libm就属于这种情况.不是每个libm都准确.据我所知,AMD libm,英特尔libimf,OpenLibm,BSD libm及其衍生产品(例如macOS中的产品)都将为您提供`pow(5,2)== 25`,您引用的示例.GNU libc是linux上使用最广泛的,但它并没有使它成为黄金标准 (3认同)

Pet*_*des 37

你已经选择了一种最慢的方法来检查

c*c == a*a + b*b   // assuming c is non-negative
Run Code Online (Sandbox Code Playgroud)

这会编译为三个整数乘法(其中一个可以从循环中提升).即使没有pow(),你仍然会转换double并采用平方根,这对吞吐量来说很糟糕.(还有延迟,但分支预测+现代CPU上的推测性执行意味着延迟不是这里的一个因素).

Intel Haswell的SQRTSD指令的吞吐量为每8-14个周期一个(源:Agner Fog的指令表),所以即使你的sqrt()版本保持FP sqrt执行单元饱和,它仍然比我发出的gcc慢4倍(下面).


b < c条件的一部分变为false时,您还可以优化循环条件以摆脱循环,因此编译器只需执行该检查的一个版本.

void foo_optimized()
{ 
  for (int a = 1; a <= SUMTOTAL; a++) {
    for (int b = a+1; b < SUMTOTAL-a-b; b++) {
        // int c = SUMTOTAL-(a+b);   // gcc won't always transform signed-integer math, so this prevents hoisting (SUMTOTAL-a) :(
        int c = (SUMTOTAL-a) - b;
        // if (b >= c) break;  // just changed the loop condition instead

        // the compiler can hoist a*a out of the loop for us
        if (/* b < c && */ c*c == a*a + b*b) {
            // Just print a newline.  std::endl also flushes, which bloats the asm
            std::cout << "a: " << a << " b: " << b << " c: "<< c << '\n';
            std::cout << a * b * c << '\n';
        }
    }
  }
}
Run Code Online (Sandbox Code Playgroud)

-O3 -mtune=haswell用内循环编译(使用gcc6.2 )代码.查看Godbolt编译器资源管理器的完整代码.

# a*a is hoisted out of the loop.  It's in r15d
.L6:
    add     ebp, 1    # b++
    sub     ebx, 1    # c--
    add     r12d, r14d        # ivtmp.36, ivtmp.43  # not sure what this is or why it's in the loop, would have to look again at the asm outside
    cmp     ebp, ebx  # b, _39
    jg      .L13    ## This is the loop-exit branch, not-taken until the end
                    ## .L13 is the rest of the outer loop.
                    ##  It sets up for the next entry to this inner loop.
.L8:
    mov     eax, ebp        # multiply a copy of the counters
    mov     edx, ebx
    imul    eax, ebp        # b*b
    imul    edx, ebx        # c*c
    add     eax, r15d       # a*a + b*b
    cmp     edx, eax  # tmp137, tmp139
    jne     .L6
 ## Fall-through into the cout print code when we find a match
 ## extremely rare, so should predict near-perfectly
Run Code Online (Sandbox Code Playgroud)

在Intel Haswell上,所有这些指令均为1 uop.(并且cmp/jcc将宏 - 融合成比较和分支的uops.)因此,这是10个融合域uops,它可以每2.5个周期在一次迭代中发出.

Haswell imul r32, r32以每个时钟一次迭代的吞吐量运行,因此内部循环内的两个乘法不会使端口1饱和,每2.5c两次乘以.这留下了空间来吸收来自ADD和SUB窃取端口1的不可避免的资源冲突.

我们甚至没有接近任何其他执行端口瓶颈,因此前端瓶颈是唯一的问题,这应该在Intel Haswell及更高版本上每2.5个周期运行一次.

循环展开可以帮助减少每次检查的uop数量.例如,用于lea ecx, [rbx+1]计算下一次迭代的b + 1,因此我们可以imul ebx, ebx不使用MOV使其具有非破坏性.


强度降低也是可能的:鉴于b*b我们可以尝试在(b-1) * (b-1)没有IMUL的情况下进行计算. (b-1) * (b-1) = b*b - 2*b + 1,也许我们可以做一个lea ecx, [rbx*2 - 1]然后从中减去b*b.(没有寻址模式可以减去而不是添加.嗯,也许我们可以保留-b在寄存器中,并向上计数到零,因此我们可以使用ECX中的lea ecx, [rcx + rbx*2 - 1]更新b*b,-b在EBX中给出).

除非你真的遇到了IMUL吞吐量的瓶颈,否则这可能最终会带来更多的麻烦,而不是一场胜利.看看编译器在C++源代码中实现这种强度降低的效果可能会很有趣.


您也可以使用SSE或AVX对其进行矢量化b,并行检查4个或8个连续值.由于命中率非常罕见,因此您只需检查8中的任何一个是否有命中,然后在极少数情况下找出匹配中的哪一个.

有关更多优化内容,另请参阅标记wiki.

  • @Fang:对于C来说绝对不是这样.所有优化编译器只看到实现代码所需的数据移动.如果有的话,限制变量的范围(并为不同的循环使用单独的变量)将有助于编译器看到它们是独立的,并且在循环之后死亡. (2认同)