title: How to write portable simd code for complex multiplicative reduction
body: <p>I want to write fast simd code to compute the multiplicative reduction of a complex array. In standard C this is:</p>
  
  <pre><code>#include &lt;complex.h&gt;
  complex float f(complex float x[], int n ) {
     complex float p = 1.0;
     for (int i = 0; i &lt; n; i++)
        p *= x[i];
     return p;
  }
  </code></pre>
  
  <p><code>n</code> will be at most 50.</p>
  
  <p>Gcc can't auto-vectorize complex multiplication but, as I am happy to assume the gcc compiler and if I knew I wanted to target sse3 I could follow <a href="https://stackoverflow.com/questions/41417061/how-to-enable-sse3-autovectorization-in-gcc">How to enable sse3 autovectorization in gcc</a> and write:</p>
  
  <pre><code>typedef float v4sf __attribute__ ((vector_size (16)));
  typedef union {
    v4sf v;
    float e[4];
  } float4
  typedef struct {
    float4 x;
    float4 y;
  } complex4;
  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 &lt; n; i++) p = complex4_mul(p, x[i]);
    return p;
  }
  </code></pre>
  
  <p>This indeed produces fast vectorized assembly code using gcc. Although you still need to pad your input to a multiple of 4. The assembly you get is:</p>
  
  <pre><code>.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
  </code></pre>
  
  <p>However, it is designed for the exact simd instruction set and is not optimal for avx2 or avx512 for example for which you need to change the code.</p>
  
  <blockquote>
    <p>How can you write C or C++ code for which gcc will produce optimal
    code when compiled for any of sse, avx2 or avx512?  That is, do you always have to write separate functions by hand for each different width of SIMD register?</p>
    
    <p>Are there any open source libraries that make this easier?</p>
  </blockquote>