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

Question

考虑这个简单的代码:

#include <complex.h>
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 

这比从两者中获得的代码简单得多gcc,clang而且比在线复制数字的代码简单得多.例如,它没有明确地用于处理复杂的NaN或无穷大.

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

Answer 1

代码不符合要求.

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

的*和/运营商满足所有实数,虚,和复杂的操作数以下无穷属性:

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

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

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

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

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

但是如果我们让Ž =∞+ 我 ∞(无穷大值)和瓦特 = 我 ∞(和无限值)的结果为Intel代码是ž*瓦特 = NaN的+ 我的NaN由于∞·0的中间体².

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

我们可以通过查看第一个引用的脚注(这里没有报告脚注)进一步证实这一点,它提到了CX_LIMITED_RANGEpragma指令.

第7.3.4节,第1段内容

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

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

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

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

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

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

f(std::complex<float>):
        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<float>):
        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

并且该__mulsc3功能实际上与标准C99推荐的复数乘法相同.
它包括上述检查.

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

²,情况得到最坏如果我们记得ž = NaN的+ 我 ∞或Z ^ =∞+ 我的NaN是有效的无限价值

Answer 2

我得到了类似但不完全相同的代码来自于clang 3.8 -O2 -march=core-avx2 -ffast-math:我对最近的x86浮点特性并不十分熟悉,但我认为它正在进行相同的计算,但是使用不同的指令在寄存器中对值进行混洗.

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时,至少.

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

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