GCC 为什么不将 a * 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)的库调用。)

我很好奇的是,当我使用 GCC 4.5.1 和选项 “ -O3 -lm -funroll-loops -msse4 ” 将a*a*a*a*a*a替换为pow(a,6)时,它会使用 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

这将乘法指令的数量减少到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)之间的所有单精度浮点数的 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( )的替代方法,后者应生成一个内联乘法树。如果您要在准确性与性能之间进行权衡,但又不想启用快速计算,请使用该选项。

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

这样可以很容易地输入幂(不需要用括号来拼写 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

当 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 读取时没有文件扩展名);并写入标准输出。

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

; 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

我正在 Ubuntu Mind 衍生版 Linux Mint 16 Petra 上使用系统 GCC。这是 gcc 版本:

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

正如其他张贴者所指出的,在浮点数中此选项是不可能的,因为浮点数算法不具有关联性。

因为 32 位浮点数(例如 1.024)不是 1.024。在计算机中,1.024 是一个间隔:从(1.024-e)到(1.024 + e),其中 “e” 表示错误。有些人没有意识到这一点,他们还认为 a * 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” 是带错误或不带错误的数字。因此,GCC 只会执行源代码的样子 - 因为这就是 GCC 的 “裸眼”。

...... 一旦你知道那种程序员的是什么,你可以使用 “-ffast - 数学” 开关告诉 GCC 说:“嘿,GCC,我知道我在做什么!”。这将使 GCC 可以将 a * a * a * a * a * a 转换为不同的文本 - 它看起来与 a * a * a * a * a * a * a 不同 - 但仍会计算错误间隔为 a * a * a * a * a * a。可以,因为您已经知道自己正在使用间隔而不是理想数字。

尚无张贴者提到浮动表达式的收缩(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,因此遵循上一段)。过去,后一种选择不会阻止转换,这意味着 GCC 在这一点上不符合要求: https//gcc.gnu.org/bugzilla/show_bug.cgi?id = 37845

我根本不希望这种情况得到优化。表达式包含子表达式的情况很少见,这些子表达式可以重新组合以删除整个操作。我希望编译器作者将时间投入到更可能导致显着改进的领域上,而不是覆盖很少遇到的边缘情况。

从其他答案中得知,使用适当的编译器开关确实可以优化此表达式,这让我感到惊讶。优化要么是微不足道的,要么是更常见的优化的边缘案例,要么是编译器编写者非常彻底。

正如您在此处所做的那样,向编译器提供提示没有错。重新排列语句和表达式,以了解它们将带来什么不同,这是微优化过程中正常且预期的部分。

尽管考虑到两个表达式传递不一致的结果(没有适当的切换)可能会证明编译器是合理的,但您不必受此限制的约束。差异将非常小,以至于如此之大,以至于如果差异对您很重要,那么您就不应首先使用标准浮点算法。

正如 Lambdageek 指出的那样,浮点乘法不是关联的,因此精度可能会降低,但是当精度更高时,您可能会反对优化,因为您需要确定性的应用程序。例如,在游戏模拟客户端 / 服务器中,每个客户端都必须模拟您希望确定点浮点计算的同一世界。

通常会精心设计诸如 “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*apow(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”。那将是编译器供应商应该做的唯一明智的事情。