我正在对一个科学应用程序进行数值优化。我注意到的一件事是,GCC将通过将调用pow(a,2)编译为a*a来优化它,但调用pov(a,6)并没有优化,实际上会调用库函数pow,这会大大降低性能。(相比之下,可执行icc的“英特尔C++编译器”将消除对pow(a,6)的库调用。)
我好奇的是,当我使用GCC 4.5.1和选项“-O3-lm-funroll-loops-msse4”将pow(a,6)替换为a*a*a*a*a*a时,它使用了5条多指令:
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也有类似的行为。
为什么编译器不认识这种优化技巧?
还没有海报提到浮动表达式的收缩(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算法)
基本上是以下操作:
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)*。
一般来说,您知道对于任意浮点值,“pow”的精度比您最终可以编写的任何函数都要高,但在某些特殊情况下,多次乘法可能具有更好的精度和性能,这取决于开发人员选择更合适的方法,最终对代码进行注释,以便其他人不会“优化”该代码。
唯一有意义的优化(个人观点,显然是GCC中没有任何特定优化或编译器标志的选择)应该是将“pow(a,2)”替换为“a*a”。这将是编译器供应商应该做的唯一明智的事情。
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-funcafe数学优化。但是,这种重新排序违反了IEEE-754,因此需要标记。
正如Peter Cordes在一篇评论中指出的,有符号整数可以在没有funsafe数学优化的情况下进行这种优化,因为它恰好在没有溢出的情况下有效,如果有溢出,则会出现未定义的行为。所以你得到
foo(long):
movq %rdi, %rax
imulq %rdi, %rax
imulq %rdi, %rax
imulq %rax, %rax
ret
只需-O。对于无符号整数,这更容易,因为它们是2的模幂,因此即使在溢出的情况下也可以自由地重新排序。