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 * ib 是无限的并且 w = c * id 是无限的，数字 z * w 必须是无限的。

同一附件，第 3 节，第 1 段定义了复数无限的含义：

具有至少一个无限部分的复数或虚数被视为无限（即使它的另一部分是 NaN）。

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

但是，如果我们让 z = ∞ + i∞（无限值）并且 w = i ∞（和无限值）英特尔代码的结果是 z * w = NaN + iNaN 由于 ∞ · 0 个中间体2.

这足以将其标记为不合格。

---

我们可以通过查看第一个引用上的脚注来进一步确认这一点（这里没有报告脚注），它提到了CX_LIMITED_RANGE pragma 指令。

第 7.3.4 节，第 1 段改为

复数乘法、除法和绝对值的常用数学公式存在问题，因为它们处理无穷大以及过度的上溢和下溢。 CX_LIMITED_RANGE pragma 可用于通知实现（在状态为“开启”的情况下）[产生 NaN] 的常用数学公式是可接受的。

标准委员会正在努力减轻复杂乘法（和除法）的巨大工作量。
In fact GCC has a flag to control this behaviour:

-fcx-limited-range
启用后，此选项表示执行复杂除法时不需要范围缩小步骤。

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

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

仅此选项makes GCC generate slow code and additional checks，如果没有它，它生成的代码与 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 function 实际上与标准 C99 推荐的复数乘法相同。
它包括上述检查。

---

¹ 其中一个数的模是从真实情况 |z| 扩展而来的到复数“z”，保持无限的定义是无限限制的结果。简单地说，在复平面中有一整周的无限值，只需要一个“坐标”就可以得到无限的模数。

² 如果我们记住 z = NaN + i∞ 或 z = ∞，情况会变得更糟+ iNaN 是有效的无限值

Answer 2

我在-O2 -march=core-avx2 -ffast-math 处获得了来自 clang 3.8 的相似但不完全相同的代码：我对现代 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。

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

这是否符合C99的条件取决于英特尔的编译器是否被记录为符合附件 F 和 G（它们指定了 C 实现提供符合 IEEE 标准的含义实数和复数算术），如果是，您必须提供哪些命令行选项才能获得一致模式。