ICC是否符合C99规范中关于复数乘法的要求？

Question

考虑这段简单的代码：

TBD

#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复数乘法的规范吗？

- Simd

---

@Zboson 这就是问题中的代码。 - Simd

GCC调用了__mulsc3https://godbolt.org/g/9rcINV。我该如何让GCC将对`__mulsc3`的调用内联？使用`-Ofast`时，GCC可以内联它，如https://godbolt.org/g/7vq9su。 - Z boson

@eleanora，我知道这是问题中的代码。我添加了链接，以防有人想玩一下它。 - Z boson

如果你的目标是效率，那么不要使用 complex。在我看来，这是一种愚蠢的类型，因为它是为了让 C 语言也有类似 Fortran 中的复数类型而添加的。我喜欢 C 语言的一个原因是它的类型是基本类型，并且通常直接映射到汇编寄存器中。但是 complex 是一种像 C++ 中的类一样的复合类型。它在 C 语言中似乎有些奇怪。我想可能有硬件可以直接实现复数类型，但我从未使用过任何我知道的这样的硬件。 - Z boson

让我们在聊天中继续这个讨论。 - Simd

Answer 1

代码不符合规范。

Annex G，Section 5.1，第4段的内容如下：

对于所有实数、虚数和复数操作数，运算符*和/满足以下无穷大性质：

— 如果一个操作数是无穷大，另一个操作数是非零有限数或无穷大，则*运算符的结果是无穷大；

因此，如果z=a * ib无限大，w=c * id无限大，则数字z * w必须无限大。

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

至少有一个无限部分的复数或虚数值被视为无穷大（即使其其他部分是NaN）。

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

然而，如果我们让z=∞ + i∞（一个无限值），w=i∞（一个无限值），则由于∞·0的中间结果，Intel代码的结果为z * w=NaN + iNaN²。

这足以将其标记为不符合规范。

---

我们可以通过查看第一引用的脚注来进一步确认此问题（未在此处报告脚注），它提到了CX_LIMITED_RANGE编译指示。

Section 7.3.4，第1段的内容如下：

复数相乘、除和绝对值的通常数学公式存在问题，因为它们处理无穷大的方式不当，并且存在过度溢出和下溢的问题。CX_LIMITED_RANGE编译指示可用于告知实现，在状态为“on”时，可以接受通常的数学公式[会生成NaN]。

在这里，标准委员会试图缓解复数乘法（和除法）的巨大工作量。
事实上，GCC有一个控制此行为的标志：

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

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

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

仅凭这个选项，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‖ 时，仍然保持无限定义为无界极限的结果。简单来说，在复平面上有一个完整的无穷圆周，并且只需要一个“坐标”变为无穷大就可以得到无穷大的模。

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

- Margaret Bloom

---

@eleanora 当一个组件是NaN而另一个是INFINITY时，你应该得到一个错误的结果。 - Margaret Bloom

@eleanora 如果你在ICC函数中使用INF + i NAN作为输入，则结果是NaN + i NaN，而标准要求它应该至少是一个无穷大分量的复数。GCC的处理方式是正确的。然而，测试这些边缘情况并不容易，因为ICC和我的GCC使用了两种不同的复数和GCC约定。我直接在汇编中测试了英特尔函数，而对于GCC，我必须生成INF + i NAN作为(INF + I * 0) * (1 + I * INF)，因为在GCC中INF + i * NAN会得到一个NAN + i NAN的结果（因为这被视为从NaN虚数和实数之和创建复数）。 - Margaret Bloom

英特尔似乎已经承认存在一个漏洞。"我可以重现您报告的问题。有两个警告显示INFINITY*I被视为复数类型超出范围的数字。这应该是不正确的。我已经升级了这个问题，并将其输入到我们的问题跟踪系统中。我会在我有关于这个问题的更新时通知您。" - Simd

NaN的逻辑数学解释通常是一个值，它可能高于最大可表示值，低于最小可表示值或介于两者之间。如果X和Y都是未知量，且高于最大可表示值，则X-Y为NaN，因为X+Y可以任意高或低。如果将（X+Yi）平方，将得到一个实部可以是任意大的正数或负数的值。如果再将其平方，那么两个部分都可以是任意大的正数或负数。 - supercat

能够将（NaN，+Inf）解释为“实部未知且虚部非常高的值”，并对于其他三种NaN和无穷大的组合也是如此，这似乎很有用，但要求这些值的乘法必须产生这样的值将会削弱这种区别。 - supercat

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标准取决于英特尔的编译器是否被记录为符合Annex F和G（这些附件指定C语言实现提供符合IEEE标准的实数和复数算术的意义），如果是，则需要给出哪些命令行选项以获得符合模式。

- zwol

---

您IP地址为143.198.54.68，由于运营成本限制，当前对于免费用户的使用频率限制为每个IP每72小时10次对话，如需解除限制，请点击左下角设置图标按钮（手机用户先点击左上角菜单按钮）。 - Simd