SIMD assinado com multiplicação sem sinal para 64 bits * 64 bits para 128 bits

Eu criei uma function que faz 64 bits * 64 bits para 128 bits usando SIMD. Atualmente eu o implementei usando o SSE2 (acertally SSE4.1). Isso significa que ele faz dois produtos 64b * 64b a 128b ao mesmo tempo. A mesma ideia poderia ser estendida para AVX2 ou AVX512, fornecendo quatro ou oito produtos 64b * 64 a 128b ao mesmo tempo. Eu baseei meu algoritmo em http://www.hackersdelight.org/hdcodetxt/muldws.c.txt

Esse algoritmo faz uma multiplicação sem sinal, uma multiplicação com sinal e duas multiplicações com sinal *. As assinadas * assinadas e não assinadas * operações não assinadas são fáceis de fazer usando _mm_mul_epi32 e _mm_mul_epu32 . Mas os produtos misturados assinados e não assinados me causaram problemas. Considere por exemplo.

 int32_t x = 0x80000000; uint32_t y = 0x7fffffff; int64_t z = (int64_t)x*y; 

O produto de palavra dupla deve ser 0xc000000080000000 . Mas como você pode obter isso se você assumir que seu compilador sabe como lidar com tipos mistos? Isto é o que eu criei:

 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 

Usando o SSE isso pode ser feito assim

 __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 

Isso dá o resultado correto. Mas eu tenho que fazer isso duas vezes (uma vez ao quadrado) e agora é uma fração significativa da minha function. Existe uma maneira mais eficiente de fazer isso com SSE4.2, AVX2 (quatro produtos de 128 bits) ou até mesmo AVX512 (oito produtos de 128 bits)?

Talvez existam maneiras mais eficientes de fazer isso do que com o SIMD? É um monte de cálculos para obter a palavra superior.

Editar: com base no comentário de @ElderBug parece que a maneira de fazer isso não é com SIMD, mas com a instrução mul . Por que vale a pena, no caso de alguém querer ver o quão complicado isso é, aqui está a function de trabalho completo (eu só tenho que trabalhar, então eu não tenho otimizado isso, mas eu não acho que vale a pena).

 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; } 

Fica pior. Não há instrução / instrução _mm_srai_epi64 até o AVX512, então tive que fazer o meu próprio.

 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); } 

Minha implementação de _mm_srai_epi64 acima está incompleta. Eu acho que eu estava usando o Vector Class Library da Agner Fog. Se você olhar no arquivo vectori128.h você encontra

 static inline Vec2q operator >> (Vec2q const & a, int32_t b) { // instruction does not exist. Split into 32-bit shifts if (b > 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); } } 

O caminho certo para pensar sobre os limites de throughput da multiplicação de inteiros usando várias instruções é em termos de quantos “bits de produto” você pode calcular por ciclo.

mulx produz um resultado de 64×64 -> 128 a cada ciclo; isto é 64×64 = 4096 “bits de produto por ciclo”

Se você juntar um multiplicador no SIMD com instruções que fazem multiplex de 32×32 a> 64 bits, será necessário obter quatro resultados a cada ciclo para corresponder ao mulx (4x32x32 = 4096). Se não houvesse aritmética além das multiplicações, você acabaria por quebrar mesmo no AVX2. Infelizmente, como você notou, há muita aritmética além das multiplicações, então isso é um total não-inicial no hardware da geração atual.

Eu encontrei uma solução SIMD que é muito mais simples e não precisa de produtos signed*unsigned . Não estou mais convencido de que o SIMD (pelo menos com AVX2 e AV512) não pode competir com o mulx . Em alguns casos, o SIMD pode competir com o mulx . O único caso em que estou ciente é na multiplicação baseada em FFT de grandes números .

O truque era fazer a multiplicação sem sinal primeiro e depois corrigir. Eu aprendi como fazer isso a partir desta resposta 32-bit-assinado-multiplicação-sem-usar-64-bit-tipo de dados . A correção é simples para (hi,lo) = x*y fazer multiplicação sem sinal primeiro e depois corrigir hi assim:

 hi -= ((x<0) ? y : 0) + ((y<0) ? x : 0) 

Isso pode ser feito com o intrínseco do 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); } 

O código para a multiplicação não assinada é mais simples, pois não precisa de produtos signed*unsigned mistos. Além disso, uma vez que não está assinado, não é necessário o direito de deslocamento aritmético, que possui apenas uma instrução para o AVX512. Na verdade, a seguinte function só precisa do 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; } 

Isso usa

 4x mul_epu32 5x add_epi64 2x shuffle_epi32 2x and 2x srli_epi64 1x slli_epi64 **************** 16 instructions 

AVX512 tem o _mm_mullo_epi64 intrínseco que pode calcular lo com uma instrução. Nesse caso, a alternativa pode ser usada (comente as linhas com o comentário (*) e descomente a linha alternativa):

 5x mul_epu32 4x add_epi64 2x shuffle_epi32 1x and 2x srli_epi64 **************** 14 instructions 

Para alterar o código para largura total AVX2, substitua _mm por _mm256 , si128 por si256 e __m128i por __m256i por AVX512 substitua-os por _mm512 , si512 e __m512i .