Знак максимального абсолютного значения в __m128, SSE4
Мне нужно знать знак значения, которое имеет максимальное абсолютное значение, хранящееся в __m128. Это решение, которое я имею сейчас:
int getMaxSign(__m128 const& vec) {
static const __m128 SIGN_BIT_MASK =
_mm_castsi128_ps(_mm_set1_epi32(0x80000000));
// This creates an int, where sign(a) is 1 if a is negative, 0 o.w.:
// sign(a3)<<3 | sign(a2)<<2 | sign(a1)<<1 | sign(a0)
const int signMask = _mm_movemask_ps(vec);
// Get the absolute value of the vector;
__m128 absValsMMX = _mm_andnot_ps(SIGN_BIT_MASK, vec);
// Figure out the horizontal max
__declspec(align(16)) float absVals[4];
_mm_store_ps(absVals, absValsMMX);
const float maxVal = std::max(std::max(absVals[0], absVals[1]), absVals[2]);
return (maxVal == absVals[0] ? signMask & 0x1 :
(maxVal == absVals[1] ? signMask & 0x2 : signMask & 0x4));
}
В этом случае знак будет равен 1, если значение с максимальным абсолютным значением было отрицательным, и 0 в противном случае, но на самом деле мне все равно, что это за соглашение. Другое дело, что я представляю однородные векторы, используя эти __m128s, поэтому я знаю, что последнее значение всегда будет 0.
Это кажется большой работой для относительно простой задачи. Как я могу сделать это быстрее?
Спасибо!
Решение
Вот одна из возможных реализаций (в C):
int getMaxSign(const __m128 v)
{
__m128 v1, vmax, vmin, vsign;
float sign;
v1 = (__m128)_mm_alignr_epi8((__m128i)v, (__m128i)v, 4); // v1 = v rotated by 1 element
vmax = _mm_max_ps(v, v1); // generate horizontal max/min
vmin = _mm_min_ps(v, v1);
vmax = _mm_max_ps(vmax, (__m128)_mm_alignr_epi8((__m128i)vmax, (__m128i)vmax, 8));
vmin = _mm_min_ps(vmin, (__m128)_mm_alignr_epi8((__m128i)vmin, (__m128i)vmin, 8));
vsign = _mm_add_ps(vmax, vmin); // add max and min to get sign of abs max
sign = _mm_extract_ps(vsign, 0);
return (int)(sign < 0.0f); // return 1 for negative
}
Хотя это выглядит как много кода, это всего лишь 9 инструкций SSE, и нет доступа к памяти, нет ветвей и очень мало скалярного кода.
Обратите внимание, что инструкции SSSE3 и SSE4.1 используются выше.
Вот вторая версия, которая требует только SSSE3:
int getMaxSign(const __m128 v)
{
__m128 v1, vmax, vmin, vsign;
int mask;
v1 = (__m128)_mm_alignr_epi8((__m128i)v, (__m128i)v, 4); // v1 = v rotated by 1 element
vmax = _mm_max_ps(v, v1); // generate horizontal max/min
vmin = _mm_min_ps(v, v1);
vmax = _mm_max_ps(vmax, (__m128)_mm_alignr_epi8((__m128i)vmax, (__m128i)vmax, 8));
vmin = _mm_min_ps(vmin, (__m128)_mm_alignr_epi8((__m128i)vmin, (__m128i)vmin, 8));
vsign = _mm_add_ps(vmax, vmin); // add max and min to get sign of abs max
mask = _mm_movemask_epi8((__m128i)vsign);
return (mask & 8) != 0; // return 1 for negative
}
Это генерирует 12 инструкций:
pshufd $57, %xmm0, %xmm1
movdqa %xmm0, %xmm2
minps %xmm1, %xmm2
pshufd $78, %xmm2, %xmm3
minps %xmm3, %xmm2
maxps %xmm1, %xmm0
pshufd $78, %xmm0, %xmm1
maxps %xmm1, %xmm0
addps %xmm2, %xmm0
pmovmskb %xmm0, %eax
shrl $3, %eax
andl $1, %eax
Обратите внимание, как хитро меняется компилятор palignr в pshufd а также реализует финальный скалярный тест, используя только shrl и andl,
Примечание для Visual Studio C / C ++: приведение между __m128 а также __m128i вам нужно будет использовать _mm_castps_si128 а также _mm_castsi128_psнапример,
mask = _mm_movemask_epi8((__m128i)vsign);
необходимо изменить на:
mask = _mm_movemask_epi8(_mm_castps_si128(vsign));
Другие решения
Если ваши числа дискретны, правильно расставлены и взяты из ограниченного подмножества, есть и другие возможности.
Если вам гарантировано, что a, b и c, например, являются целыми числами, то вы можете умножить вектор на себя, чтобы получить нечетную степень, а затем поставить точки <1, 1, 1>. Например, если мы умножим это на 4 раза, это даст вам < а ^ 5, б ^ 5, с ^ 5>. Если | a | является наибольшим и | a | = 2, тогда мы знаем, что b и c будут равны 1 или 0, поэтому значение a ^ 3 будет доминировать, а скалярное произведение будет иметь свой знак. Например, если X = < a = -2, b = 1, c = 0>, тогда X ^ 5 = <-32, 1, 0>. Когда вы укажете это <1, 1, 1> вы получите -31, знак которого отражает знак наибольшего абсолютного значения. По мере увеличения абсолютного значения наибольшего числа расхождение между ним и другими терминами будет стремиться сходиться — например, если мы имеем <-8, 7, 7>, то приведенный выше алгоритм дает X ^ 5 =<-32768, 16807, 16807>, вы ставите это <1, 1, 1> и получим 846, поэтому алгоритм завершится неудачно с показателем 5. Если мы увеличим показатель до 7, получим <-2097152, 823543, 823543>, усеяно <1, 1, 1> дает нам -450066, что является правильным ответом. В конечном счете ошибки округления также сломают этот метод. Но я надеюсь, что это может дать некоторое представление о других альтернативах, если вы знаете пределы своего набора данных.
В качестве сноски, помните, что X ^ 5 = (X * X) * (X * X) * X, поэтому вы делаете одно умножение, чтобы получить X ^ 2, умножаете это на себя, чтобы получить X ^ 4, а затем умножаете на X — Всего три умножения. Вам нужен нечетный показатель степени, чтобы сохранить знак.
m = min(a,b,c);
M = max(a,b,c);
// return abs(m)>abs(M) ? sign(m): sign(M); // was
return sign(m+M);
Как правильно заметил Paul_R, знак получается просто из суммы значений min и max. Который когда-либо имеет большее (противоположно подписанное) абсолютное значение, выигрывает.
Но идея может быть использована более: сумма мин / макс такая же, как сумма всех элементов минус средний, который можно найти при макс. 3 сравнениях.
return sign(a+b+c - middle(a,b,c)); // or
return sign(a*aw + b*bw + c*cw); // where aw,bw,cw = [0,1]
aw, bw, cw могут быть получены из числа выигранных сравнений (которые, я думаю, должны быть тщательно спланированы для случая, когда есть 2 или 3 равных значения).
И далее:
x = abs(b)>abs(a)?b:a;
return sign(x+c);
Возможно, даже дальше:
s = sign(a + b); // store the sign of larger of a or b
a = abs(a); b=abs(b);
a = max(a,b) | s; // somehow copy the sign.
return sign(a+c);