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

Question

我一直在进行一些项目 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) 是否有一些明显的实现细节我遗漏了，这使得第一个表达式的计算速度如此之慢？

Answer 1

pow() 使用实浮点数并在后台使用公式

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

计算x^y。 int 在调用 pow 之前转换为 double。 （log是自然对数，基于e）

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

根据相关cmets进行编辑

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

Answer 2

您选择了一种最慢的检查方法

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

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

英特尔 Haswell 的 SQRTSD 指令每 8-14 个周期 (source: Agner Fog's instruction tables) 的吞吐量为 1，因此即使您的 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）以使用此内部循环进行编码。完整代码见the Godbolt compiler explorer。

# 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

在 Intel Haswell 上，所有这些指令都是 1 uop。 （并且 cmp/jcc 将宏融合成比较和分支微指令。）所以这是 10 个融合域微指令，which can issue at one iteration per 2.5 cycles。

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

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

循环展开可以帮助减少每次检查的微指令数。例如使用lea ecx, [rbx+1] 为下一次迭代计算 b+1，因此我们可以在不使用 MOV 的情况下使用imul ebx, ebx 使其成为非破坏性的。

---

强度降低也是可能的：给定b*b，我们可以尝试在没有IMUL 的情况下计算(b-1) * (b-1)。 (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，给定@987654344 @ 在 EBX 中）。

除非您实际上在 IMUL 吞吐量上遇到瓶颈，否则这最终可能会占用更多微指令而不是胜利。看看编译器在 C++ 源代码中强度降低的效果如何可能会很有趣。

---

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

有关更多优化内容，另请参阅x86 标签 wiki。