ICC是否满足复数的乘法C99规范?

考虑这个简单的代码:

#include  complex float f(complex float x) { return x*x; } 

如果使用英特尔编译器使用-O3 -march=core-avx2 -fp-model strict进行编译,则可以得到:

 f: vmovsldup xmm1, xmm0 #3.12 vmovshdup xmm2, xmm0 #3.12 vshufps xmm3, xmm0, xmm0, 177 #3.12 vmulps xmm4, xmm1, xmm0 #3.12 vmulps xmm5, xmm2, xmm3 #3.12 vaddsubps xmm0, xmm4, xmm5 #3.12 ret 

这比你从gccclang获得的代码简单得多,而且比你在网上找到的用于乘以复数的代码简单得多。 例如,它没有明确地用于处理复杂的NaN或无穷大。

这个组件是否符合C99复数乘法的规范?

代码不符合要求。

附件G第5.1节第4段内容如下

对于所有实数,虚数和复数操作数, */运算符满足以下无穷大属性:

– 如果一个操作数是无穷大而另一个操作数是非零有限数或无穷大,那么*运算符的结果是无穷大;

因此,如果z = a * i b是无穷大且w = c * i d是无穷大,则数字z * w必须是无穷大。

同一附件第3节第1段定义了复数对无穷大的意义:

具有至少一个无限部分的复数或虚数值被视为无穷大(即使其另一部分是NaN)。

因此,如果a或b是z,z是无穷大。
这确实是一个明智的选择,因为它反映了数学框架1

然而,如果我们让z =∞+i∞(无穷大值)和w =i∞(和无穷大值),由于∞·0中间体2 ,英特尔代码的结果是z * w = NaN + i NaN。

这足以将其标记为不符合要求。


我们可以通过查看第一个引用的脚注(此处未报告脚注)进一步确认这一点,它提到了CX_LIMITED_RANGE编译指令。

第7.3.4节,第1段内容如下

复数乘法,除法和绝对值的通常数学公式是有问题的,因为它们对无穷大的处理以及由于不适当的溢出和下溢。 CX_LIMITED_RANGE pragma可用于通知实现(状态为”on”)通常的数学公式[产生NaN]是可接受的。

在这里,标准委员会正试图减轻复杂乘法(和除法)的巨大工作量。
实际上GCC有一个标志来控制这种行为 :

-fcx-limited-range
启用时,此选项表示执行复杂分割时不需要范围缩减步骤。

此外,没有检查复数乘法或除法的结果是否是NaN + I * NaN,试图在这种情况下挽救情况。

默认值为-fno-cx-limited-range ,但-ffast-math启用
此选项控制ISO C99 CX_LIMITED_RANGE pragma的默认设置。

只有这个选项才能使GCC生成慢速代码和额外的检查 ,没有它生成的代码与Intel的代码有相同的缺陷(我将源代码转换为C ++)

 f(std::complex): movq QWORD PTR [rsp-8], xmm0 movss xmm0, DWORD PTR [rsp-8] movss xmm2, DWORD PTR [rsp-4] movaps xmm1, xmm0 movaps xmm3, xmm2 mulss xmm1, xmm0 mulss xmm3, xmm2 mulss xmm0, xmm2 subss xmm1, xmm3 addss xmm0, xmm0 movss DWORD PTR [rsp-16], xmm1 movss DWORD PTR [rsp-12], xmm0 movq xmm0, QWORD PTR [rsp-16] ret 

如果没有它,代码就是

 f(std::complex): sub rsp, 40 movq QWORD PTR [rsp+24], xmm0 movss xmm3, DWORD PTR [rsp+28] movss xmm2, DWORD PTR [rsp+24] movaps xmm1, xmm3 movaps xmm0, xmm2 call __mulsc3 movq QWORD PTR [rsp+16], xmm0 movss xmm0, DWORD PTR [rsp+16] movss DWORD PTR [rsp+8], xmm0 movss xmm0, DWORD PTR [rsp+20] movss DWORD PTR [rsp+12], xmm0 movq xmm0, QWORD PTR [rsp+8] add rsp, 40 ret 

并且__mulsc3function与C99推荐的复数乘法实际上相同。
它包括上述检查。


1数字的模数从实例扩展| z | 复杂的“‖z”,保持无限的定义是无限制的结果。 简单地说,在复杂的平面中,存在无限值的整个圆周,并且只需要一个“坐标”就可以获得无限模数。

2如果我们记住z = NaN +i∞或z =∞+ i NaN是有效的无限值,情况会变得最糟

我得到类似但不完全相同的代码来自clang 3.8 at -O2 -march=core-avx2 -ffast-math :我对后期的x86浮点function并不十分熟悉,但我认为它的计算方法相同但是使用不同的指令来重排寄存器中的值。

 f: vmovshdup %xmm0, %xmm1 # xmm1 = xmm0[1,1,3,3] vaddss %xmm0, %xmm0, %xmm2 vmulss %xmm2, %xmm1, %xmm2 vmulss %xmm1, %xmm1, %xmm1 vfmsub231ss %xmm0, %xmm0, %xmm1 vinsertps $16, %xmm2, %xmm1, %xmm0 # xmm0 = xmm1[0],xmm2[0],xmm1[2,3] retq 

具有相同选项的GCC 6.3似乎再次进行相同的计算,但是以第三种方式改变了值:

 f: vmovq %xmm0, -8(%rsp) vmovss -4(%rsp), %xmm2 vmovss -8(%rsp), %xmm0 vmulss %xmm2, %xmm2, %xmm1 vfmsub231ss %xmm0, %xmm0, %xmm1 vmulss %xmm2, %xmm0, %xmm0 vmovss %xmm1, -16(%rsp) vaddss %xmm0, %xmm0, %xmm0 vmovss %xmm0, -12(%rsp) vmovq -16(%rsp), %xmm0 ret 

如果没有-ffast-math ,两个编译器都会生成截然不同的代码,至少可以检查NaN。

我从中得出结论,即使使用-fp-model strict ,英特尔的编译器也不会生成完全符合IEEE标准的复数乘法。 可能还有一些其他命令行开关,使其生成完全符合IEEE标准的代码。

这是否符合违反C99的要求取决于英特尔的编译器是否符合附件F和G(其中规定了实现C以提供符合IEEE标准的实数和复数运算的含义),如果是,您必须提供什么命令行选项才能获得符合模式。