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


我正在对科学应用进行一些数值优化。我注意到的一件事是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 有类似的行为。

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


1965
2018-06-21 18:49


起源


“识别pow(a,6)”是什么意思? - Varun Madiath
嗯......你知道吗一个一个一个一个a和(a一个一)*(一a * a)与浮点数不一样,不是吗?你必须使用-funsafe-math或-ffast-math或其他东西。 - Damon
我建议你阅读David Goldberg撰写的“每个计算机科学家应该知道的关于浮点运算的知识”: download.oracle.com/docs/cd/E19957-01/806-3568/... 之后,您将对您刚刚走进的焦油坑有更全面的了解! - Phil Armstrong
一个非常合理的问题。 20年前,我问过同样的一般性问题,通过粉碎这个单一的瓶颈,将蒙特卡罗模拟的执行时间从21小时减少到7小时。内循环中的代码在此过程中执行了13万亿次,但它将模拟带入了一个过夜的窗口。 (见下面的答案)
也许扔 (a*a)*(a*a)*(a*a) 也参与其中。相同数量的乘法,但可能更准确。 - Rok Kralj


答案:


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

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


2567
2018-06-22 15:32



是。使用-ffast-math,它正在进行这样的优化。好主意!但由于我们的代码比速度更准确,因此最好不要通过它。 - xis
IIRC C99允许编译器进行这种“不安全”的FP优化,但是GCC(除了x87之外的任何东西)合理地尝试遵循IEEE 754 - 它不是“错误界限”; 只有一个正确的答案。 - tc.
实施细节 pow 既不在这里也不在那里;这个答案甚至没有参考 pow。 - Stephen Canon
@nedR:ICC默认允许重新关联。如果要获得符合标准的行为,则需要进行设置 -fp-model precise 与ICC。 clang 和 gcc 默认为严格一致性w.r.t.重新关联。 - Stephen Canon
@xis,实际上不是那样的 -fassociative-math 将是不准确的;就是这样 a*a*a*a*a*a 和 (a*a*a)*(a*a*a) 是不同的。这与准确性无关;它是关于标准一致性和严格可重复的结果,例如任何编译器的结果都相同。浮点数已经不准确了。编译时很少是不合适的 -fassociative-math。 - Paul Draper


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 而不是乘法树减少了由a约束的误差 因子4。编译器不应(并且通常不)进行“优化”以增加错误,除非用户许可(例如通过 -ffast-math)。

请注意GCC提供 __builtin_powi(x,n) 作为替代 pow( ),应生成内联乘法树。如果您想要牺牲性能的准确性,但又不想启用快速数学运算,请使用它。


614
2018-06-22 22:39



另请注意,Visual C ++提供了pow()的“增强”版本。通过电话 _set_SSE2_enable(<flag>) 同 flag=1,如果可能,它将使用SSE2。这会稍微降低精度,但会提高速度(在某些情况下)。 MSDN: _set_SSE2_enable() 和 POW() - TkTech
@TkTech:任何降低的准确性都是由于Microsoft的实现,而不是所用寄存器的大小。有可能提供一个 正确的舍入  pow 如果库编写器如此激励,则仅使用32位寄存器。有基于SSE的 pow 实现是 更多 比大多数基于x87的实现更准确,并且还有一些实现在速度上折衷一些准确性。 - Stephen Canon
@TkTech:当然,我只想说明准确性的降低是由于图书馆作者做出的选择,而不是SSE使用所固有的。 - Stephen Canon
我很想知道你在这里使用什么作为计算相对误差的“黄金标准” - 我通常会预料到它会是 a*a*a*a*a*a,但显然不是这样! :) - j_random_hacker
@j_random_hacker:因为我在比较单精度结果,双精度就足以满足黄金标准 - 来自的错误一个一个一个一个以double计算的是* 小于任何单精度计算的误差。 - Stephen Canon


另一个类似的案例:大多数编译器都不会优化 a + b + c + d 至 (a + b) + (c + d) (这是一个优化,因为第二个表达式可以更好地流水线化)并将其评估为给定(即as (((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


152
2018-06-23 11:44



这不完全一样。 Changin乘法/除法的顺序(除以除以0)比和/减法的changin顺序更安全。在我看来,编译器应该尝试关联mults./divs。因为这样做会减少操作总数,除了性能增益之外,还可以获得精确的增益。 - GameDeveloper
@DarioOO:这不安全。乘法和除法与指数的加法和减法相同,并且改变顺序很容易导致临时值超过指数的可能范围。 (不完全相同,因为指数不会损失精度......但是表示仍然非常有限,并且重新排序可能导致无法代表的值) - Ben Voigt
我想你错过了一些微积分背景。对2个数进行乘法和除法会引入相同的误差量。虽然减去/增加2个数字可能会引入更大的误差,特别是当2个数字的数量级不同时,因此它比sub / add更安全的重新分配mul / divide,因为它引入了最终错误的微小变化。 - GameDeveloper
@DarioOO:风险与mul / div不同:重新排序要么在最终结果中产生微不足道的变化,要么指数在某个时刻溢出(之前不会有),结果大不相同(可能是+ inf或0)。 - Peter Cordes


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)

这样可以轻松输入功率(无需拼出6 as with parens),让你无需进行这种优化 -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

74
2018-06-23 10:07



找到最佳功率树可能很难,但由于它只对小功率感兴趣,显而易见的答案是预先计算一次(Knuth提供一个表达100)并使用该硬编码表(这就是gcc内部为powi做的) 。 - Marc Glisse
在现代处理器上,速度受到延迟的限制。例如,乘法的结果可能在五个周期后可用。在那种情况下,找到创造一些力量的最快方法可能会更棘手。 - gnasher729
您还可以尝试查找给出相对舍入误差的最低上限或最低平均相对舍入误差的幂树。 - gnasher729
Boost也支持这一点,例如升压::数学:: POW <6>(N);我认为它甚至试图通过提取常见因子来减少乘法次数。 - gast128
好主意 !我已经为析因预计算做了这个。 - Caduchon


因为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-math”开关告诉GCC“嘿,海湾合作委员会,我知道我在做什么!”。这将允许GCC将* a * a * a * a * a转换为不同的文本 - 它看起来与a * a * a * a * a * a不同 - 但仍然计算错误间隔内的数字A * A * A * A * A * A。这没关系,因为你已经知道你正在使用间隔,而不是理想的数字。


49
2018-03-29 06:51



浮点数是准确的。它们不一定完全符合您的预期。此外,epsilon技术本身就是如何解决现实中的事物的近似,因为真实的预期误差是相对于尾数的比例,即,你通常高达大约1 LSB,但是这可能会增加如果你不小心的话,每次操作都要执行,所以在做一些浮点数非常重要的事情之前,请咨询数值分析师。如果可能,请使用合适的库。 - Donal Fellows
@DonalFellows:IEEE标准要求浮点计算产生的结果最精确地匹配源操作数是精确值时的结果,但这并不意味着它们实际上是 代表 确切的价值。在许多情况下,将0.1f视为(1,677,722 +/- 0.5)/ 16,777,216更有帮助,这应该与该不确定性所暗示的小数位数一起显示,而不是将其视为精确数量(1,677,722 +/- 0.5)/ 16,777,216(应显示为24位小数)。 - supercat
@supercat:IEEE-754在浮点数据方面非常明确 做 代表确切的值;第3.2至3.4条是相关部分。当然,您可以选择另外解释它们,就像您可以选择解释一样 int x = 3 就是这个意思 x 是3 +/- 0.5。 - Stephen Canon
@supercat:我完全同意,但这并不意味着 Distance 并不完全等于其数值;这意味着数值只是建模的某些物理量的近似值。 - Stephen Canon
对于数值分析,如果您将浮点数解释为不是间隔,而是作为精确值(恰好不是您想要的值),您的大脑会感谢您。例如,如果x在4.5左右,误差小于0.1,并且计算(x + 1) - x,则“间隔”解释会给你一个0.8到1.2的间隔,而“精确值”解释告诉你你的结果将是1,双精度误差最多为2 ^( - 50)。 - gnasher729


当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读取时没有文件扩展名);并写信给stdout。

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

    ; 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

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


49
2018-06-27 21:03



这对于整数乘法是合法的,因为两个补码溢出是未定义的行为。如果会出现溢出,无论重新排序操作如何,它都会发生在某个地方。因此,没有溢出的表达式评估相同,溢出的表达式是未定义的行为,因此编译器可以更改溢出发生的点。 gcc这样做 unsigned int也是。 - Peter Cordes


没有海报提到浮动表达的收缩(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,因此遵循上面的段落)。在过去,后一种选择并没有阻止转型,这意味着海湾合作委员会在这一点上不符合要求: https://gcc.gnu.org/bugzilla/show_bug.cgi?id=37845


27
2018-06-23 12:44



长期存在的热门问题有时会显示其年龄。这个问题在2011年得到了回答,当时海湾合作委员会可以原因是不尊重当时最近的C99标准。当然现在是2014年,所以GCC ......唉。 - Pascal Cuoq
不过,你不应该在没有接受答案的情况下回答相对较新的浮点问题吗?咳嗽 stackoverflow.com/questions/23703408 咳嗽 - Pascal Cuoq
我发现它......令人不安的是gcc没有实现C99浮点pragma。 - David Monniaux


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


26
2018-06-21 18:52



浮点总是确定性的。 - Alice
@Alice似乎相当清楚Bjorn在这里使用'确定性'在代码意义上给出了不同平台和不同编译器版本等相同的结果(外部变量可能超出了程序员的控制) - 而不是缺乏运行时的实际数字随机性。如果你指出这不是对这个词的正确使用,我不打算与之争论。 - greggo
@greggo即使在你对他所说的内容的解释中,它仍然是错的;这就是IEEE 754的全部要点,为跨平台的大多数(如果不是全部)操作提供相同的特性。现在,他没有提到平台或编译器版本,如果你希望每个远程服务器/客户端上的每一个操作都是相同的,这将是一个有效的问题....但是从他的陈述中这并不明显。一个更好的词可能是“可靠相似”或其他东西。 - Alice
@Alice你通过争论语义来浪费每个人的时间,包括你自己的时间。他的意思很清楚。 - Lanaru
@Lanaru标准的全部意义是语义学;他的意思显然不明确。 - Alice


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

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

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

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


26
2018-01-03 16:40



正如另一位评论者所指出的那样,这是荒谬的,这是不正确的;差异可能是成本的一​​半到10%,如果在紧密循环中运行,这将转化为许多浪费的指令,以获得可能无关紧要的额外精度。当你做monte carlo时说你不应该使用标准FP就像是说你应该总是用飞机穿越国家;它忽略了许多外部性。最后,这不是一个不常见的优化;死代码分析和代码缩减/重构很常见。 - Alice


像“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)”。这将是编译器供应商应该做的唯一理智的事情。


22
2017-10-01 19:33



downvoters应该意识到这个答案非常好。我可以引用几十个来源和文档来支持我的答案,而且我可能比任何downvoter更多地参与浮点精度。在StackOverflow中添加其他答案未涵盖的缺失信息是完全合理的,因此要礼貌并解释原因。 - GameDeveloper
在我看来,斯蒂芬佳能的答案涵盖了你所说的话。您似乎坚持使用样条函数实现libms:它们通常使用参数减少(取决于正在实现的函数)加上单个多项式,其系数通过Remez算法的或多或少复杂变体获得。连接点的平滑度不被认为是值得追求的libm函数的目标(如果它们最终足够准确,无论路径分成多少块,它们都会自动完全平滑)。 - Pascal Cuoq
你的答案的后半部分完全忽略了编译器应该生成实现源代码所说的代码的时间点。当你的意思是“准确性”时,你也使用“精确”一词。 - Pascal Cuoq
感谢您的输入,我稍微纠正了答案,最后2行^^仍然存在新的内容 - GameDeveloper