为什么GCC不能优化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)的库调用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具有相似的行为。

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

因为浮点数学不是联想的 。 用浮点乘法对操作数进行分组的方式对答案的数字精度有影响。

因此,大多数编译器对重新计算浮点计算都非常保守,除非他们能确定答案保持不变,或者除非您告诉他们您不关心数值精度。 例如:gcc 的-fassociative-math选项允许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]之间所有单精度浮点数的^ 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( )的替代,它应该生成一个内联乘法树。 如果你想要牺牲性能的准确性,但不想启用快速数学,那么使用它。

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

澄清好奇:这并没有找到最佳的方式来计算权力,但因为找到最佳的解决方案是一个NP完全问题 ,这只是值得做的小功率(而不是使用pow ),没有理由与细节大惊小怪。

那么就把它当作power<6>(a)

这样可以很容易地输入权力(不需要用parens拼出6 a s),并且可以在不需要-ffast-math情况下进行这种优化,例如补偿求和 (例子的操作是必不可少的)。

你也许可以忘记这是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 

因为32位浮点数(例如1.024)不是1.024。 在计算机中,1.024是一个时间间隔:从(1.024-e)到(1.024 + e),其中“e”代表错误。 有些人没有意识到这一点,也相信* 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”是一个有或没有错误的数字。 所以海湾合作委员会只会做源代码的样子 – 因为这是海湾合作委员会用“裸眼”所看到的。

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

当a是一个整数时,GCC实际上优化了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 -xc - 

有很多海湾合作委员会的标志,但没有什么幻想。 他们的意思是:从标准输入读取; 使用O2优化级别; 输出汇编语言列表而不是二进制文件; 该列表应使用英特尔汇编语言语法; 输入是C语言(通常语言是从输入文件扩展名推断出来的,但从stdin读取时没有文件扩展名); 并写入标准输出。

这是输出的重要部分。 我用注释说明了汇编语言中发生了什么:

  ; x is in edi to begin with. eax will be used as a temporary register. mov eax, edi ; temp1 = x imul eax, edi ; temp2 = x * temp1 imul eax, edi ; temp3 = x * temp2 imul eax, eax ; temp4 = temp3 * temp3 

我在Linux Mint 16 Petra上使用了系统GCC,这是一个Ubuntu衍生产品。 这是gcc版本:

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

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

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

我很惊讶地从其他答案中学习,这个表达式确实可以用适当的编译器开关来优化。 要么是优化是微不足道的,要么是更为普遍的优化的边缘情况,或者编译器作者是非常彻底的。

正如你在这里所做的那样,向编译器提供提示没有任何问题。 微观优化过程中的正常和期望的部分是重新排列语句和表达式,以查看它们将带来的差异。

虽然编译器在考虑这两个表达式来提供不一致的结果(没有正确的开关)时可能是合理的,但是不需要受限制。 这个差别是非常小的,所以如果差别对你有影响,你不应该首先使用标准的浮点运算。

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

没有海报提到浮动表达式的收缩(ISO C标准,6.5p8和7.12.2)。 如果FP_CONTRACT杂注被设置为“on”,则编译器允许将诸如a a a * a之类的表达看作是单个操作,就好像通过单个舍入来精确评估。 例如,编译器可以用更快更准确的内部电源功能代替它。 这是特别有趣的,因为行为是由程序员直接在源代码中部分控制的,而由最终用户提供的编译器选项有时可能被错误地使用。

FP_CONTRACT编译指示的默认状态是实现定义的,所以默认情况下允许编译器进行这种优化。 因此,需要严格遵循IEEE 754规则的便携式代码应明确地将其设置为“关闭”。

如果编译器不支持这个编译指示,那么在开发人员选择将其设置为“off”的情况下,它必须避免任何这样的优化。

海湾合作委员会不支持这个编译指示。 但是,对于具有硬件FMA的目标,它仍然可以(有时是无效的)转换a * b + c到FMA(a,b,c): https ://gcc.gnu.org/bugzilla/show_bug.cgi ? id =37845

这个问题已经有了一些很好的答案,但是为了完整起见,我想指出C标准的适用部分是5.1.2.2.3 / 15(与第1.9 / 9节中的相同C ++ 11标准)。 本节规定,如果运营商是真正的联想性或交换性的,则只能进行重组。

像“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. 如果对pow(a,6)优化a*a*a*a*a*a pow(a,6)实际上可能会降低准确性,因为“a”是允许乘法没有错误的一些特殊值(2的幂或某个小整数)
  3. 如果优化pow(a,6)(a*a*a)*(a*a*a)(a*a)*(a*a)*(a*a) ,仍然会有精度损失相比pow功能。

一般来说,对于任意浮点值,“pow”比任何你可能写的函数都有更好的准确性,但是在一些特殊情况下,多次乘法可能会有更好的精度和性能,这取决于开发者选择什么更合适,最终对代码进行评论,以便其他人不会“优化”该代码。

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

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 optimizations的情况下做这种优化,因为它恰好在没有溢出的情况下持有,并且如果有溢出,你会得到未定义的行为。 所以你明白了

 foo(long): movq %rdi, %rax imulq %rdi, %rax imulq %rdi, %rax imulq %rax, %rax ret 

只有-O 。 对于无符号整数来说,它更容易,因为它们工作在2的mod权,所以即使面对溢出,也可以自由地重新排序。

Interesting Posts