为什么GCC不优化a*a*a*a*a*a到(a*a*a)*(a*a*a)?

我正在对一个科学应用程序进行一些数值优化。我注意到的一件事是,GCC将通过编译a*a来优化调用pow(a,2),但调用pow(a,6)没有优化,实际上会调用库函数pow,这大大降低了性能。(相比之下,英特尔C++编译器,可执行文件icc,将消除pow(a,6)的库调用。)

我好奇的是,当我使用GCC 4.5.1和选项“-O3 -lm -funroll-loops -msse4”将pow(a,6)替换为a*a*a*a*a*a时,它使用了5mulsd指令:

movapd  %xmm14, %xmm13mulsd   %xmm14, %xmm13mulsd   %xmm14, %xmm13mulsd   %xmm14, %xmm13mulsd   %xmm14, %xmm13mulsd   %xmm14, %xmm13

如果我写(a*a*a)*(a*a*a),它会产生

movapd  %xmm14, %xmm13mulsd   %xmm14, %xmm13mulsd   %xmm14, %xmm13mulsd   %xmm13, %xmm13

这将乘法指令的数量减少到3。icc也有类似的行为。

为什么编译器不能识别这种优化技巧?

227589 次浏览

我根本不希望这种情况得到优化。表达式包含可以重新组合以删除整个操作的子表达式的情况不太常见。我希望编译器作者将时间投入到更有可能导致显着改进的领域,而不是覆盖很少遇到的边缘情况。

我很惊讶地从其他答案中得知,这个表达式确实可以通过适当的编译器开关进行优化。要么优化是微不足道的,要么是更常见优化的边缘情况,要么编译器作者非常彻底。

像你在这里所做的那样向编译器提供提示没有错。重新排列语句和表达式以查看它们会带来什么差异是微优化过程中正常且预期的部分。

虽然编译器可以合理地考虑这两个表达式以提供不一致的结果(没有适当的开关),但您无需受该限制的约束。差异将非常小-以至于如果差异对您很重要,您不应该首先使用标准浮点运算。

因为浮点数不是关联的。浮点乘法中对操作数进行分组的方式会影响答案的数值精度。

因此,大多数编译器对重新排序浮点计算非常保守,除非他们可以确保答案保持不变,或者除非你告诉他们你不关心数值精度。例如:gcc的#0选项允许gcc重新关联浮点运算,甚至-ffast-math选项允许更积极的精度与速度权衡。

Lambda极客正确地指出,因为结合性不适用于浮点数,a*a*a*a*a*a(a*a*a)*(a*a*a)的“优化”可能会改变值。这就是为什么C99不允许它(除非用户特别允许,通过编译器标志或Pragma)。通常,假设程序员写她所做的是有原因的,编译器应该尊重这一点。如果你想要(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-08worst relative error using (a*a*a)*(a*a*a): 2.94e-07worst relative error using     a*a*a*a*a*a: 2.58e-07

使用pow而不是乘法树可以减少因子4绑定的错误。编译器不应该(通常也不会)进行增加错误的“优化”,除非用户允许这样做(例如通过-ffast-math)。

请注意,GCC提供__builtin_powi(x,n)作为pow( )的替代方案,它应该生成一个内联乘法树。如果您想权衡准确性以获得性能,但不想启用快速数学,请使用它。

另一个类似的情况:大多数编译器不会优化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

因为32位浮点数——比如1.024——不是1.024。在计算机中,1.024是一个区间:从(1.024-e)到(1.024+e),其中“e”代表错误。有些人没有意识到这一点,还认为*在*a中代表任意精度数字的乘法,这些数字没有任何错误。有些人没有意识到这一点的原因可能是他们在小学时练习的数学计算:只处理没有错误的理想数字,并认为在进行乘法时简单地忽略“e”是可以的。他们没有看到“浮点数a=1.2”、“a*a*a”和类似的C代码中隐含的“e”。

如果大多数程序员认识到(并且能够执行)C表达式a*a*a*a*a*a*a实际上不适用于理想数字的想法,那么GCC编译器可以自由优化“a*a*a*a*a*a”成“t=(a*a); t*t*t”,这需要较少的乘法次数。但不幸的是,GCC编译器不知道编写代码的程序员是否认为“a”是一个有错误的数字还是没有错误的数字。因此GCC只会做源代码的样子——因为这是GCC用它的“肉眼”看到的。

…一旦你知道了是什么样的程序员,你就可以使用“-ffast-max”开关告诉GCC“嘿,GCC,我知道我在做什么!”。这将允许GCC将a*a*a*a*a*a转换为不同的文本片段——它看起来与a*a*a*a*a*a不同——但仍然在a*a*a*a*a*a*a的错误区间内计算一个数字。这没问题,因为你已经知道你正在使用的是区间,而不是理想数字。

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

为好奇的人澄清:这并没有找到计算幂的最佳方法,但是因为找到最佳解是一个NP完全问题而且这只值得为小幂做(而不是使用pow),没有理由对细节大惊小怪。

然后将其用作power<6>(a)

这使得键入幂变得容易(不需要用括号拼出6a),并且可以让您在没有-ffast-math的情况下进行这种优化,以防您有精度依赖的东西,例如补偿求和(操作顺序至关重要的示例)。

您可能还会忘记这是C++,只需在C程序中使用它(如果它使用C++编译器编译)。

希望这能有用。

编辑:

这是我从我的编译器中得到的:

对于a*a*a*a*a*a

    movapd  %xmm1, %xmm0mulsd   %xmm1, %xmm0mulsd   %xmm1, %xmm0mulsd   %xmm1, %xmm0mulsd   %xmm1, %xmm0mulsd   %xmm1, %xmm0

对于(a*a*a)*(a*a*a)

    movapd  %xmm1, %xmm0mulsd   %xmm1, %xmm0mulsd   %xmm1, %xmm0mulsd   %xmm0, %xmm0

对于power<6>(a)

    mulsd   %xmm0, %xmm0movapd  %xmm0, %xmm1mulsd   %xmm0, %xmm1mulsd   %xmm0, %xmm1

正如Lambdageek所指出的,浮点乘法不是关联的,你可以获得更低的精度,但是当获得更好的精度时,你可以反对优化,因为你想要一个确定性的应用程序。例如在游戏模拟客户端/服务器中,每个客户端都必须模拟同一个世界,你希望浮点计算是确定性的。

这个问题已经有了一些很好的答案,但为了完整起见,我想指出C标准的适用部分是5.1.2.2.3/15(与C++11标准中的第1.9/9节相同)。本节指出,只有当操作符真正是关联或交换的时,才能重新组合。

当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标志,但没有什么花哨的。它们的意思是:从stdin读取;使用O2优化级别;输出汇编语言列表而不是二进制;列表应使用英特尔汇编语言语法;输入是C语言(通常从输入文件扩展名推断语言,但从stdin读取时没有文件扩展名);并写入标准输出。

这是输出的重要部分。我用一些注释对其进行了注释,说明了汇编语言中发生的事情:

; x is in edi to begin with.  eax will be used as a temporary register.mov  eax, edi  ; temp = ximul eax, edi  ; temp = x * tempimul eax, edi  ; temp = x * tempimul eax, eax  ; temp = temp * temp

我在Ubuntu衍生LinuxMint 16 Petra上使用系统GCC。这是gcc版本:

$ gcc --versiongcc (Ubuntu/Linaro 4.8.1-10ubuntu9) 4.8.1

正如其他海报所指出的,这个选项在浮点数中是不可能的,因为浮点数算术不是关联的。

目前还没有海报提到浮动表达式的收缩(ISO C标准,6.5p8和7.12.2)。如果FP_CONTRACT pragma设置为ON,则允许编译器将a*a*a*a*a*a等表达式视为单个操作,就好像使用单个舍入精确计算一样。例如,编译器可能会用更快更准确的内部幂函数替换它。这特别有趣,因为这种行为部分由程序员直接在源代码中控制,而最终用户提供的编译器选项有时可能会被错误使用。

FP_CONTRACT Pragma的默认状态是实现定义的,因此默认情况下允许编译器进行此类优化。因此,需要严格遵循IEEE 754规则的可移植代码应显式将其设置为OFF

如果编译器不支持此语用,则必须通过避免任何此类优化来保持保守,以防开发人员选择将其设置为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个乘法)。

编译器应该非常小心它正在进行的优化类型:

  1. 如果优化pow(a,6)a*a*a*a*a*a,则可能可以提高性能,但会大大降低浮点数的准确性。
  2. 如果优化a*a*a*a*a*apow(a,6),它实际上可能会降低精度,因为“a”是一些允许无错误乘法的特殊值(2的幂或一些小整数)
  3. 如果优化pow(a,6)(a*a*a)*(a*a*a)(a*a)*(a*a)*(a*a),与pow函数相比,仍然可能会丢失精度。

一般来说,您知道对于任意浮点值“pow”比您最终可以编写的任何函数都具有更好的准确性,但在某些特殊情况下,多次乘法可能具有更好的准确性和性能,这取决于开发人员选择更合适的内容,最终注释代码,以便其他人不会“优化”该代码。

唯一有意义的优化(个人意见,显然是GCC中没有任何特定优化或编译器标志的选择)应该将“pow(a,2)”替换为“a*a”。这是编译器供应商唯一应该做的理智的事情。

GCC实际上可以进行这种优化,即使对于浮点数也是如此。例如,

double foo(double a) {return a*a*a*a*a*a;}

成为

foo(double):mulsd   %xmm0, %xmm0movapd  %xmm0, %xmm1mulsd   %xmm0, %xmm1mulsd   %xmm1, %xmm0ret

但是,这种重新排序违反了IEEE-754,因此它需要标志。

正如Peter Cordes在评论中指出的那样,有符号整数可以在没有-funsafe-math-optimizations的情况下进行这种优化,因为它在没有溢出的时候准确地保持不变,如果有溢出,你会得到未定义的行为。所以你得到

foo(long):movq    %rdi, %raximulq   %rdi, %raximulq   %rdi, %raximulq   %rax, %raxret

对于无符号整数,它甚至更容易,因为它们的mod幂为2,因此即使在溢出的情况下也可以自由地重新排序。