为什么 pow(int, int) 函数如此缓慢？

Question

我一直在做一些Project 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;
        }
    }
}

这需要17毫秒来计算。

然而，如果我更改这行代码

if (c == sqrt(pow(a,2)+pow(b,2)) && b < c)

到

if (c == sqrt((a*a)+(b*b)) && b < c)

计算需要2毫秒的时间。我是否错过了一些明显的pow(int，int)实现细节，导致第一个表达式计算速度如此之慢？

- Fang

---

"a*a" 可能是一条指令。"pow" 至少需要进行函数调用，以及函数本身的计算工作。 - jtbandes

这需要在17毫秒内计算。首先，pow是一个浮点函数。其次，只有在运行经过优化的版本时，发布函数的运行时间才有意义。如果你运行未经优化的“调试”版本，则时间没有意义。最后但同样重要的是，如果指数是整数，请勿使用pow。 - PaulMcKenzie

这个评论 或许会对您有所帮助。正如ringo所说，它既是一个库调用，也是一个“过度强大”的函数。 - Zeta

如果不考虑c*c溢出的情况，使用c*c = a*a + b*b会更快，因为乘法（特别是整数乘法）比开方运算更快。 - Roel Schroeven

另一个提示：如果您反转条件语句的顺序，可以加快速度：if (b < c && c*c == a*a + b*b) 而不是 if (c*c == a*a + b*b && b < c)：b < c 是一项快速操作，使程序在不需要时跳过相对较慢的计算。 - Roel Schroeven

Answer 1

pow() 使用实数浮点数，并在幕后使用以下公式：

pow(x,y) = e^(y log(x))

计算x^y时，调用pow前将int转换为double。（log是以e为底的自然对数）

即使使用整数指数，使用pow()也比x*x慢。

• 使用double类型的数学函数pow，即使是整数指数，可能会产生不正确的结果（PaulMcKenzie）
• 除了使用带double参数的数学函数外，pow是一个函数调用（而x*x则不是）（jtbandes）
• 许多现代编译器实际上会优化带有常量整数参数的pow，但不能依赖于此。

- Déjà vu

---

不仅速度较慢，即使对于整数的底数和指数，也可能得到不精确的答案。 - PaulMcKenzie

@YanZhou -- 它并不总是会给出确切的结果，否则这个问题就不会被提出。 - PaulMcKenzie

正如我所说，这是有声誉的libm的情况。并非所有的libm都能给出精确结果。据我所知，AMD libm、Intel libimf、OpenLibm、BSD libm以及其衍生品（例如macOS中的那个）都可以给出你引用的pow(5, 2) == 25的结果。GNU libc在Linux上被广泛使用，但这并不意味着它是黄金标准。 - Yan Zhou

@PaulMcKenzie 有一个更正，GNU libc 库也会给出相同的结果。你引用的那篇文章可能是在 Windows 上使用 GCC+Code Blocks，而 MS libm 库则不太可靠。 - Yan Zhou

我的观点是，即使使用整数指数，C++标准中也没有保证pow函数会给出精确的结果。为了让生活更愉快，可以使用查找表或其他能够保证没有舍入误差的方法来计算pow函数的值，而不受所使用的库的影响。 - PaulMcKenzie

Answer 2

你选择了最慢的方式来进行检查。

c*c == a*a + b*b   // assuming c is non-negative

这段代码编译后会产生三个整数乘法（其中一个可以提出循环）。即使没有使用pow()函数，你仍然需要将其转换为double类型并进行平方根运算，这对于吞吐量来说非常糟糕。（对于延迟也是如此，但现代CPU上的分支预测和推测执行意味着延迟在这里不是一个因素）。

Intel Haswell的SQRTSD指令的吞吐量为每8-14个周期执行一次（source: Agner Fog's instruction tables），因此即使你的sqrt()版本可以保持FP sqrt执行单元饱和，它仍然比我得到的gcc编译版本慢4倍（下面附有代码）。

---

您还可以优化循环条件，当条件中的b < c不再成立时跳出循环，这样编译器只需要进行一次检查版本。

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

这段代码（使用gcc6.2的-O3 -mtune=haswell编译）生成了以下内部循环的代码。请在 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

在英特尔Haswell架构中，所有这些指令都是1个微操作。 (而cmp / jcc对宏融合成比较和分支微操作。) 因此，这是10个融合域微操作，每2.5个周期可以迭代一次。

Haswell以每个时钟周期的速度运行imul r32，r32，因此内部循环中的两个乘法不会以每2.5c的速度饱和端口1。这为吸收ADD和SUB从端口1窃取的不可避免的资源冲突留出了空间。

我们甚至没有接近任何其他执行端口瓶颈，因此前端瓶颈是唯一的问题，在英特尔Haswell及更高版本上应该以每2.5个周期的速度运行。

循环展开可以帮助减少每次检查的uops数量。例如，使用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，并向零计数，这样我们就可以使用lea ecx，[rcx + rbx*2-1]来更新ECX中的b*b，给定EBX中的-b）。

除非您实际上受到IMUL吞吐量的限制，否则这可能会导致更多uops，并且不会获胜。 在C ++源代码中使用此强度降低技术，可能会很有趣。

---

你可能也可以使用SSE或AVX将其向量化，同时并行检查4或8个连续的b值。由于命中非常罕见，只需检查这8个值中是否有任何一个命中，然后在极少数情况下确定是哪一个命中。

另请参阅x86标签wiki以获取更多优化内容。

- Peter Cordes

---

@étale-cohomology：自动向量化之所以有效，是因为if()很少被执行，而编译器并不知道这一点。如果完整向量中没有多个元素需要执行if()，那么进行向量化计算可能会带来损失。(虽然你可以通过扫描PMOVMSKB位图而不是重新计算来查看要打印的元素)。总之，是的，自动向量化可以工作，但是由于条件行为，这个问题对于编译器来说并不简单。编译器技术还没有达到我希望他们能够解决这个问题的地步。 - Peter Cordes

@PeterCordes：我们的教授总是告诉我们，在循环外声明变量可以提高性能。当然，我读过不同的观点，这取决于编译器和编程语言。 - Fang

@Fang：对于C语言来说，这绝对是不正确的。所有优化编译器都只会看到实现代码所需的数据移动。如果有的话，限制变量的作用域（并为不同的循环使用单独的变量）将有助于编译器确定它们是独立的，并且在循环结束后就不再被使用了。 - Peter Cordes

@Zboson：使用O(log2(n))方法，pow(x,4)只需要2步，而不是4步。您到底建议使用哪种瓶颈吞吐量的O(n)方法？肯定不是为了让乘数保持忙碌而重复执行相同的乘法吧？此外，您是否打算在讨论pow(double,int)的其他线程中留下该评论？ - Peter Cordes

让我们在聊天中继续这个讨论。 - Z boson