Multiply 128-bit vectors of signed 16-bit integers, widening to 32-bit elements

I have 2 __m128i. Each contains 8x int16_t (signed)

__mm128i a = {a0,...,a7}
__mm128i b = {b0,...,b7}

I want to multiply 8 elements. The result of each multiply is int32_t so each register will hold only 4 results:

__mm128i c0 = {a0*b0,...,a3*b3}
__mm128i c1 = {a4*b4,...,a7*b7}

I did not find such intrinsic.

SSE2 has a 16-bit "mul_hi" that returns the high half of the product. The high 16-bits and low 16-bits are combined into 32-bits using unpack.

    __m128i lo = _mm_mullo_epi16(a, b);
    __m128i hi = _mm_mulhi_epi16(a, b);
    __m128i c0 = _mm_unpacklo_epi16(lo, hi);
    __m128i c1 = _mm_unpackhi_epi16(lo, hi);

@aqrit's answer points out that it's as easy as pmullw + pmulhw, and unpack lo / hi to combine 16-bit halves of each 32-bit product. Accept that one.

Mostly obsolete answer - only useful if odd/even interleaved is useful

(Or perhaps as an alternate building block for an AVX2 version that wants to widen to __m256i results - it's non-obvious how to do that efficiently for __m256i inputs. With __m128i inputs you should just widen each input to 8x 32-bit with vpmovzxwd (_mm256_cvtepu16_epi32) and use _mm256_madd_epi16.)

x86's only widening 16x16 -> 32-bit SIMD multiply is pmaddwd, which horizontally adds pairs of products. So it's not directly usable, but if one of the products produces zero, then you have a 32-bit product of two signed 16-bit integers within that dword element.

You could unpack (with zero-extension) both inputs (4 shuffles) to feed 2 pmaddwd instruction.
Unpacking with zero-extension is cheaper for the high half, since _mm_unpackhi_epi16 against _mm_setzero_si128() is cheaper than manually doing sign extension with an arithmetic right shift, or doing extra work to feed the high half to _mm_cvtepi16_epi32 (pmovsxwd).

You'd only need or want sign extension if you were using pmulld (32x32 -> 32-bit multiply), but you want to avoid that since it's slower on Intel CPUs (2 uops per instruction). https://uops.info/

Multiplying 16-bit signed integers zero-extended to 32-bit with pmaddwd is doing a0*b0 + 0*0 and so on in each element. If you sign-extended the inputs, you could be doing a0*b0 + -1*-1, which is why you need zero-extension if you're going to do this optimization.

If unpacking to odd/even works for you, that's more efficient since 0 x garbage = 0 so we only need 2 instructions (plus some register copies) to prepare 2 pairs of inputs, and none of them are shuffles. (Shifts within 32-bit words like psrld would require shifting both inputs to keep the multiplicands lined up correctly.)

   __m128i a, b;  // inputs
__m128i even_mask = _mm_set1_epi32(0x0000FFFF);
__m128i a_even = _mm_and_si128(a, even_mask);
__m128i b_odd = _mm_andnot_si128(even_mask, b);  
// only need to mask one input to each multiply: 0 x garbage = 0
__m128i prod_even = _mm_madd_epi16(a_even, b);  // { a0*b0, a2*b2, a4*b4, a6*b6 }
__m128i prod_odd  = _mm_madd_epi16(a, b_odd);   // { a1*b1, a3*b3, a5*b5, a7*b7 } 

We could just mask a both ways; both ways have end up having the same amount of movdqa instructions to copy a register, at least in a function that can destroy the input registers holding a and b. It does have the fun advantage that whichever input was ready last (in out-of-order exec), the other one had time to already get masked and be ready to exec one of the multiplies. This is more interesting on old CPUs with only 1/clock throughput for pmaddwd (e.g. Intel before Skylake; https://uops.info/).

That odd/even product might be a better starting point for what you want, since it only requires 2 more shuffles to interleave prod_even and prod_odd into prod_low and prod_high

  __m128i prod_low  = _mm_unpacklo_epi32(prod_even, prod_odd); // { a0*b0, a1*b1, a2*b2, a3*b3 }
__m128i prod_high = _mm_unpackhi_epi32(prod_even, prod_odd);

This is 2 bitwise ANDs, 2 pmaddwd, and 2x punpckdq.

Unpacking both inputs first would have cost 2x pmovzxwd, 2x punpckhwd, and 2x pmaddwd, so all 4 of the non-multiply instructions would be shuffles. If the surrounding code also needs any shuffles, less pressure on port 5 is good on Intel CPUs with their limited shuffle throughput (especially before Ice Lake).

Using PAND we need a vector constant, but the 4x shuffle way only needs a zeroed register. I didn't look at the asm to see if one way needs more movdqa vector register copies. The 4x shuffle way might have the advantage there thanks to pmovzxwd being a copy-and-shuffle.

The 4x shuffle way looks like this, and is probably more efficient unless surrounding code also uses a lot of shuffles that would create a bottleneck on shuffle execution unit throughput.

   __m128i a, b;  // inputs
__m128i a_lo = _mm_cvtepu16_epi32(a);
__m128i b_lo = _mm_cvtepu16_epi32(b);
__m128i a_hi = _mm_unpackhi_epi16(a, _mm_setzero_si128());  // unpack(zero, a) would work equally well but compile less efficiently
__m128i b_hi = _mm_unpackhi_epi16(b, _mm_setzero_si128());
__m128i prod_lo = _mm_madd_epi16(a_lo, b_lo);  // { a0*b0, a1*b1, a2*b2, a3*b3 }
__m128i prod_hi = _mm_madd_epi16(a_hi, b_hi);

I put them on Godbolt inside functions. That doesn't test how they'd inline into loops, e.g. the compilers can overwrite the constant they load. But even with that, the 4x shuffle way is fewer instructions. Clang "optimizes" the and/andnot to 2x pblendw, which is 1 uop for port 5 only on Intel CPUs before Ice Lake. That might still be better than loading a constant for one call, unlike in a loop.

But anyway, the 4x shuffle way has no extra movdqa register-copy instructions, so it costs less front-end bandwidth than the other way that distributes the real work across different ports.

@chtz suggested prod_odd = _mm_sub_epi32(_mm_madd_epi16(a, b), prod_even); instead of pandn before the multiplies. That has longer critical-path latency, but does save 2 movdqa register-copy instructions when we have to compile with AVX. https://godbolt.org/z/7Gxn67d8h