我正在对科学应用程序进行一些数值优化。我注意到的一件事是,GCC 会通过将调用 pow(a,2)
编译成 a*a
来优化它,但调用 pow(a,6)
并没有优化,实际上会调用库函数 pow
,这会大大降低性能。 (相比之下,可执行文件 icc
的 Intel C++ Compiler 将消除对 pow(a,6)
的库调用。)
我很好奇的是,当我使用 GCC 4.5.1 和选项“-O3 -lm -funroll-loops -msse4
”将 pow(a,6)
替换为 a*a*a*a*a*a
时,它使用了 5 个 mulsd
指令:
movapd %xmm14, %xmm13
mulsd %xmm14, %xmm13
mulsd %xmm14, %xmm13
mulsd %xmm14, %xmm13
mulsd %xmm14, %xmm13
mulsd %xmm14, %xmm13
而如果我写(a*a*a)*(a*a*a)
,它会产生
movapd %xmm14, %xmm13
mulsd %xmm14, %xmm13
mulsd %xmm14, %xmm13
mulsd %xmm13, %xmm13
这将乘法指令的数量减少到 3。icc
具有类似的行为。
为什么编译器无法识别这种优化技巧?
(a*a)*(a*a)*(a*a)
加入其中。相同数量的乘法,但可能更准确。
因为Floating Point Math is not Associative。在浮点乘法中对操作数进行分组的方式会影响答案的数值准确性。
因此,大多数编译器对重新排序浮点计算非常保守,除非他们可以确定答案将保持不变,或者除非您告诉他们您不关心数值准确性。例如:gcc 的 the -fassociative-math
option 允许 gcc 重新关联浮点操作,甚至 -ffast-math
选项允许在精度与速度之间进行更激进的权衡。
Lambdageek 正确指出,因为关联性不适用于浮点数,a*a*a*a*a*a
到 (a*a*a)*(a*a*a)
的“优化”可能会改变值。这就是 C99 不允许它的原因(除非用户特别允许,通过编译器标志或编译指示)。通常,假设是程序员编写了她所做的事情是有原因的,编译器应该尊重这一点。如果您想要 (a*a*a)*(a*a*a)
,请写下它。
不过,写起来可能会很痛苦。当您使用 pow(a,6)
时,为什么编译器不能 [您认为是] 正确的事情?因为这样做是错误的事情。在具有良好数学库的平台上,pow(a,6)
比 a*a*a*a*a*a
或 (a*a*a)*(a*a*a)
准确得多。只是为了提供一些数据,我在我的 Mac Pro 上进行了一个小实验,测量 [1,2) 之间的所有单精度浮点数评估 a^6 时的最严重错误:
worst relative error using powf(a, 6.f): 5.96e-08
worst relative error using (a*a*a)*(a*a*a): 2.94e-07
worst relative error using a*a*a*a*a*a: 2.58e-07
使用 pow
而不是乘法树可以减少 4 倍的误差。编译器不应该(通常也不应该)进行会增加错误的“优化”,除非得到用户的许可(例如通过 -ffast-math
)。
请注意,GCC 提供 __builtin_powi(x,n)
作为 pow( )
的替代,它应该生成一个内联乘法树。如果您想在性能上权衡准确性,但又不想启用快速数学,请使用它。
flag=1
调用 _set_SSE2_enable(<flag>)
,它将尽可能使用 SSE2。这会稍微降低准确性,但会提高速度(在某些情况下)。 MSDN:_set_SSE2_enable() 和 pow()
pow
。有一些基于 SSE 的 pow
实现比大多数基于 x87 的实现更准确,还有一些实现会牺牲一些准确性来换取速度。
a*a*a*a*a*a
,但显然情况并非如此! :)
另一个类似的情况:大多数编译器不会将 a + b + c + d
优化为 (a + b) + (c + d)
(这是一种优化,因为可以更好地流水线化第二个表达式)并将其评估为给定(即作为 (((a + b) + c) + d)
)。这也是因为极端情况:
float a = 1e35, b = 1e-5, c = -1e35, d = 1e-5;
printf("%e %e\n", a + b + c + d, (a + b) + (c + d));
这会输出 1.000000e-05 0.000000e+00
Fortran(专为科学计算而设计)具有内置的幂运算符,据我所知,Fortran 编译器通常会以与您描述的方式类似的方式优化整数幂。遗憾的是,C/C++ 没有幂运算符,只有库函数 pow()
。这并不妨碍智能编译器特别处理 pow
并在特殊情况下以更快的方式计算它,但似乎他们不太常见......
几年前,我试图让以最佳方式计算整数幂更方便,并提出了以下建议。它是 C++,而不是 C,并且仍然取决于编译器在如何优化/内联事物方面有点聪明。无论如何,希望你会发现它在实践中很有用:
template<unsigned N> struct power_impl;
template<unsigned N> struct power_impl {
template<typename T>
static T calc(const T &x) {
if (N%2 == 0)
return power_impl<N/2>::calc(x*x);
else if (N%3 == 0)
return power_impl<N/3>::calc(x*x*x);
return power_impl<N-1>::calc(x)*x;
}
};
template<> struct power_impl<0> {
template<typename T>
static T calc(const T &) { return 1; }
};
template<unsigned N, typename T>
inline T power(const T &x) {
return power_impl<N>::calc(x);
}
对好奇的澄清:这并没有找到计算幂的最佳方法,但是因为 finding the optimal solution is an NP-complete problem 并且无论如何这只值得为小幂做(而不是使用 pow
) ,没有理由对细节大惊小怪。
然后只需将其用作 power<6>(a)
。
这使得键入幂变得容易(无需用括号拼出 6 个 a
),并允许您在没有 -ffast-math
的情况下进行这种优化,以防您有一些依赖于精度的东西,例如 compensated summation(一个示例,其中操作顺序是必不可少的)。
您可能还忘记了这是 C++,而只是在 C 程序中使用它(如果它使用 C++ 编译器编译)。
希望这会很有用。
编辑:
这是我从编译器中得到的:
对于 a*a*a*a*a*a
,
movapd %xmm1, %xmm0
mulsd %xmm1, %xmm0
mulsd %xmm1, %xmm0
mulsd %xmm1, %xmm0
mulsd %xmm1, %xmm0
mulsd %xmm1, %xmm0
对于 (a*a*a)*(a*a*a)
,
movapd %xmm1, %xmm0
mulsd %xmm1, %xmm0
mulsd %xmm1, %xmm0
mulsd %xmm0, %xmm0
对于 power<6>(a)
,
mulsd %xmm0, %xmm0
movapd %xmm0, %xmm1
mulsd %xmm0, %xmm1
mulsd %xmm0, %xmm1
当 a 为整数时,GCC 确实将 a*a*a*a*a*a
优化为 (a*a*a)*(a*a*a)
。我试过这个命令:
$ echo 'int f(int x) { return x*x*x*x*x*x; }' | gcc -o - -O2 -S -masm=intel -x c -
有很多 gcc 标志,但没有什么花哨的。它们的意思是:从标准输入读取;使用 O2 优化级别;输出汇编语言列表而不是二进制;清单应使用英特尔汇编语言语法;输入是 C 语言(通常语言是从输入文件扩展名推断出来的,但是从标准输入读取时没有文件扩展名);并写入标准输出。
这是输出的重要部分。我用一些注释来注释它,说明汇编语言中发生了什么:
; x is in edi to begin with. eax will be used as a temporary register.
mov eax, edi ; temp = x
imul eax, edi ; temp = x * temp
imul eax, edi ; temp = x * temp
imul eax, eax ; temp = temp * temp
我在 Linux Mint 16 Petra(Ubuntu 衍生产品)上使用系统 GCC。这是 gcc 版本:
$ gcc --version
gcc (Ubuntu/Linaro 4.8.1-10ubuntu9) 4.8.1
正如其他海报所指出的,此选项在浮点中是不可能的,因为浮点运算不是关联的。
unsigned int
进行了此操作。
a*a*a*a*a*a
和 (a*a*a)*(a*a*a)
得到相同的结果。 (当然对于无符号类型,溢出不是UB。)
因为 32 位浮点数(例如 1.024)不是 1.024。在计算机中,1.024 是一个区间:从 (1.024-e) 到 (1.024+e),其中“e”表示错误。有些人没有意识到这一点,并且还认为 * in a*a 代表任意精度数字的乘法,而这些数字没有任何错误。有些人没有意识到这一点的原因可能是他们在小学时进行的数学计算:只使用没有附加错误的理想数字,并认为在执行乘法时可以简单地忽略“e”。他们看不到“float a=1.2”、“a*a*a”和类似的 C 代码中隐含的“e”。
如果大多数程序员认识到(并能够执行)C 表达式 a*a*a*a*a*a 实际上不是使用理想数字的想法,那么 GCC 编译器将可以免费优化 "a*a *a*a*a*a" 变成说 "t=(a*a); t*t*t" 这需要更少的乘法。但不幸的是,GCC 编译器不知道编写代码的程序员是否认为“a”是一个有错误或没有错误的数字。所以 GCC 只会做源代码的样子——因为那是 GCC 用它的“肉眼”看到的。
...一旦您知道自己是哪种程序员,就可以使用“-ffast-math”开关告诉GCC“嘿,GCC,我知道我在做什么!”。这将允许 GCC 将 a*a*a*a*a*a 转换为不同的文本 - 它看起来与 a*a*a*a*a*a 不同 - 但仍会在错误间隔内计算一个数字a*a*a*a*a*a。这没关系,因为您已经知道您正在使用间隔,而不是理想的数字。
int x = 3
解释为意味着 x
是 3+/-0.5。
Distance
不完全等于它的数值;这意味着数值只是对某些正在建模的物理量的近似值。
还没有海报提到浮动表达式的收缩(ISO C 标准,6.5p8 和 7.12.2)。如果 FP_CONTRACT
pragma 设置为 ON
,则允许编译器将诸如 a*a*a*a*a*a
的表达式视为单个操作,就像使用单个舍入精确计算一样。例如,编译器可以用更快、更准确的内部幂函数代替它。这一点特别有趣,因为行为部分是由程序员直接在源代码中控制的,而最终用户提供的编译器选项有时可能会被错误地使用。
FP_CONTRACT
pragma 的默认状态是实现定义的,因此默认情况下允许编译器进行此类优化。因此,需要严格遵循 IEEE 754 规则的可移植代码应将其显式设置为 OFF
。
如果编译器不支持此 pragma,则必须避免任何此类优化,以防开发人员选择将其设置为 OFF
。
GCC 不支持这个 pragma,但是使用默认选项,它假定它是 ON
;因此对于具有硬件 FMA 的目标,如果想要阻止将 a*b+c
转换为 fma(a,b,c),则需要提供诸如 -ffp-contract=off
之类的选项(将 pragma 显式设置为 OFF
)或 -std=c99
(告诉 GCC 符合某些 C 标准版本,这里是 C99,因此遵循上述段落)。过去,后一个选项不会阻止转换,这意味着 GCC 不符合这一点:https://gcc.gnu.org/bugzilla/show_bug.cgi?id=37845
像“pow”这样的库函数通常经过精心设计,以产生最小的可能错误(在一般情况下)。这通常是通过样条曲线逼近函数来实现的(根据 Pascal 的评论,最常见的实现似乎是使用 Remez algorithm)
基本上是以下操作:
pow(x,y);
具有与任何单次乘法或除法中的误差大致相同大小的固有误差。
同时进行以下操作:
float a=someValue;
float b=a*a*a*a*a*a;
具有大于单个乘法或除法误差的 5 倍以上的固有误差(因为您正在组合 5 次乘法)。
编译器应该非常小心它正在执行的优化类型:
如果将 pow(a,6) 优化为 a*a*a*a*a*a 它可能会提高性能,但会大大降低浮点数的准确性。如果将 a*a*a*a*a*a 优化为 pow(a,6) 它实际上可能会降低精度,因为“a”是一些允许乘法无误的特殊值(2 的幂或一些小整数)如果将 pow(a,6) 优化为 (a*a*a)*(a*a*a) 或 (a*a)*(a*a)*(a*a) 仍然可能会损失准确性与 pow 函数相比。
一般来说,您知道对于任意浮点值,“pow”比您最终编写的任何函数具有更好的精度,但在某些特殊情况下,多次乘法可能具有更好的精度和性能,这取决于开发人员选择更合适的方法,最终评论代码,以便其他人不会“优化”该代码。
唯一有意义的事情(个人意见,显然是 GCC 中的选择,没有任何特定的优化或编译器标志)要优化应该用“a*a”替换“pow(a,2)”。那将是编译器供应商应该做的唯一明智的事情。
正如 Lambdageek 指出的那样,浮点乘法不是关联的,您可以获得较低的准确性,但当获得更高的准确性时,您可以反对优化,因为您想要一个确定性的应用程序。例如,在游戏模拟客户端/服务器中,每个客户端都必须模拟同一个世界,您希望浮点计算具有确定性。
我根本没想到这种情况会被优化。表达式包含可以重新组合以删除整个操作的子表达式的情况并不常见。我希望编译器编写者将他们的时间投入到更有可能带来显着改进的领域,而不是涵盖很少遇到的边缘情况。
我很惊讶地从其他答案中得知,这个表达式确实可以通过适当的编译器开关进行优化。要么优化是微不足道的,要么是更常见优化的边缘案例,或者编译器编写者非常彻底。
正如您在此处所做的那样,向编译器提供提示并没有错。重新排列语句和表达式以查看它们会带来什么差异是微优化过程中正常且预期的部分。
虽然编译器可能有理由考虑这两个表达式来提供不一致的结果(没有适当的开关),但您无需受该限制的约束。差异将非常小 - 以至于如果差异对您很重要,那么您首先不应该使用标准浮点运算。
这个问题已经有一些很好的答案,但为了完整起见,我想指出 C 标准的适用部分是 5.1.2.2.3/15(与第 1.9/9 节相同) C++11 标准)。本节指出,只有当它们真正具有关联性或可交换性时,才能重新组合运算符。
gcc 实际上可以进行这种优化,即使对于浮点数也是如此。例如,
double foo(double a) {
return a*a*a*a*a*a;
}
变成
foo(double):
mulsd %xmm0, %xmm0
movapd %xmm0, %xmm1
mulsd %xmm0, %xmm1
mulsd %xmm1, %xmm0
ret
与 -O -funsafe-math-optimizations
。但是,这种重新排序违反了 IEEE-754,因此它需要该标志。
正如 Peter Cordes 在评论中指出的那样,有符号整数可以在没有 -funsafe-math-optimizations
的情况下进行此优化,因为它恰好在没有溢出时保持不变,如果有溢出,您会得到未定义的行为。所以你得到
foo(long):
movq %rdi, %rax
imulq %rdi, %rax
imulq %rdi, %rax
imulq %rax, %rax
ret
只有 -O
。对于无符号整数,它甚至更容易,因为它们可以使用 2 的模幂,因此即使在溢出的情况下也可以自由重新排序。
-ffast-math
)
pow
的实现细节既不存在也不存在;这个答案甚至没有引用pow
。-fp-model precise
。clang
和gcc
默认为严格一致性 wrt 重新关联。-fassociative-math
不准确;只是a*a*a*a*a*a
和(a*a*a)*(a*a*a)
不同。这与准确性无关。它是关于标准一致性和严格可重复的结果,例如在任何编译器上的相同结果。浮点数已经不准确了。使用-fassociative-math
编译很少不合适。