SIMD signé avec multiplication non signée pour 64 bits * 64 bits à 128 bits

J’ai créé une fonction qui exécute SIMD de 64 bits * 64 bits à 128 bits. Actuellement, je l’ai implémenté en utilisant SSE2 (en réalité SSE4.1). Cela signifie qu’il fabrique deux produits 64b * 64b à 128b en même temps. La même idée pourrait être étendue à AVX2 ou AVX512 en donnant quatre ou huit produits 64b * 64 à 128b en même temps. J’ai basé mon algorithme sur http://www.hackersdelight.org/hdcodetxt/muldws.c.txt

Cet algorithme effectue une multiplication non signée, une multiplication signée et deux multiplications signées * non signées. Les opérations signées * signées et non signées * sont faciles à effectuer avec _mm_mul_epi32 et _mm_mul_epu32 . Mais les produits mixtes signés et non signés m’ont causé des problèmes. Considérons par exemple.

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

Le produit de mot double doit être 0xc000000080000000 . Mais comment pouvez-vous l’obtenir si vous supposez que votre compilateur sait comment gérer des types mixtes? Voici ce que je suis venu avec:

 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 

En utilisant SSE cela peut être fait comme ceci

 __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 

Cela donne le résultat correct. Mais je dois le faire deux fois (une fois au carré) et c’est maintenant une fraction importante de ma fonction. Existe-t-il un moyen plus efficace de le faire avec SSE4.2, AVX2 (quatre produits 128 bits) ou même AVX512 (huit produits 128 bits)?

Peut-être existe-t-il des moyens plus efficaces de le faire qu’avec SIMD? C’est beaucoup de calculs pour obtenir le mot supérieur.

Edit: basé sur le commentaire de @ElderBug, il semble que la manière de procéder ne soit pas avec SIMD mais avec l’instruction mul . Pour ce que ça vaut, au cas où quelqu’un voudrait voir à quel point c’est compliqué, voici la fonction qui fonctionne (je viens de le faire fonctionner, donc je ne l’ai pas optimisé mais je ne pense pas que ça en vaut la peine).

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

Ça a empiré. Il n’y a pas d’instructions / d’instructions _mm_srai_epi64 avant AVX512, je devais donc créer le mien.

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

Mon implémentation de _mm_srai_epi64 ci-dessus est incomplète. Je pense que j’utilisais la bibliothèque de classes vectorielles d’ Agner Fog. Si vous regardez dans le fichier vectori128.h vous trouvez

 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 bonne façon de penser aux limites de débit de la multiplication d’entiers à l’aide de diverses instructions est de définir le nombre de «produits» que vous pouvez calculer par cycle.

mulx produit un résultat 64×64 -> 128 à chaque cycle; c’est 64×64 = 4096 “bits de produit par cycle”

Si vous assemblez un multiplicateur sur SIMD à partir d’instructions faisant des multiplications 32×32 -> 64 bits, vous devez pouvoir obtenir quatre résultats à chaque cycle pour correspondre à mulx (4x32x32 = 4096). S’il n’y avait pas d’arithmétique autre que les multiplications, vous auriez juste à équilibrer sur AVX2. Malheureusement, comme vous l’avez remarqué, il y a beaucoup d’arithmétiques autres que les multiplications, il s’agit donc d’un non-démarreur total sur le matériel de la génération actuelle.

J’ai trouvé une solution SIMD beaucoup plus simple et ne nécessitant signed*unsigned produits signed*unsigned . Je ne suis plus convaincu que SIMD (du moins avec AVX2 et AV512) ne peut rivaliser avec mulx . Dans certains cas, SIMD peut concurrencer mulx . Le seul cas que je connaisse concerne la multiplication basée sur la FFT de grands nombres .

L’astuce consistait à faire d’abord une multiplication non signée, puis à corriger. J’ai appris comment faire cela grâce à cette réponse de type multiplication 32 bits signée sans utilisation de données 64 bits . La correction est simple car (hi,lo) = x*y fait d’abord la multiplication non signée, puis corrige hi comme ceci:

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

Cela peut être fait avec le _mm_cmpgt_epi64 insortingnsèque à SSE4.2

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

Le code de la multiplication des signes non signés est plus simple car il n’est pas nécessaire de signed*unsigned produits signed*unsigned . De plus, comme il n'est pas signé, il n'a pas besoin de décalage arithmétique à droite qui ne contient qu'une instruction pour AVX512. En fait, la fonction suivante nécessite uniquement 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; } 

Cette utilise

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

AVX512 a le _mm_mullo_epi64 insortingnsèque qui peut calculer lo avec une instruction. Dans ce cas, l'alternative peut être utilisée (commentez les lignes avec le commentaire (*) et décommentez la ligne alternative):

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

Pour modifier le code de pleine largeur AVX2, remplacez _mm par _mm256 , si128 par si256 et __m128i par __m256i pour AVX512, remplacez-les par _mm512 , si512 et __m512i .