我正在对一个科学应用程序进行数值优化。我注意到的一件事是,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也有类似的行为。
为什么编译器不认识这种优化技巧?
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为整数时,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 ; 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
正如其他海报所指出的,在浮点运算中,这个选项是不可能的,因为浮点运算是不相关的。
像“pow”这样的库函数通常经过精心设计,以产生最小可能的错误(在一般情况下)。这通常是用样条逼近函数实现的(根据Pascal的评论,最常见的实现似乎是使用Remez算法)
基本上是以下操作:
pow(x,y);
具有与任何单个乘法或除法中的误差大致相同大小的固有误差。
执行以下操作时:
float a=someValue;
float b=a*a*a*a*a*a;
其固有误差大于单个乘法或除法的误差的5倍(因为您组合了5个乘法)。
编译器应该非常小心它正在进行的优化:
如果将pow(a,6)优化为a*a*a*a*a*a,可能会提高性能,但会大大降低浮点数的精度。如果将a*a*a*a*a*a优化为pow(a,6),实际上可能会降低精度,因为“a”是一个特殊的值,它允许无误差的乘法(2的幂或一些小整数)如果将pow(a,6)优化为(a*a*a)*(a*a*a)或(a*a)*。
一般来说,您知道对于任意浮点值,“pow”的精度比您最终可以编写的任何函数都要高,但在某些特殊情况下,多次乘法可能具有更好的精度和性能,这取决于开发人员选择更合适的方法,最终对代码进行注释,以便其他人不会“优化”该代码。
唯一有意义的优化(个人观点,显然是GCC中没有任何特定优化或编译器标志的选择)应该是将“pow(a,2)”替换为“a*a”。这将是编译器供应商应该做的唯一明智的事情。