为什么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具有相似的行为。

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

Answer 1

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

因此，大多数编译器对重新计算浮点计算都非常保守，除非他们能确定答案保持不变，或者除非您告诉他们您不关心数值精度。 例如：gcc 的-fassociative-math选项允许gcc重新关联浮点运算，甚至是-ffast-math选项，它允许更准确地反对速度。

Answer 2

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( )的替代，它应该生成一个内联乘法树。 如果你想要牺牲性能的准确性，但不想启用快速数学，那么使用它。

Answer 3

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

Answer 4

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 

Answer 5

因为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。 这是可以的，因为你已经知道你正在使用间隔，而不是理想的数字。

Answer 6

当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 

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

Answer 7

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

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

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

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

Answer 8

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

Answer 9

没有海报提到浮动表达式的收缩（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

Answer 10

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

Answer 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）”。 这将是编译器厂商应该做的唯一理智的事情。

Answer 12

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权，所以即使面对溢出，也可以自由地重新排序。