+/*
+ * Provide an SSE implementation of _mm_mul_epi32() in terms of
+ * _mm_mul_epu32().
+ *
+ * Basically, albeit surprising at first (and second, and third...) look
+ * if a * b is done signed instead of unsigned, can just
+ * subtract b from the high bits of the result if a is negative
+ * (and the same for a if b is negative). Modular arithmetic at its best!
+ *
+ * So for int32 a,b in crude pseudo-code ("*" here denoting a widening mul)
+ * fixupb = (signmask(b) & a) << 32ULL
+ * fixupa = (signmask(a) & b) << 32ULL
+ * a * b = (unsigned)a * (unsigned)b - fixupb - fixupa
+ * = (unsigned)a * (unsigned)b -(fixupb + fixupa)
+ *
+ * This does both lo (dwords 0/2) and hi parts (1/3) at the same time due
+ * to some optimization potential.
+ */
+static inline __m128i
+mm_mullohi_epi32(const __m128i a, const __m128i b, __m128i *res13)
+{
+ __m128i a13, b13, mul02, mul13;
+ __m128i anegmask, bnegmask, fixup, fixup02, fixup13;
+ a13 = _mm_shuffle_epi32(a, _MM_SHUFFLE(2,3,0,1));
+ b13 = _mm_shuffle_epi32(b, _MM_SHUFFLE(2,3,0,1));
+ anegmask = _mm_srai_epi32(a, 31);
+ bnegmask = _mm_srai_epi32(b, 31);
+ fixup = _mm_add_epi32(_mm_and_si128(anegmask, b),
+ _mm_and_si128(bnegmask, a));
+ mul02 = _mm_mul_epu32(a, b);
+ mul13 = _mm_mul_epu32(a13, b13);
+ fixup02 = _mm_slli_epi64(fixup, 32);
+ fixup13 = _mm_and_si128(fixup, _mm_set_epi32(-1,0,-1,0));
+ *res13 = _mm_sub_epi64(mul13, fixup13);
+ return _mm_sub_epi64(mul02, fixup02);
+}