我正在对科学应用程序进行一些数值优化。 我注意到的一件事是,GCC将通过将其编译为a*a
来优化调用pow(a,2)
,但是调用pow(a,6)
并未经过优化,实际上会调用库函数pow
,这大大降低了速度表演。 (相反, 英特尔C ++编译器 (可执行文件icc
)将消除对pow(a,6)
的库调用。)
我很好奇的是,当我使用GCC 4.5.1和选项“-O3 -lm -funroll-loops -msse4
”将pow(a,6)
替换为a*a*a*a*a*a
,它使用了5mulsd
说明:
movapd %xmm14, %xmm13mulsd %xmm14, %xmm13mulsd %xmm14, %xmm13mulsd %xmm14, %xmm13mulsd %xmm14, %xmm13mulsd %xmm14, %xmm13
而如果我写(a*a*a)*(a*a*a)
,它将产生
movapd %xmm14, %xmm13mulsd %xmm14, %xmm13mulsd %xmm14, %xmm13mulsd %xmm13, %xmm13
这将乘法指令的数量减少到icc
具有相似的行为。
为什么编译器无法识别此优化技巧?
#1楼
这个问题已经有了一些很好的答案,但是为了完整起见,我想指出,C标准的适用部分是5.1.2.2.3 / 15(与C.1中的1.9 / 9相同)。 C ++ 11标准)。 本节指出,只有运算符确实是关联的或可交换的,才可以重新组合。
#2楼
当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 = ximul eax, edi ; temp = x * tempimul eax, edi ; temp = x * tempimul eax, eax ; temp = temp * temp
我正在Ubuntu Mind衍生版Linux Mint 16 Petra上使用系统GCC。 这是gcc版本:
$ gcc --versiongcc (Ubuntu/Linaro 4.8.1-10ubuntu9) 4.8.1
正如其他张贴者所指出的,在浮点数中此选项是不可能的,因为浮点数算法不具有关联性。
#3楼
尚无张贴者提到浮动表达式的收缩(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 : ///bugzilla/show_bug.cgi?id = 37845
#4楼
通常会精心设计诸如“ 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)*(a*a)*(a*a)
则仍然会损失精度与pow
功能相比。
通常,您知道对于任意浮点值,“ pow”的精度要比您最终可以编写的任何函数更好,但是在某些特殊情况下,多次乘法可能具有更好的精度和性能,这取决于开发人员选择更合适的值,最终对代码进行注释,以使其他人都无法“优化”该代码。
唯一有意义的事情(个人观点,显然是在GCC中没有任何特定的优化或编译器标记的选择)以进行优化,应该将“ pow(a,2)”替换为“ a * a”。 那将是编译器供应商应该做的唯一明智的事情。
#5楼
gcc实际上可以进行此优化,即使对于浮点数也是如此。 例如,
double foo(double a) {return a*a*a*a*a*a;}
变成
foo(double):mulsd %xmm0, %xmm0movapd %xmm0, %xmm1mulsd %xmm0, %xmm1mulsd %xmm1, %xmm0ret
-O -funsafe-math-optimizations
。 但是,此重新排序违反了IEEE-754,因此需要该标志。
正如彼得·科德斯(Peter Cordes)在评论中指出的那样,带符号整数可以在没有-funsafe-math-optimizations
optimizations的情况下进行此优化,因为它在没有溢出时以及在有溢出的情况下都具有不确定的行为。 所以你得到
foo(long):movq %rdi, %raximulq %rdi, %raximulq %rdi, %raximulq %rax, %raxret
与-O
。 对于无符号整数,这甚至更容易,因为它们的mod幂为2,因此即使面对溢出也可以自由地重新排序。
#6楼
我根本不希望这种情况得到优化。 表达式包含子表达式的情况很少见,这些子表达式可以重新组合以删除整个操作。 我希望编译器作者将时间投入到更可能导致显着改进的领域上,而不是覆盖很少遇到的边缘情况。
从其他答案中得知,使用适当的编译器开关确实可以优化此表达式,这让我感到惊讶。 优化要么是微不足道的,要么是更常见的优化的边缘案例,要么是编译器编写者非常彻底。
正如您在此处所做的那样,向编译器提供提示没有错。 重新排列语句和表达式,以了解它们将带来什么不同,这是微优化过程中正常且预期的部分。
尽管考虑到两个表达式传递不一致的结果(没有适当的切换)可能会证明编译器是合理的,但您不必受此限制的约束。 差异将非常小,以至于如此之大,以至于如果差异对您很重要,那么您就不应首先使用标准浮点算法。
#7楼
因为浮点数学不是关联的 。 以浮点乘法将操作数分组的方式会影响答案的数值精度。
结果,大多数编译器在对浮点计算进行重新排序时非常保守,除非他们可以确保答案保持不变,或者除非您告诉他们您不关心数值精度。 例如:gcc 的-fassociative-math
选项允许gcc重新关联浮点运算,甚至-ffast-math
选项允许更加精确地权衡速度。
#8楼
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-08worst relative error using (a*a*a)*(a*a*a): 2.94e-07worst relative error usinga*a*a*a*a*a: 2.58e-07
使用pow
而不是乘法树可将错误范围限制为4。 除非经过用户许可(例如,通过-ffast-math
),否则编译器不应(并且通常不会)进行“优化”以增加错误。
请注意,GCC提供了__builtin_powi(x,n)
作为pow( )
的替代方法,后者应生成一个内联乘法树。 如果您要在准确性与性能之间进行权衡,但又不想启用快速计算,请使用该选项。
#9楼
另一个类似的情况:大多数编译器不会将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
#10楼
因为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。 可以,因为您已经知道自己正在使用间隔而不是理想数字。
#11楼
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)
。
这样可以很容易地输入幂(不需要用括号来拼写6a
s),并且可以在不-ffast-math
的情况下进行这种优化,以防万一您对精度有依赖性,例如补偿求和 (例如操作至关重要)。
您可能还会忘记这是C ++,并且仅在C程序中使用它(如果它使用C ++编译器进行编译)。
希望这会有用。
编辑:
这是我从编译器得到的:
对于a*a*a*a*a*a
,
movapd %xmm1, %xmm0mulsd %xmm1, %xmm0mulsd %xmm1, %xmm0mulsd %xmm1, %xmm0mulsd %xmm1, %xmm0mulsd %xmm1, %xmm0
对于(a*a*a)*(a*a*a)
,
movapd %xmm1, %xmm0mulsd %xmm1, %xmm0mulsd %xmm1, %xmm0mulsd %xmm0, %xmm0
对于power<6>(a)
,
mulsd %xmm0, %xmm0movapd %xmm0, %xmm1mulsd %xmm0, %xmm1mulsd %xmm0, %xmm1
#12楼
正如Lambdageek指出的那样,浮点乘法不是关联的,因此精度可能会降低,但是当精度更高时,您可能会反对优化,因为您需要确定性的应用程序。 例如,在游戏模拟客户端/服务器中,每个客户端都必须模拟您希望确定点浮点计算的同一世界。