SIMD подписано с умножением без знака для 64-битных * 64-битных на 128-битных
Я создал функцию, которая выполняет 64-битные * 64-битные и 128-битные с использованием SIMD. В настоящее время я реализовал это с использованием SSE2 (точнее, SSE4.1). Это означает, что он одновременно работает с двумя продуктами от 64b*64b до 128b. Эта же идея может быть расширена до AVX2 или AVX512, дающих четыре или восемь продуктов от 64b * 64 до 128b одновременно. Я основал свой алгоритм на http://www.hackersdelight.org/hdcodetxt/muldws.c.txt
Этот алгоритм выполняет одно умножение без знака, одно умножение со знаком и два умножения со знаком * без знака. Подписанные * подписанные и неподписанные * неподписанные операции легко сделать с помощью _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;
}
Становится хуже. Здесь нет _mm_srai_epi64 Instinsic/ инструкция до AVX512, поэтому я должен был сделать свой собственный.
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 выше не является полным. Я думаю, что использовал Библиотеку векторных классов Агнера Фога. Если вы посмотрите в файле 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);
}
}
2 ответа
Правильный способ определения пределов пропускной способности целочисленного умножения с использованием различных инструкций заключается в том, сколько "битов продукта" можно вычислить за цикл.
mulx производит один результат 64x64 -> 128 каждый цикл; это 64x64 = 4096 "бит продукта за цикл"
Если вы соберете множитель на SIMD из инструкций, которые выполняют умножения 32x32 -> 64 бита, вы должны иметь возможность получать четыре результата в каждом цикле, чтобы соответствовать mulx (4x32x32 = 4096). Если бы не было арифметики, кроме умножений, вы бы просто взломали AVX2. К сожалению, как вы заметили, существует много арифметических операций, кроме умножений, так что это абсолютно не пусковое устройство для оборудования текущего поколения.
Я нашел SIMD решение, которое намного проще и не нуждается signed*unsigned товары. Я больше не уверен, что SIMD (по крайней мере, с AVX2 и AV512) не может конкурировать с mulx , В некоторых случаях SIMD может конкурировать с mulx, Единственный известный мне случай - это умножение больших чисел на основе БПФ.
Хитрость заключалась в том, чтобы сначала сделать беззнаковое умножение, а затем исправить. Я узнал, как это сделать, из этого ответа 32-битный со знаком-умножение-без-использования-64-битного типа данных. Коррекция проста для (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,