如何在gcc中启用sse3自动矢量化

问题是该complex类型不是SIMD友好的。我从来都不喜欢这种complex类型的风扇，因为它是一个复合对象，通常不会映射到原始类型或硬件中的单个操作（某些情况下不是x86硬件）。

为了使复杂的算术SIMD友好，您需要同时对多个复数进行运算。对于SSE，您需要一次处理四个复数。

我们可以使用GCC的向量扩展名来简化语法。

typedef float v4sf __attribute__ ((vector_size (16)));

然后我们可以删除数组和向量扩展的并集

typedef union {
  v4sf v;
  float e[4];
} float4

最后，我们定义一个由四个复数组成的块，如下所示

typedef struct {
  float4 x;
  float4 y;
} complex4;

其中x是四个实部，y是四个虚部。

一旦有了这个，我们就可以像这样一次多个4个复数

static complex4 complex4_mul(complex4 a, complex4 b) {
  return (complex4){a.x.v*b.x.v -a.y.v*b.y.v, a.y.v*b.x.v + a.x.v*b.y.v};
}

最后，我们修改了您的函数，使其一次可以处理四个复数。

complex4 f4(complex4 x[], int n) {
  v4sf one = {1,1,1,1};
  complex4 p = {one,one};
  for (int i = 0; i < n; i++) p = complex4_mul(p, x[i]);
  return p;
}

让我们看一下汇编（Intel语法），看它是否最佳

.L3:
    movaps  xmm4, XMMWORD PTR [rsi]
    add     rsi, 32
    movaps  xmm1, XMMWORD PTR -16[rsi]
    cmp     rdx, rsi
    movaps  xmm2, xmm4
    movaps  xmm5, xmm1
    mulps   xmm1, xmm3
    mulps   xmm2, xmm3
    mulps   xmm5, xmm0
    mulps   xmm0, xmm4
    subps   xmm2, xmm5
    addps   xmm0, xmm1
    movaps  xmm3, xmm2
    jne     .L3

这就是四个4宽乘法，一个4宽加法和一个4宽减法。变量p保留在寄存器中，并且仅从数组x中加载数组，就像我们想要的那样。

让我们看一下复数乘积的代数

{a, bi}*{c, di} = {(ac - bd),(bc + ad)i}

那就是四个乘法，一个加法和一个减法。

正如我在此答案中解释的那样，高效的SIMD代数通常与标量算法相同。因此，我们用四个4宽乘法，加法和减法代替了四个1宽乘法，加法和减法。这是您使用4宽SIMD可以做到的最好：四个以一个价格出售。

请注意，这不需要SSE以外的任何指令，并且没有其他SSE指令（FMA4除外）会更好。因此，在64位系统上，您可以使用进行编译-O3。

将其扩展到具有AVX的8宽SIMD并不容易。

使用GCC向量扩展的一个主要优点是，无需任何额外工作即可获得FMA。例如，如果您使用-O3 -mfma4主循环进行编译，

.L3:
    vmovaps xmm0, XMMWORD PTR 16[rsi]
    add     rsi, 32
    vmulps  xmm1, xmm0, xmm2
    vmulps  xmm0, xmm0, xmm3
    vfmsubps        xmm1, xmm3, XMMWORD PTR -32[rsi], xmm1
    vmovaps xmm3, xmm1
    vfmaddps        xmm2, xmm2, XMMWORD PTR -32[rsi], xmm0
    cmp     rdx, rsi
    jne     .L3