Масштабирование значений байтовых пикселей (y = ax + b) с помощью SSE2 (как поплавки)?

a и b разные для каждого пикселя? Это затруднит векторизацию, если нет шаблона или вы можете его сгенерировать.

Есть ли способ эффективно генерировать a и b в векторах как с фиксированной, так и с плавающей запятой? Если нет, вставка 4 значений FP или 8 16-битных целых чисел может быть хуже, чем просто скалярные операции.

Фиксированная точка

Если a и b вообще можно повторно использовать или сгенерировать с фиксированной точкой, это может быть хорошим вариантом использования для математики с фиксированной точкой. (т.е. целые числа, представляющие значение * 2 ^ шкала). SSE/AVX не имеет умножения 8b*8b->16b; самые маленькие элементы - это слова, поэтому вам нужно распаковывать байты в слова, но не полностью до 32 бит. Это означает, что вы можете обрабатывать в два раза больше данных за одну инструкцию.

Есть инструкция _mm_maddubs_epi16, которая может быть полезна, если b и a изменяются достаточно редко, или вы можете легко сгенерировать вектор с чередованием байтов a*2^4 и b*2^1. Очевидно, это действительно удобно для билинейной интерполяции, но все же работа будет выполнена за нас с минимальным перетасовкой, если мы сможем подготовить вектор a и b.

float a, b;
const int logascale = 4, logbscale=1;
const int ascale = 1<<logascale;  // fixed point scale for a: 2^4
const int bscale = 1<<logbscale;  // fixed point scale for b: 2^1

const __m128i brescale = _mm_set1_epi8(1<<(logascale-logbscale));  // re-scale b to match a in the 16bit temporary result

for (i=0 ; i<n; i+=16) {
    //__m128i avec = get_scaled_a(i);  
    //__m128i bvec = get_scaled_b(i);
    //__m128i ab_lo = _mm_unpacklo_epi8(avec, bvec);
    //__m128i ab_hi = _mm_unpackhi_epi8(avec, bvec);

    __m128i abvec = _mm_set1_epi16( ((int8_t)(bscale*b) << 8) | (int8_t)(ascale*a) );  // integer promotion rules might do sign-extension in the wrong place here, so check this if you actually write it this way.

    __m128i block = _mm_load_si128(&buf[i]);  // call this { v[0] .. v[15] }

    __m128i lo = _mm_unpacklo_epi8(block, brescale);  // {v[0], 8, v[1], 8, ...}
    __m128i hi = _mm_unpackhi_epi8(block, brescale);  // {v[8], 8, v[9], 8, ...
    lo = _mm_maddubs_epi16(lo, abvec);  // first arg is unsigned bytes, 2nd arg is signed bytes
    hi = _mm_maddubs_epi16(hi, abvec);
    // lo = { v[0]*(2^4*a) + 8*(2^1*b), ... }

    lo = _mm_srli_epi16(lo, logascale);  // truncate from scaled fixed-point to integer
    hi = _mm_srli_epi16(hi, logascale);

    // and re-pack.  Logical, not arithmetic right shift means sign bits can't be set
    block = _mm_packuswb(lo, hi);  
    _mm_store_si128(&buf[i], block);
}
// then a scalar cleanup loop

2^4 — произвольный выбор. Он оставляет 3 бита без знака для целой части a и 4 бита дроби. Таким образом, он эффективно округляет a до ближайшего 16-го и переполняется, если его величина превышает 8 и 15/16-е. 2^6 даст больше дробных битов и позволит a от -2 до +1 и 63/64-х.

Поскольку b складывается, а не умножается, его полезный диапазон намного больше, а дробная часть гораздо менее полезна. Чтобы представить его в 8 битах, округление до ближайшей половины по-прежнему сохраняет немного дробной информации, но позволяет ему быть [-64 : 63,5] без переполнения.

Для большей точности хорошим выбором будет фиксированная точка 16b. Вы можете масштабировать a и b на 2 ^ 7 или что-то вроде того, чтобы иметь дробную точность 7b и по-прежнему разрешать целочисленной части быть [-256 .. 255]. Для этого случая нет инструкции по умножению и сложению, поэтому вам придется делать это отдельно. Хорошие варианты выполнения умножения включают:

• _mm_mulhi_epu16: без знака 16b*16b->high16 (биты [31:16]). Полезно, если a не может быть отрицательным
• _mm_mulhi_epi16: со знаком 16b*16b->high16 (биты [31:16]).
• _mm_mulhrs_epi16: подписанные 16b*16b->бит [30:15] из 32b временные, с округлением. С хорошим выбором коэффициента масштабирования для a это должно быть лучше. Насколько я понимаю, SSSE3 ввела эту инструкцию именно для такого использования.
• _mm_mullo_epi16: подписано 16b*16b->low16 (биты [15:0]). Это позволяет только 8 значащих битов для a до того, как результат low16 переполнится, поэтому я думаю, что все, что вы получаете по сравнению с 8-битным решением _mm_maddubs_epi16, - это большая точность для b.

Чтобы использовать их, вы должны получить масштабированные 16-битные векторы значений a и b, а затем:

• распакуйте свои байты с нулем (или pmovzx byte->word), чтобы слова со знаком все еще находились в диапазоне [0..255]
• сдвиньте слова влево на 7.
• умножьте на ваш вектор a из 16b слов, взяв верхнюю половину каждого результата 16*16->32. (например, мул
• сдвиньте вправо здесь, если вам нужны разные масштабы для a и b, чтобы получить более дробную точность для a
• добавьте к этому b.
• сдвиг вправо, чтобы выполнить окончательное усечение от фиксированной точки до [0..255].

При правильном выборе масштаба с фиксированной запятой он должен обрабатывать более широкий диапазон a и b, а также более дробную точность, чем 8-битная фиксированная точка.

Если вы не сдвинете свои байты влево после распаковки их в слова, a должно быть полнодиапазонным, чтобы получить 8 битов, установленных в старшем 16-м разряде результата. Это означало бы очень ограниченный диапазон a, который вы могли бы поддерживать, не усекая свое временное значение менее чем до 8 бит во время умножения. Даже _mm_mulhrs_epi16 не оставляет много места, так как начинается с 30-го бита.

расширить байты до чисел с плавающей запятой

Если вы не можете эффективно генерировать значения a и b с фиксированной точкой для каждого пикселя, может быть лучше преобразовать ваши пиксели в числа с плавающей запятой. Это требует больше распаковки/переупаковки, поэтому задержка и пропускная способность хуже. Стоит изучить создание a и b с фиксированной точкой.

Чтобы сжатое число с плавающей запятой заработало, вам все равно нужно эффективно построить вектор из a значений для 4 соседних пикселей.

Это хороший вариант использования для pmovzx (SSE4.1), потому что он может напрямую переходить от 8b элементов к 32b. Другими вариантами являются SSE2 punpck[l/h]bw/punpck[l/h]wd с несколькими шагами или SSSE3 pshufb для имитации pmovzx. (Вы можете выполнить одну загрузку 16 байт и перетасовать ее 4 различными способами, чтобы распаковать в четыре вектора по 32 байта целых чисел.)

char *buf;

// const __m128i zero = _mm_setzero_si128();
for (i=0 ; i<n; i+=16) {
    __m128 a = get_a(i);
    __m128 b = get_b(i);

    // IDK why there isn't an intrinsic for using `pmovzx` as a load, because it takes a m32 or m64 operand, not m128.  (unlike punpck*)
    __m128i unsigned_dwords = _mm_cvtepu8_epi32((__m128i)(buf+i));  // load 4B at once.
    __m128 floats = _mm_cvtepi32_ps(unsigned_dwords);

    floats = _mm_fmadd_ps(floats, a, b);  // with FMA available, this might as well be 256b vectors, even with the inconvenience of the different lane-crossing semantics of pmovzx vs. punpck
    // or without FMA, do this with _mm_mul_ps and _mm_add_ps
    unsigned_dwords = _mm_cvtps_epi32(floats);

    // repeat 3 more times for buf+4, buf+8, and buf+12, then:

    __m128i packed01 = _mm_packss_epi32(dwords0, dwords1);  // SSE2
    __m128i packed23 = _mm_packss_epi32(dwords2, dwords3);
    // packuswb wants SIGNED input, so do signed saturation on the first step

    // saturate into [0..255] range
    __m12i8 packedbytes=_mm_packus_epi16(packed01, packed23);  // SSE2
    _mm_store_si128(buf+i, packedbytes);  // or storeu if buf isn't aligned.
}

// cleanup code to handle the odd up-to-15 leftover bytes, if n%16 != 0

Предыдущая версия этого ответа исходила из векторов float-> uint8 с packusdw/packuswb и содержала целый раздел, посвященный обходным путям без SSE4.1. Ни один из этих битов маскировки знака после неподписанного пакета не требуется, если вы просто остаетесь в подписанном целочисленном домене до последнего пакета< /а>. Я предполагаю, что это причина, по которой SSE2 включает только подписанный пакет от двойного слова к слову, но и подписанный, и неподписанный пакет от слова к байту. packuswd полезен только в том случае, если вашей конечной целью является uint16_t, а не дальнейшая упаковка.

Я предполагаю, что вы ищете составной __m128 _mm_cvtpi8_ps(__m64 a ) встроенный.

Вот минимальный пример:

#include <xmmintrin.h>
#include <stdio.h>

int main() {
  unsigned char  a[4] __attribute__((aligned(32)))= {1,2,3,4};
  float b[4] __attribute__((aligned(32)));
  _mm_store_ps(b, _mm_cvtpi8_ps(*(__m64*)a));
  printf("%f %f, %f, %f\n", b[0], b[1], b[2], b[3]);
  return 0;
}