为什么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).)

我很好奇的是,当我更换pow(a,6)与a*a*a*a*a*a使用GCC 4.5.1和选项" -O3 -lm -funroll-loops -msse4",它采用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

因为浮点数学不是关联的.在浮点乘法中对操作数进行分组的方式会影响答案的数值准确性.

因此,大多数编译器对浮点计算的重新排序非常保守,除非他们能够确定答案保持不变,或者除非你告诉他们你不关心数值精度.例如:该-fassociative-math选项的gcc允许GCC重新关联浮点运算,或甚至-ffast-math其允许甚至对速度精度的更积极的折衷选择.

Answer 2

Lambdageek正确地指出,因为关联性不适用于浮点数,所以a*a*a*a*a*ato的"优化"(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++没有power运算符,只有库函数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).

这样可以很容易地输入功率(不需要a用parens 拼出6 秒),并且允许你进行这种优化,而不-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

当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 -x c -

有很多gcc标志,但没有什么花哨的.他们的意思是:从stdin读取; 使用O2优化级别; 输出汇编语言列表而不是二进制; 列表应使用英特尔汇编语言语法; 输入是用C语言编写的(通常是从输入文件扩展名推断语言,但从stdin读取时没有文件扩展名); 并写信给stdout.

这是输出的重要部分.我用一些评论来注释它,表明汇编语言中发生了什么:

; 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"是可以的.他们没有看到"浮动a = 1.2","a*a*a"和类似的C代码中隐含的"e".

如果大多数程序员认识到(并且能够执行)C表达式a*a*a*a*a*a实际上并不适用于理想数字的想法,那么GCC编译器将可以自由地优化"a*a*a*a*a*a"to say"t =(a*a); t*t*t",它需要较少的乘法次数.但不幸的是,GCC编译器不知道编写代码的程序员是否认为"a"是带有或不带错误的数字.所以GCC只会做源代码的样子 - 因为这就是GCC用"肉眼"看到的东西.

......一旦你知道什么样的程序员的你是,你可以使用"-ffast -数学"开关告诉GCC说:"嘿,GCC,我知道我在做什么!".这将允许GCC将*a*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_CONTRACTpragma设置为ON,则允许编译器将表达式a*a*a*a*a*a视为单个操作,就像使用单个舍入精确计算一样.例如,编译器可以用更快和更准确的内部功率函数代替它.这一点特别有趣,因为行为部分由程序员直接在源代码中控制,而最终用户提供的编译器选项有时可能会被错误地使用.

FP_CONTRACTpragma 的默认状态是实现定义的,因此默认情况下允许编译器执行此类优化.因此,需要严格遵循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,因此遵循上段).在过去,后一种选择并没有阻止转型,这意味着海湾合作委员会不符合这一点:https://gcc.gnu.org/bugzilla/show_bug.cgi？id = 37845

Answer 8

我不希望这个案例得到优化.在表达式包含可以重新分组以删除整个操作的子表达式的情况下,通常不会这样.我希望编译器编写者将时间投入到更有可能带来明显改进的领域,而不是覆盖很少遇到的边缘情况.

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

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

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

Answer 9

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

Answer 10

像"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的选择,不论是任何特定的优化或编译器标志)要优化应该用"a*a"替换"pow(a,2)".这将是编译器供应商应该做的唯一理智的事情.

Answer 11

这个问题已经有了一些很好的答案,但为了完整起见,我想指出C标准的适用部分是5.1.2.2.3/15(与第1.9/9节中的相同) C++ 11标准).本节规定,如果运算符实际上是关联的或可交换的,则只能重新分组.

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.对于无符号整数,它更容易,因为它们的mod功率为2,因此即使面对溢出也可以自由重新排序.