SIMD signed with unsigned multiplication for 64-bit * 64-bit to 128-bit
SIMD signed with unsigned multiplication for 64-bit * 64-bit to 128-bit
我创建了一个使用 SIMD 执行 64 位 * 64 位到 128 位的函数。目前我已经使用 SSE2（实际上是 SSE4.1）实现了它。这意味着它同时处理两个 64b*64b 到 128b 的产品。同样的想法可以扩展到 AVX2 或 AVX512，同时提供四个或八个 64b*64 到 128b 产品。
我的算法基于 http://www.hackersdelight.org/hdcodetxt/muldws.c.txt
该算法执行一次无符号乘法、一次带符号乘法和两次带符号 * 无符号乘法。 signed * signed 和 unsigned * unsigned 操作很容易使用 _mm_mul_epi32 和 _mm_mul_epu32 来完成。但是混合签名和未签名的产品给我带来了麻烦。
举个例子。
int32_t x = 0x80000000;
uint32_t y = 0x7fffffff;
int64_t z = (int64_t)x*y;
双字积应该是0xc000000080000000。但是如果你假设你的编译器确实知道如何处理混合类型，你怎么能得到这个呢？这是我想出的：
int64_t sign = x<0; sign*=-1; //get the sign and make it all ones
uint32_t t = abs(x); //if x<0 take two's complement again
uint64_t prod = (uint64_t)t*y; //unsigned product
int64_t z = (prod ^ sign) - sign; //take two's complement based on the sign
使用 SSE 可以像这样完成
__m128i xh; //(xl2, xh2, xl1, xh1) high is signed, low unsigned
__m128i yl; //(yh2, yl2, yh2, yl2)
__m128i xs = _mm_cmpgt_epi32(_mm_setzero_si128(), xh); // get sign
xs = _mm_shuffle_epi32(xs, 0xA0); // extend sign
__m128i t = _mm_sign_epi32(xh,xh); // abs(xh)
__m128i prod = _mm_mul_epu32(t, yl); // unsigned (xh2*yl2,xh1*yl1)
__m128i inv = _mm_xor_si128(prod,xs); // invert bits if negative
__m128i z = _mm_sub_epi64(inv,xs); // add 1 if negative
这给出了正确的结果。但是我必须这样做两次（平方时一次），它现在是我功能的重要部分。使用 SSE4.2、AVX2（四个 128 位产品）甚至 AVX512（八个 128 位产品）是否有更有效的方法？
也许有比使用 SIMD 更有效的方法来做到这一点？得到上位词需要大量计算
编辑：根据@ElderBug 的评论，看起来这样做的方法不是使用 SIMD，而是使用 mul 指令。对于它的价值，如果有人想看看它有多复杂，这里是完整的工作功能（我刚刚开始工作，所以我没有优化它，但我认为它不值得）。
void muldws1_sse(__m128i x, __m128i y, __m128i *lo, __m128i *hi) {
__m128i lomask = _mm_set1_epi64x(0xffffffff);
__m128i xh = _mm_shuffle_epi32(x, 0xB1); // x0l, x0h, x1l, x1h
__m128i yh = _mm_shuffle_epi32(y, 0xB1); // y0l, y0h, y1l, y1h
__m128i xs = _mm_cmpgt_epi32(_mm_setzero_si128(), xh);
__m128i ys = _mm_cmpgt_epi32(_mm_setzero_si128(), yh);
xs = _mm_shuffle_epi32(xs, 0xA0);
ys = _mm_shuffle_epi32(ys, 0xA0);
__m128i w0 = _mm_mul_epu32(x, y); // x0l*y0l, y0l*y0h
__m128i w3 = _mm_mul_epi32(xh, yh); // x0h*y0h, x1h*y1h
xh = _mm_sign_epi32(xh,xh);
yh = _mm_sign_epi32(yh,yh);
__m128i w1 = _mm_mul_epu32(x, yh); // x0l*y0h, x1l*y1h
__m128i w2 = _mm_mul_epu32(xh, y); // x0h*y0l, x1h*y0l
__m128i yinv = _mm_xor_si128(w1,ys); // invert bits if negative
w1 = _mm_sub_epi64(yinv,ys); // add 1
__m128i xinv = _mm_xor_si128(w2,xs); // invert bits if negative
w2 = _mm_sub_epi64(xinv,xs); // add 1
__m128i w0l = _mm_and_si128(w0, lomask);
__m128i w0h = _mm_srli_epi64(w0, 32);
__m128i s1 = _mm_add_epi64(w1, w0h); // xl*yh + w0h;
__m128i s1l = _mm_and_si128(s1, lomask); // lo(wl*yh + w0h);
__m128i s1h = _mm_srai_epi64(s1, 32);
__m128i s2 = _mm_add_epi64(w2, s1l); //xh*yl + s1l
__m128i s2l = _mm_slli_epi64(s2, 32);
__m128i s2h = _mm_srai_epi64(s2, 32); //arithmetic shift right
__m128i hi1 = _mm_add_epi64(w3, s1h);
hi1 = _mm_add_epi64(hi1, s2h);
__m128i lo1 = _mm_add_epi64(w0l, s2l);
*hi = hi1;
*lo = lo1;
}
情况变得更糟。在 AVX512 之前没有 _mm_srai_epi64 instrinsic/instruction 所以我不得不自己制作。
static inline __m128i _mm_srai_epi64(__m128i a, int b) {
__m128i sra = _mm_srai_epi32(a,32);
__m128i srl = _mm_srli_epi64(a,32);
__m128i mask = _mm_set_epi32(-1,0,-1,0);
__m128i out = _mm_blendv_epi8(srl, sra, mask);
}
我上面 _mm_srai_epi64 的实现不完整。我想我使用的是 Agner Fog 的 Vector Class Library。如果您查看文件 vectori128.h，您会发现
static inline Vec2q operator >> (Vec2q const & a, int32_t b) {
// instruction does not exist. Split into 32-bit shifts
if (b <= 32) {
__m128i bb = _mm_cvtsi32_si128(b); // b
__m128i sra = _mm_sra_epi32(a,bb); // a >> b signed dwords
__m128i srl = _mm_srl_epi64(a,bb); // a >> b unsigned qwords
__m128i mask = _mm_setr_epi32(0,-1,0,-1); // mask for signed high part
return selectb(mask,sra,srl);
}
else { // b > 32
__m128i bm32 = _mm_cvtsi32_si128(b-32); // b - 32
__m128i sign = _mm_srai_epi32(a,31); // sign of a
__m128i sra2 = _mm_sra_epi32(a,bm32); // a >> (b-32) signed dwords
__m128i sra3 = _mm_srli_epi64(sra2,32); // a >> (b-32) >> 32 (second shift unsigned qword)
__m128i mask = _mm_setr_epi32(0,-1,0,-1); // mask for high part containing only sign
return selectb(mask,sign,sra3);
}
}
考虑使用各种指令进行整数乘法的吞吐量限制的正确方法是根据每个周期可以计算多少"product bits"。
mulx 每个周期产生一个 64x64 -> 128 的结果；那是 64x64 = 4096 "product bits per cycle"
如果您在 SIMD 上从执行 32x32 -> 64 位乘法的指令中拼凑出一个乘法器，则每个周期需要能够获得四个结果才能匹配 mulx (4x32x32 = 4096)。如果除了乘法之外没有其他算术，您将在 AVX2 上实现收支平衡。不幸的是，正如您所注意到的，除了乘法运算之外还有很多算术运算，因此这在当前这一代硬件上完全无法启动。
我找到了一个简单得多且不需要 signed*unsigned 产品的 SIMD 解决方案。 我不再相信 SIMD（至少对于 AVX2 和 AV512）无法与 mulx 竞争。 在某些情况下 SIMD 可以与 mulx 竞争。我知道的唯一情况是 FFT based multiplication of large numbers.
诀窍是先做无符号乘法，然后再校正。我从这个答案 32-bit-signed-multiplication-without-using-64-bit-data-type 中学会了如何做到这一点。更正很简单 (hi,lo) = x*y 先做无符号乘法然后更正 hi 像这样：
hi -= ((x<0) ? y : 0) + ((y<0) ? x : 0)
这可以通过 SSE4.2 内部函数完成 _mm_cmpgt_epi64
void muldws1_sse(__m128i x, __m128i y, __m128i *lo, __m128i *hi) {
muldwu1_sse(x,y,lo,hi);
//hi -= ((x<0) ? y : 0) + ((y<0) ? x : 0);
__m128i xs = _mm_cmpgt_epi64(_mm_setzero_si128(), x);
__m128i ys = _mm_cmpgt_epi64(_mm_setzero_si128(), y);
__m128i t1 = _mm_and_si128(y,xs);
__m128i t2 = _mm_and_si128(x,ys);
*hi = _mm_sub_epi64(*hi,t1);
*hi = _mm_sub_epi64(*hi,t2);
}
无符号乘法的代码更简单，因为它不需要混合 signed*unsigned 乘积。此外，由于它是无符号的，因此不需要只有 AVX512 指令的算术右移。实际上下面的函数只需要SSE2:
void muldwu1_sse(__m128i x, __m128i y, __m128i *lo, __m128i *hi) {
__m128i lomask = _mm_set1_epi64x(0xffffffff);
__m128i xh = _mm_shuffle_epi32(x, 0xB1); // x0l, x0h, x1l, x1h
__m128i yh = _mm_shuffle_epi32(y, 0xB1); // y0l, y0h, y1l, y1h
__m128i w0 = _mm_mul_epu32(x, y); // x0l*y0l, x1l*y1l
__m128i w1 = _mm_mul_epu32(x, yh); // x0l*y0h, x1l*y1h
__m128i w2 = _mm_mul_epu32(xh, y); // x0h*y0l, x1h*y0l
__m128i w3 = _mm_mul_epu32(xh, yh); // x0h*y0h, x1h*y1h
__m128i w0l = _mm_and_si128(w0, lomask); //(*)
__m128i w0h = _mm_srli_epi64(w0, 32);
__m128i s1 = _mm_add_epi64(w1, w0h);
__m128i s1l = _mm_and_si128(s1, lomask);
__m128i s1h = _mm_srli_epi64(s1, 32);
__m128i s2 = _mm_add_epi64(w2, s1l);
__m128i s2l = _mm_slli_epi64(s2, 32); //(*)
__m128i s2h = _mm_srli_epi64(s2, 32);
__m128i hi1 = _mm_add_epi64(w3, s1h);
hi1 = _mm_add_epi64(hi1, s2h);
__m128i lo1 = _mm_add_epi64(w0l, s2l); //(*)
//__m128i lo1 = _mm_mullo_epi64(x,y); //alternative
*hi = hi1;
*lo = lo1;
}
这使用
4x mul_epu32
5x add_epi64
2x shuffle_epi32
2x and
2x srli_epi64
1x slli_epi64
****************
16 instructions
AVX512 具有 _mm_mullo_epi64 内在函数，可以用一条指令计算 lo。在这种情况下，可以使用替代方法（用 (*) 注释注释行并取消注释替代行）：
5x mul_epu32
4x add_epi64
2x shuffle_epi32
1x and
2x srli_epi64
****************
14 instructions
要更改全角 AVX2 的代码，请将 _mm 替换为 _mm256，将 si128 替换为 si256，将 __m128i 替换为 __m256i对于 AVX512，将它们替换为 _mm512、si512 和 __m512i。
我创建了一个使用 SIMD 执行 64 位 * 64 位到 128 位的函数。目前我已经使用 SSE2（实际上是 SSE4.1）实现了它。这意味着它同时处理两个 64b*64b 到 128b 的产品。同样的想法可以扩展到 AVX2 或 AVX512，同时提供四个或八个 64b*64 到 128b 产品。 我的算法基于 http://www.hackersdelight.org/hdcodetxt/muldws.c.txt
该算法执行一次无符号乘法、一次带符号乘法和两次带符号 * 无符号乘法。 signed * signed 和 unsigned * unsigned 操作很容易使用 _mm_mul_epi32 和 _mm_mul_epu32 来完成。但是混合签名和未签名的产品给我带来了麻烦。
举个例子。
int32_t x = 0x80000000;
uint32_t y = 0x7fffffff;
int64_t z = (int64_t)x*y;
双字积应该是0xc000000080000000。但是如果你假设你的编译器确实知道如何处理混合类型，你怎么能得到这个呢？这是我想出的：
int64_t sign = x<0; sign*=-1; //get the sign and make it all ones
uint32_t t = abs(x); //if x<0 take two's complement again
uint64_t prod = (uint64_t)t*y; //unsigned product
int64_t z = (prod ^ sign) - sign; //take two's complement based on the sign
使用 SSE 可以像这样完成
__m128i xh; //(xl2, xh2, xl1, xh1) high is signed, low unsigned
__m128i yl; //(yh2, yl2, yh2, yl2)
__m128i xs = _mm_cmpgt_epi32(_mm_setzero_si128(), xh); // get sign
xs = _mm_shuffle_epi32(xs, 0xA0); // extend sign
__m128i t = _mm_sign_epi32(xh,xh); // abs(xh)
__m128i prod = _mm_mul_epu32(t, yl); // unsigned (xh2*yl2,xh1*yl1)
__m128i inv = _mm_xor_si128(prod,xs); // invert bits if negative
__m128i z = _mm_sub_epi64(inv,xs); // add 1 if negative
这给出了正确的结果。但是我必须这样做两次（平方时一次），它现在是我功能的重要部分。使用 SSE4.2、AVX2（四个 128 位产品）甚至 AVX512（八个 128 位产品）是否有更有效的方法？
也许有比使用 SIMD 更有效的方法来做到这一点？得到上位词需要大量计算
编辑：根据@ElderBug 的评论，看起来这样做的方法不是使用 SIMD，而是使用 mul 指令。对于它的价值，如果有人想看看它有多复杂，这里是完整的工作功能（我刚刚开始工作，所以我没有优化它，但我认为它不值得）。
void muldws1_sse(__m128i x, __m128i y, __m128i *lo, __m128i *hi) {
__m128i lomask = _mm_set1_epi64x(0xffffffff);
__m128i xh = _mm_shuffle_epi32(x, 0xB1); // x0l, x0h, x1l, x1h
__m128i yh = _mm_shuffle_epi32(y, 0xB1); // y0l, y0h, y1l, y1h
__m128i xs = _mm_cmpgt_epi32(_mm_setzero_si128(), xh);
__m128i ys = _mm_cmpgt_epi32(_mm_setzero_si128(), yh);
xs = _mm_shuffle_epi32(xs, 0xA0);
ys = _mm_shuffle_epi32(ys, 0xA0);
__m128i w0 = _mm_mul_epu32(x, y); // x0l*y0l, y0l*y0h
__m128i w3 = _mm_mul_epi32(xh, yh); // x0h*y0h, x1h*y1h
xh = _mm_sign_epi32(xh,xh);
yh = _mm_sign_epi32(yh,yh);
__m128i w1 = _mm_mul_epu32(x, yh); // x0l*y0h, x1l*y1h
__m128i w2 = _mm_mul_epu32(xh, y); // x0h*y0l, x1h*y0l
__m128i yinv = _mm_xor_si128(w1,ys); // invert bits if negative
w1 = _mm_sub_epi64(yinv,ys); // add 1
__m128i xinv = _mm_xor_si128(w2,xs); // invert bits if negative
w2 = _mm_sub_epi64(xinv,xs); // add 1
__m128i w0l = _mm_and_si128(w0, lomask);
__m128i w0h = _mm_srli_epi64(w0, 32);
__m128i s1 = _mm_add_epi64(w1, w0h); // xl*yh + w0h;
__m128i s1l = _mm_and_si128(s1, lomask); // lo(wl*yh + w0h);
__m128i s1h = _mm_srai_epi64(s1, 32);
__m128i s2 = _mm_add_epi64(w2, s1l); //xh*yl + s1l
__m128i s2l = _mm_slli_epi64(s2, 32);
__m128i s2h = _mm_srai_epi64(s2, 32); //arithmetic shift right
__m128i hi1 = _mm_add_epi64(w3, s1h);
hi1 = _mm_add_epi64(hi1, s2h);
__m128i lo1 = _mm_add_epi64(w0l, s2l);
*hi = hi1;
*lo = lo1;
}
情况变得更糟。在 AVX512 之前没有 _mm_srai_epi64 instrinsic/instruction 所以我不得不自己制作。
static inline __m128i _mm_srai_epi64(__m128i a, int b) {
__m128i sra = _mm_srai_epi32(a,32);
__m128i srl = _mm_srli_epi64(a,32);
__m128i mask = _mm_set_epi32(-1,0,-1,0);
__m128i out = _mm_blendv_epi8(srl, sra, mask);
}
我上面 _mm_srai_epi64 的实现不完整。我想我使用的是 Agner Fog 的 Vector Class Library。如果您查看文件 vectori128.h，您会发现
static inline Vec2q operator >> (Vec2q const & a, int32_t b) {
// instruction does not exist. Split into 32-bit shifts
if (b <= 32) {
__m128i bb = _mm_cvtsi32_si128(b); // b
__m128i sra = _mm_sra_epi32(a,bb); // a >> b signed dwords
__m128i srl = _mm_srl_epi64(a,bb); // a >> b unsigned qwords
__m128i mask = _mm_setr_epi32(0,-1,0,-1); // mask for signed high part
return selectb(mask,sra,srl);
}
else { // b > 32
__m128i bm32 = _mm_cvtsi32_si128(b-32); // b - 32
__m128i sign = _mm_srai_epi32(a,31); // sign of a
__m128i sra2 = _mm_sra_epi32(a,bm32); // a >> (b-32) signed dwords
__m128i sra3 = _mm_srli_epi64(sra2,32); // a >> (b-32) >> 32 (second shift unsigned qword)
__m128i mask = _mm_setr_epi32(0,-1,0,-1); // mask for high part containing only sign
return selectb(mask,sign,sra3);
}
}
考虑使用各种指令进行整数乘法的吞吐量限制的正确方法是根据每个周期可以计算多少"product bits"。
mulx 每个周期产生一个 64x64 -> 128 的结果；那是 64x64 = 4096 "product bits per cycle"
如果您在 SIMD 上从执行 32x32 -> 64 位乘法的指令中拼凑出一个乘法器，则每个周期需要能够获得四个结果才能匹配 mulx (4x32x32 = 4096)。如果除了乘法之外没有其他算术，您将在 AVX2 上实现收支平衡。不幸的是，正如您所注意到的，除了乘法运算之外还有很多算术运算，因此这在当前这一代硬件上完全无法启动。
我找到了一个简单得多且不需要 signed*unsigned 产品的 SIMD 解决方案。 我不再相信 SIMD（至少对于 AVX2 和 AV512）无法与 在某些情况下 SIMD 可以与 mulx 竞争。mulx 竞争。我知道的唯一情况是 FFT based multiplication of large numbers.
诀窍是先做无符号乘法，然后再校正。我从这个答案 32-bit-signed-multiplication-without-using-64-bit-data-type 中学会了如何做到这一点。更正很简单 (hi,lo) = x*y 先做无符号乘法然后更正 hi 像这样：
hi -= ((x<0) ? y : 0) + ((y<0) ? x : 0)
这可以通过 SSE4.2 内部函数完成 _mm_cmpgt_epi64
void muldws1_sse(__m128i x, __m128i y, __m128i *lo, __m128i *hi) {
muldwu1_sse(x,y,lo,hi);
//hi -= ((x<0) ? y : 0) + ((y<0) ? x : 0);
__m128i xs = _mm_cmpgt_epi64(_mm_setzero_si128(), x);
__m128i ys = _mm_cmpgt_epi64(_mm_setzero_si128(), y);
__m128i t1 = _mm_and_si128(y,xs);
__m128i t2 = _mm_and_si128(x,ys);
*hi = _mm_sub_epi64(*hi,t1);
*hi = _mm_sub_epi64(*hi,t2);
}
无符号乘法的代码更简单，因为它不需要混合 signed*unsigned 乘积。此外，由于它是无符号的，因此不需要只有 AVX512 指令的算术右移。实际上下面的函数只需要SSE2:
void muldwu1_sse(__m128i x, __m128i y, __m128i *lo, __m128i *hi) {
__m128i lomask = _mm_set1_epi64x(0xffffffff);
__m128i xh = _mm_shuffle_epi32(x, 0xB1); // x0l, x0h, x1l, x1h
__m128i yh = _mm_shuffle_epi32(y, 0xB1); // y0l, y0h, y1l, y1h
__m128i w0 = _mm_mul_epu32(x, y); // x0l*y0l, x1l*y1l
__m128i w1 = _mm_mul_epu32(x, yh); // x0l*y0h, x1l*y1h
__m128i w2 = _mm_mul_epu32(xh, y); // x0h*y0l, x1h*y0l
__m128i w3 = _mm_mul_epu32(xh, yh); // x0h*y0h, x1h*y1h
__m128i w0l = _mm_and_si128(w0, lomask); //(*)
__m128i w0h = _mm_srli_epi64(w0, 32);
__m128i s1 = _mm_add_epi64(w1, w0h);
__m128i s1l = _mm_and_si128(s1, lomask);
__m128i s1h = _mm_srli_epi64(s1, 32);
__m128i s2 = _mm_add_epi64(w2, s1l);
__m128i s2l = _mm_slli_epi64(s2, 32); //(*)
__m128i s2h = _mm_srli_epi64(s2, 32);
__m128i hi1 = _mm_add_epi64(w3, s1h);
hi1 = _mm_add_epi64(hi1, s2h);
__m128i lo1 = _mm_add_epi64(w0l, s2l); //(*)
//__m128i lo1 = _mm_mullo_epi64(x,y); //alternative
*hi = hi1;
*lo = lo1;
}
这使用
4x mul_epu32
5x add_epi64
2x shuffle_epi32
2x and
2x srli_epi64
1x slli_epi64
****************
16 instructions
AVX512 具有 _mm_mullo_epi64 内在函数，可以用一条指令计算 lo。在这种情况下，可以使用替代方法（用 (*) 注释注释行并取消注释替代行）：
5x mul_epu32
4x add_epi64
2x shuffle_epi32
1x and
2x srli_epi64
****************
14 instructions
要更改全角 AVX2 的代码，请将 _mm 替换为 _mm256，将 si128 替换为 si256，将 __m128i 替换为 __m256i对于 AVX512，将它们替换为 _mm512、si512 和 __m512i。