为什么 GCC 不将 a*a*a*a*a*a 优化为 (a*a*a)*(a*a*a)？

Question

我正在对科学应用程序进行一些数值优化。我注意到的一件事是GCC会通过将调用pow(a,2)编译成a*a来优化调用，但是调用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时，它使用了 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] 之间所有单精度浮点数的 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( )的替代品，它应该生成一个内联乘法树。 如果您想在准确性和性能之间进行权衡，但又不想启用快速数学运算，请使用它。

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 %en", 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).

这使得键入幂变得容易(无需用参数拼写出 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

GCC实际上确实优化了a*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语言(通常语言是从输入文件扩展名推断出来的，但从stdin读取时没有文件扩展名(;并写信给标准输出。

这是输出的重要部分。 我用一些注释注释了它，指出汇编语言中发生了什么：

; 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上使用系统GCC，这是Ubuntu的衍生产品。 这是 gcc 版本：

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

正如其他海报所指出的，此选项在浮点运算中是不可能的，因为浮点运算不是关联的。

Answer 6

因为 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*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 的错误区间内计算一个数字。这没关系，因为您已经知道您正在使用间隔，而不是理想的数字。

Answer 7

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

FP_CONTRACT杂注的默认状态是实现定义的，因此默认情况下允许编译器执行此类优化。因此，需要严格遵守IEEE 754规则的可移植代码应明确设置为OFF。

如果编译器不支持此编译指示，则必须通过避免任何此类优化来保守，以防开发人员选择将其设置为 OFF 。

GCC 不支持此编译指示，但使用默认选项，它假定它是ON的;因此，对于具有硬件 FMA 的目标，如果要a*b+c防止转换为 fma(a，b，c(，则需要提供一个选项，例如 -ffp-contract=off(明确将编译指示设置为 OFF(或-std=c99(告诉 GCC 符合某些 C 标准版本， 此处为C99，因此遵循上一段(。过去，后一种选择不会阻止转型，这意味着海湾合作委员会在这一点上不符合规定：https://gcc.gnu.org/bugzilla/show_bug.cgi?id=37845

Answer 8

像"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*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中任何特定优化或编译器标志的选择(进行优化应该是将"pow(a，2("替换为"a*a"。这将是编译器供应商应该做的唯一理智的事情。

Answer 9

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

Answer 10

我根本没想到这种情况会得到优化。表达式包含可以重新分组以删除整个操作的子表达式的情况并不常见。我希望编译器编写者将时间投入到更有可能带来明显改进的领域，而不是涵盖很少遇到的边缘情况。

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

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

虽然编译器可能有理由认为这两个表达式提供不一致的结果(没有适当的开关(，但您没有必要受该限制的约束。差异将非常小 - 以至于如果差异对您很重要，那么您首先不应该使用标准的浮点运算。

Answer 11

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

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的情况下进行此优化，因为它恰好在没有溢出时成立，如果有溢出，您将获得未定义的行为。所以你得到

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

只有-O.对于无符号整数，这甚至更容易，因为它们的工作功率为 2，因此即使面对溢出也可以自由重新排序。