SIMD firmado con multiplicación sin signo para 64 bits * 64 bits a 128 bits

Creé una función que usa 64 bits * 64 bits a 128 bits usando SIMD. Actualmente lo he implementado usando SSE2 (acutalmente SSE4.1). Esto significa que tiene dos productos de 64b * 64b a 128b al mismo tiempo. La misma idea podría extenderse a AVX2 o AVX512 dando cuatro u ocho productos de 64b * 64 a 128b al mismo tiempo. Basé mi algoritmo en http://www.hackersdelight.org/hdcodetxt/muldws.c.txt

Ese algoritmo realiza una multiplicación sin signo, una multiplicación firmada y dos multiplicaciones firmadas * sin signo. Las operaciones * firmadas y sin firmar * firmadas son fáciles de hacer usando _mm_mul_epi32 y _mm_mul_epu32 . Pero los productos mixtos, firmados y no firmados me causaron problemas. Considera por ejemplo.

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

El producto de doble palabra debe ser 0xc000000080000000 . ¿Pero cómo puede obtener esto si supone que su comstackdor sabe cómo manejar tipos mixtos? Esto es lo que se me ocurrió:

 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 SSE esto se puede hacer así

 __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 

Esto da el resultado correcto. Pero tengo que hacer esto dos veces (una al cuadrar) y ahora es una fracción significativa de mi función. ¿Hay alguna forma más eficiente de hacerlo con SSE4.2, AVX2 (cuatro productos de 128 bits) o incluso AVX512 (ocho productos de 128 bits)?

Tal vez hay formas más eficientes de hacer esto que con SIMD? Son muchos cálculos para obtener la palabra superior.

Editar: basado en el comentario de @ElderBug parece que la forma de hacerlo no es con SIMD sino con la instrucción mul . Por lo que vale, en caso de que alguien quiera ver lo complicado que es esto, aquí está la función de trabajo completa (lo hice funcionar, así que no lo he optimizado, pero no creo que valga la 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; } 

Se pone peor. No hay _mm_srai_epi64 / instruction hasta el AVX512, así que tuve que hacer el mío.

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

Mi implementación de _mm_srai_epi64 anterior está incompleta. Creo que estaba usando Agner Fog’s Vector Class Library . Si miras en el archivo vectori128.h encuentras

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

La manera correcta de pensar en los límites de rendimiento de la multiplicación de enteros usando varias instrucciones es en términos de cuántos “bits de producto” puede calcular por ciclo.

mulx produce un resultado de 64×64 -> 128 cada ciclo; eso es 64×64 = 4096 “bits de producto por ciclo”

Si junta un multiplicador en el SIMD sin instrucciones que multipliquen 32×32 -> 64 bits, debe poder obtener cuatro resultados en cada ciclo para que coincida con mulx (4x32x32 = 4096). Si no hubiera más aritmética que las multiplicaciones, solo alcanzaría el punto de equilibrio en AVX2. Desafortunadamente, como habrás notado, hay mucha aritmética además de las multiplicaciones, por lo que este es un total no arrancador en el hardware de generación actual.

Encontré una solución SIMD que es mucho más simple y no necesita productos signed*unsigned . Ya no estoy convencido de que SIMD (al menos con AVX2 y AV512) no pueda competir con mulx . En algunos casos SIMD puede competir con mulx . El único caso que conozco es en la multiplicación basada en FFT de números grandes .

El truco consistía en hacer una multiplicación sin firmar primero y luego corregir. Aprendí cómo hacer esto a partir de esta respuesta 32-bit-signed-multiplication-without-using-64-bit-data-type . La corrección es simple para (hi,lo) = x*y hacer la multiplicación sin signo primero y luego corregirla de esta manera:

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

Esto se puede hacer con el SSE4.2 intrínseco _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); } 

El código para la multiplicación sin firmar es más simple ya que no necesita productos mixtos signed*unsigned . Además, como no tiene firma, no necesita desplazamiento aritmético, que solo tiene una instrucción para AVX512. De hecho, la siguiente función solo necesita 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; } 

Esto usa

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

AVX512 tiene el _mm_mullo_epi64 intrínseco que puede calcular lo con una instrucción. En este caso, se puede usar la alternativa (comente las líneas con el comentario (*) y elimine el comentario de la línea alternativa):

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

Para cambiar el código para AVX2 de ancho completo, reemplace _mm con _mm256 , si128 con si256 , y __m128i con __m256i para AVX512 reemplácelos con _mm512 , si512 y __m512i .