我正在对一个科学应用程序进行数值优化。我注意到的一件事是,GCC将通过将调用pow(a,2)编译为a*a来优化它,但调用pov(a,6)并没有优化,实际上会调用库函数pow,这会大大降低性能。(相比之下,可执行icc的“英特尔C++编译器”将消除对pow(a,6)的库调用。)

我好奇的是,当我使用GCC 4.5.1和选项“-O3-lm-funroll-loops-msse4”将pow(a,6)替换为a*a*a*a*a*a时,它使用了5条多指令:

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也有类似的行为。

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


当前回答

因为浮点数学不是关联的。浮点乘法中操作数的分组方式会影响答案的数值精度。

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

其他回答

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-funcafe数学优化。但是,这种重新排序违反了IEEE-754,因此需要标记。

正如Peter Cordes在一篇评论中指出的,有符号整数可以在没有funsafe数学优化的情况下进行这种优化,因为它恰好在没有溢出的情况下有效,如果有溢出,则会出现未定义的行为。所以你得到

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

只需-O。对于无符号整数,这更容易,因为它们是2的模幂,因此即使在溢出的情况下也可以自由地重新排序。

Fortran(为科学计算而设计)有一个内置的幂运算符,据我所知,Fortran编译器通常会以与您描述的方式类似的方式优化整数幂的提升。不幸的是,C/C++没有幂运算符,只有库函数pow()。这并不妨碍智能编译器专门处理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),因此没有理由大惊小怪。

然后将其用作功率<6>(a)。

这样可以很容易地输入幂(不需要像用括号一样拼出6),并且可以在不使用数学的情况下进行这种优化,以防出现精度相关的情况,例如补偿求和(这是一个操作顺序至关重要的示例)。

您可能也会忘记这是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

对于功率<6>(a),

    mulsd   %xmm0, %xmm0
    movapd  %xmm0, %xmm1
    mulsd   %xmm0, %xmm1
    mulsd   %xmm0, %xmm1

另一个类似的情况是:大多数编译器不会将a+b+c+d优化为(a+b)+(c+d)(这是一种优化,因为第二个表达式可以更好地进行流水线处理),并按照给定的方式对其求值(即(((a+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.00000e-05 0.000000e+00

因为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*a”优化为“t=(a*a);t*t*t”,这需要更少的乘法运算。但不幸的是,GCC编译器不知道编写代码的程序员是否认为“a”是一个有或没有错误的数字。所以GCC只会做源代码的样子——因为这是GCC用“肉眼”看到的。

…一旦你知道自己是什么样的程序员,你就可以使用“-fast math”开关告诉GCC“嘿,GCC,我知道我在做什么!”。这将允许GCC将a*a*a*a*a*a转换为一段不同的文本-它看起来与a*a*a*a*a*a*a*b*a不同-但仍在a*a a*a a*a*a a*的错误间隔内计算一个数字。这是可以的,因为你已经知道你使用的是时间间隔,而不是理想的数字。

这个问题已经有了一些很好的答案,但为了完整起见,我想指出C标准的适用部分是5.1.2.2.3/15(与C++11标准中的1.9/9节相同)。本节指出,只有当运算符真的是结合的或可交换的时,才能重新组合它们。