В отличие от 32x32 => 64, здесь нет инструкции умножения SSE с расширением 16x16 -> 32.
Вместо этого есть _mm_mulhi_epi16 и _mm_mulhi_epu16, которые дают только подписанные илибеззнаковая верхняя половина полного результата.
(и _mm_mullo_epi16, что делает упакованное 16x16 => 16-битовое умножение на половину усечения с половиной, то же самое для знака или без знака).
Вы можете использовать _mm_unpacklo/hi_epi16 для чередования половинок низких / высоких частот в пару векторов с 32-битными элементами, но это будет довольно медленно.Но да, вы можете _mm_srai_epi32(v, 8+4) арифметически сдвинуть вправо на 12, а затем переупаковать, возможно, с помощью _mm_packs_epi32 (насыщение со знаком обратно в 16-бит).Тогда я думаю, проверить на насыщенность?
Ваш случай использования необычен.Есть _mm_mulhrs_epi16, который дает старшие 17 бит, округленные и затем усеченные до 16 бит.(См. Описание).Это полезно для некоторых алгоритмов с фиксированной запятой, где входы масштабируются, чтобы поместить результат в верхнюю половину, и вы хотите округлить, включая младшую половину вместо усечения.
Вы могли бы фактически использовать _mm_mulhrs_epi16 или _mm_mulhi_epi16 как лучшая ставка для сохранения максимальной точности, может быть, сдвинув v0 влево перед тем, как возвести в квадрат до той точки, где верхняя половина даст вам (v0*v0) >> (8+4).
Как вы думаете,проще не допустить переполнения результата, а просто создать маску с _mm_cmpge_epi16(v1, vThreshold), как это делает автор в оригинальной статье?
Черт, да!получение еще одного или двух битов точности обойдется вам, возможно, в 2 раза в производительности, потому что вам придется вычислять другой результат умножения, чтобы проверить на переполнение, или эффективно увеличить до 32-бит (сокращение числа элементов на вектор в два раза)), как описано выше.
С результатом сравнения v0 = ( v1 & ~m ) | (vR & m ); становится смесью SSE4.1: _mm_blendv_epi8.
Есливаш vThreshold имеет 2 неустановленных бита вверху, у вас есть место для сдвига влево, не теряя ни одного из наиболее значимых битов .Поскольку mulhi дает вам (v0*v0) >> 16, то вы можете сделать это:
// losing the high 2 bits of v0
__m128i v0_lshift2 = _mm_slli_epi16(v0, 2); // left by 2 before squaring
__m128i v0_sqr_asr12 = _mm_mulhi_epi16(v0_lshift2, v0_lshift2);
__m128i v1 = _mm_add_epi16(v0, I);
v1 = _mm_add_epi16(v1, v0_sqr_asr12);
// v1 = ((v0<<2)* (int)(v0<<2))) >> 16) + v0 + I
// v1 = ((v0*(int)v0) >> 12) + v0 + I
Сдвиг влево на 2 до возведения в квадрат такой же, как сдвиг влево на 4 после возведения в квадрат (из полного 32-разрядного результата),Он помещает 16 бит, которые мы хотим, в верхние 16 точно.
Но это непригодно, если ваш v0 настолько близок к полному диапазону, что вы могли бы переполниться при сдвиге влево.
В противном случае вы можете потерять 6 младших битов v0 перед умножением
Округление до -Infinity с арифметическим сдвигом вправо теряет 6 бит точности, но переполнение невозможно.
// losing the low 6 bits of v0
__m128i v0_asr6 = _mm_srai_epi16(v0, 6);
__m128i v0_sqr_asr12 = _mm_mullo_epi16(v0_asr6, v0_asr6);
__m128i v1 = _mm_add_epi16(v0, I);
v1 = _mm_add_epi16(v1, v0_sqr_asr12);
// v1 = (v0>>6) * (int)(v0>>6)) + v0 + I
// v1 ~= ((v0*(int)v0) >> 12) + v0 + I
Я думаю, что таким образом вы теряете больше точности, поэтому, вероятно, лучше установить vThreshold достаточно маленьким, чтобы у вас было достаточно накладных расходов, чтобы использовать умножения на половину.Этот способ включает в себя, возможно, худшее округление.
pmulhrsw округлять вместо усечения может быть даже лучше, если мы можем настроить его так же эффективно.Но я не думаю, что мы можем, потому что сдвиг вправо на 1 является нечетным числом.Я думаю, что нам нужно было бы сделать 2 отдельных входа, один v0_lshift2 и один только влево, смещенный на 1.