为什么 pow(int, int) 这么慢?
Why is pow(int, int) so slow?
我一直在做一些项目 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)
实现细节我遗漏了,这使得第一个表达式的计算速度如此之慢?
pow()
适用于实浮点数,并在幕后使用公式
pow(x,y) = e^(y log(x))
计算x^y
。在调用 pow
之前,int
被转换为 double
。 (log
为自然对数,以e为基础)
x^2
使用 pow()
因此比 x*x
慢。
根据相关评论编辑
- 即使使用整数指数,使用
pow
也可能会产生不正确的结果 (PaulMcKenzie)
- 除了使用 double 类型的数学函数外,
pow
是一个函数调用(而 x*x
不是)(jtbandes)
- 许多现代编译器实际上会优化带有常量整数参数的 pow,但不应依赖这一点。
您选择了一种最慢的检查方法
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执行单元饱和,它仍然是4倍左右比我让 gcc 发出的慢(下)。
您还可以优化循环条件以在条件的 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 个融合域微指令,。
Haswell 运行s imul r32, r32
具有每个时钟一次迭代的吞吐量,因此内部循环内的两个乘法不会以每 2.5c 的两个乘法使端口 1 饱和。这为吸收 ADD 和 SUB 窃取端口 1 不可避免的资源冲突留下了空间。
我们甚至没有接近任何其他执行端口瓶颈,因此 前端瓶颈是唯一的问题,这应该 运行 每 2.5 个周期进行一次迭代 在 Intel Haswell 及更高版本上。
循环展开可以帮助减少每次检查的微指令数。例如使用 lea ecx, [rbx+1]
计算下一次迭代的 b+1,因此我们可以 imul ebx, ebx
不使用 MOV 使其无损。
强度降低也是可能的:给定 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]
来更新 b*b
在 ECX 中,给定 -b
在 EBX 中)。
除非你真的在 IMUL 吞吐量上遇到瓶颈,否则这可能最终会花费更多的 uops 而不是胜利。看看编译器在 C++ 源代码中如何处理这种强度降低可能会很有趣。
您也可以使用 SSE 或 AVX 对其进行矢量化,并行检查 4 或 8 个连续的 b
值。由于命中非常罕见,您只需检查 8 个中是否有任何一个命中,然后在极少数情况下找出匹配的那个。
另请参阅 x86 标签 wiki 了解更多优化内容。
我一直在做一些项目 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)
实现细节我遗漏了,这使得第一个表达式的计算速度如此之慢?
pow()
适用于实浮点数,并在幕后使用公式
pow(x,y) = e^(y log(x))
计算x^y
。在调用 pow
之前,int
被转换为 double
。 (log
为自然对数,以e为基础)
x^2
使用 pow()
因此比 x*x
慢。
根据相关评论编辑
- 即使使用整数指数,使用
pow
也可能会产生不正确的结果 (PaulMcKenzie) - 除了使用 double 类型的数学函数外,
pow
是一个函数调用(而x*x
不是)(jtbandes) - 许多现代编译器实际上会优化带有常量整数参数的 pow,但不应依赖这一点。
您选择了一种最慢的检查方法
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执行单元饱和,它仍然是4倍左右比我让 gcc 发出的慢(下)。
您还可以优化循环条件以在条件的 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 个融合域微指令,
Haswell 运行s imul r32, r32
具有每个时钟一次迭代的吞吐量,因此内部循环内的两个乘法不会以每 2.5c 的两个乘法使端口 1 饱和。这为吸收 ADD 和 SUB 窃取端口 1 不可避免的资源冲突留下了空间。
我们甚至没有接近任何其他执行端口瓶颈,因此 前端瓶颈是唯一的问题,这应该 运行 每 2.5 个周期进行一次迭代 在 Intel Haswell 及更高版本上。
循环展开可以帮助减少每次检查的微指令数。例如使用 lea ecx, [rbx+1]
计算下一次迭代的 b+1,因此我们可以 imul ebx, ebx
不使用 MOV 使其无损。
强度降低也是可能的:给定 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]
来更新 b*b
在 ECX 中,给定 -b
在 EBX 中)。
除非你真的在 IMUL 吞吐量上遇到瓶颈,否则这可能最终会花费更多的 uops 而不是胜利。看看编译器在 C++ 源代码中如何处理这种强度降低可能会很有趣。
您也可以使用 SSE 或 AVX 对其进行矢量化,并行检查 4 或 8 个连续的 b
值。由于命中非常罕见,您只需检查 8 个中是否有任何一个命中,然后在极少数情况下找出匹配的那个。
另请参阅 x86 标签 wiki 了解更多优化内容。