在SSE2上如何签名32位加宽乘法？

问题描述

在审查WebAssembly SIMD扩展扩展提案时出现了这个问题。

要支持较旧的硬件，我们需要支持SSE2，并且32位整数的唯一矢量乘法运算是pmuludq。 （仅在SSE4.1中添加了签名的pmuldq）

（非扩展相对容易；随机播放2倍pmuludq，并取4个结果的低半部分来模拟SSE4.1 pmulld）。

解决方法

mulhs(a,b) = mulhu(a,b) - (a < 0 ? b : 0) - (b < 0 ? a : 0)

使用这种方法，可以计算出两个带符号的全角乘积，

__m128i mul_epi32(__m128i a,__m128i b) {
    a = _mm_shuffle_epi32(a,_MM_SHUFFLE(3,1,0));
    b = _mm_shuffle_epi32(b,0));
    __m128i unsignedProduct = _mm_mul_epu32(a,b);
    __m128i threshold = _mm_set_epi32(INT_MIN,INT_MIN,0);
    __m128i signA = _mm_cmplt_epi32(a,threshold);
    __m128i signB = _mm_cmplt_epi32(b,threshold);
    __m128i x = _mm_shuffle_epi32(_mm_and_si128(signA,b),_MM_SHUFFLE(2,3,1));
    __m128i y = _mm_shuffle_epi32(_mm_and_si128(signB,a),1));
    return _mm_sub_epi32(_mm_sub_epi32(unsignedProduct,x),y);
}

与其他提案相比，这节省了一些操作，但它非常接近，现在包含了一个负载，如果此代码很冷，则可能会很糟糕。

,

在## asm上向@GeDaMo大喊大叫，以帮助提出此解决方案。

Godbolt

C / C ++：

#include <xmmintrin.h>
#include <stdint.h>
#include <tmmintrin.h>
#include <smmintrin.h>
#include <cstdio>

typedef int32_t int32x4_t __attribute__((vector_size(16))) __attribute__((aligned(16)));
typedef int64_t int64x2_t __attribute__((vector_size(16))) __attribute__((aligned(16)));

int64x2_t multiply32_low_s(int32x4_t a,int32x4_t b) {
    auto aSigns = a >> 31;
    auto bSigns = b >> 31;
    auto aInt = a ^ aSigns;
    aInt -= aSigns;
    auto bInt = b ^ bSigns;
    bInt -= bSigns;
    const auto shuffleMask = _MM_SHUFFLE(1,0);
    auto absProd = _mm_mul_epu32(_mm_shuffle_epi32((__m128i)aInt,shuffleMask),_mm_shuffle_epi32((__m128i)bInt,shuffleMask));
    auto aSignsInt = _mm_shuffle_epi32((__m128i)aSigns,shuffleMask);
    auto bSignsInt = _mm_shuffle_epi32((__m128i)bSigns,shuffleMask);
    auto prodSigns = aSignsInt ^ bSignsInt;
    absProd ^= prodSigns;
    absProd -= prodSigns;
    return (int64x2_t)absProd;
}

int64x2_t multiply32_high_s(int32x4_t a,int32x4_t b) {
    auto aSigns = a >> 31;
    auto bSigns = b >> 31;
    auto aInt = a ^ aSigns;
    aInt -= aSigns;
    auto bInt = b ^ bSigns;
    bInt -= bSigns;
    const auto shuffleMask = _MM_SHUFFLE(3,2,2);
    auto absProd = _mm_mul_epu32(_mm_shuffle_epi32((__m128i)aInt,shuffleMask);
    auto prodSigns = aSignsInt ^ bSignsInt;
    absProd ^= prodSigns;
    absProd -= prodSigns;
    return (int64x2_t)absProd;
}


int main(int argc,char* argv[]) {
    int32x4_t a{-5,500,-5000,50000};
    int32x4_t b{10,-100,500000000};
    auto c = multiply32_low_s(a,b);
    auto d = multiply32_high_s(a,b);
    printf("%ld %ld\n",c[0],c[1]);
    printf("%ld %ld\n",d[0],d[1]);
}

组装

multiply32_low_s(int __vector(4),int __vector(4)):
 movdqa xmm3,xmm0
 movdqa xmm2,xmm1
 psrad  xmm3,0x1f
 psrad  xmm2,0x1f
 pxor   xmm0,xmm3
 pxor   xmm1,xmm2
 psubd  xmm1,xmm2
 psubd  xmm0,xmm3
 pshufd xmm2,xmm2,0x50
 pshufd xmm1,xmm1,0x50
 pshufd xmm0,xmm0,0x50
 pshufd xmm3,xmm3,0x50
 pmuludq xmm0,xmm1
 pxor   xmm2,xmm3
 pxor   xmm0,xmm2
 psubq  xmm0,xmm2
 ret    
 nop    WORD PTR [rax+rax*1+0x0]
multiply32_high_s(int __vector(4),0xfa
 pshufd xmm1,0xfa
 pshufd xmm0,0xfa
 pshufd xmm3,0xfa
 pmuludq xmm0,xmm2
 ret    
 nop    WORD PTR [rax+rax*1+0x0]