你如何在 SSE2 上进行带符号的 32 位扩展乘法？

Posted 2023-02-16

【中文标题】你如何在 SSE2 上进行带符号的 32 位扩展乘法？

【英文标题】：How do you do signed 32bit widening multiplication on SSE2?

【发布时间】：2020-10-09 02:32:24

【问题描述】：

在查看 WebAssembly SIMD 扩展乘法提案时提出了这个问题。

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

（非扩展相对容易；随机输入 2x pmuludq 并取 4 个结果的下半部分来模拟 SSE4.1 pmulld）。

【问题讨论】：

请注意，SSE4.1 添加了 pmuldq 扩大了 signed 乘法。所以它存在，但你是正确的，它不是 SSE2 的一部分，因此不是 x86-64 的基线：/

【参考方案1】：

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, 1, 0));
    b = _mm_shuffle_epi32(b, _MM_SHUFFLE(3, 1, 1, 0));
    __m128i unsignedProduct = _mm_mul_epu32(a, b);
    __m128i threshold = _mm_set_epi32(INT_MIN, 0, 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, 0, 1));
    __m128i y = _mm_shuffle_epi32(_mm_and_si128(signB, a), _MM_SHUFFLE(2, 3, 0, 1));
    return _mm_sub_epi32(_mm_sub_epi32(unsignedProduct, x), y);

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

【讨论】：

我认为我们的要求是它以与原始相同的顺序输出。 pmuludq 做同样的事情，甚至你正在谈论的东西。我的版本中的洗牌是为了适应标准的要求而存在的。如果您为自己添加随机播放，它的表现会更好还是不同？

@DanWeber 它仍然减少了两条指令，但我没有对任何东西进行基准测试

如果你看一下 Clang 的汇编版本（它在神螺栓上），它设法减少了我的 1 个 shuffle op，尽管我不完全确定我理解它做了什么。

threshold 可以在 2 条指令中即时生成：pcmpeqd xmm0, xmm0（全为）/psllq xmm0, 63 仅在每个 qword 一半中设置符号位。编译器仍然只会加载一个常量，但是如果您正在 JITing 并且想要避免加载，您可以这样做，@DanWeber。如果您想从循环中提升出来，那就太好了，否则可能会加载一个常量，而不是每次都花费更多的向量 ALU 指令来重新生成它。

你能用psrad算术右移得到a<0掩码吗？

【参考方案2】：

在##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,1,0,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,3,2,2);
    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;



int main(int argc, char* argv[]) 
    int32x4_t a-5,500,-5000,50000;
    int32x4_t b10,-100,-5000,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), 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,0xfa
 pshufd xmm1,xmm1,0xfa
 pshufd xmm0,xmm0,0xfa
 pshufd xmm3,xmm3,0xfa
 pmuludq xmm0,xmm1
 pxor   xmm2,xmm3
 pxor   xmm0,xmm2
 psubq  xmm0,xmm2
 ret    
 nop    WORD PTR [rax+rax*1+0x0]

【讨论】：

这值得吗？标量后备看起来相对有希望..

这是个好问题。我很想听听有关它的反馈。

@harold：可能做这么多向量工作的唯一原因（即使使用您的优化版本）是如果将结果放在向量中真的很有用，并且输入以向量 regs 开始。即作为 SIMD 操作之间的处理步骤，而不是单独的。也许那时还没有； godbolt.org/z/MWY8hq 显示它有 16 条指令，对于 SSE4.1 来说太旧的 CPU 没有 mov-elimination of xor-elimination。 K10 和 Core2 至少有 128 位宽的向量 ALU；早期的 CPU 没有，这使得这些向量 insn 更加昂贵。

所以，是的，movq/shuffle/movq 或存储/加载以提供一些标量 imul r64,r64（或者在具有慢速 64 位整数乘法的 CPU 上可能是 imul r32）可能很好，即使我们有吃一个商店转发的摊位以重新进入向量。 （如果我们使用 32 位 imul，则 4x 标量存储 + 矢量加载可能比 4x movd + 3x punpck 具有更好的吞吐量，因此结果被分成 4 个 regs。）

@PeterCordes 我很想知道为什么我们支持 sse2 的 wasm。我的猜测是 v8 在适用时需要标量变体的程序集。