Преобразовать подписанный короткий, чтобы плавать в C ++ SIMD

Question

У меня есть массив с коротким знаком, который я хочу разделить на 2048, и в результате получить массив с плавающей точкой.

Я нашел SSE: преобразовать короткое целое в число с плавающей точкой , что позволяетконвертировать беззнаковые шорты в плавающие, но я хочу также обрабатывать подписанные шорты.

Код ниже работает, но только для положительных шорт.

// We want to divide some signed short by 2048 and get a float.
const auto floatScale = _mm256_set1_ps(2048);

short* shortsInput = /* values from somewhere */;
float* floatsOutput = /* initialized */;

__m128i* m128iInput = (__m128i*)&shortsInput[0];

// Converts the short vectors to 2 float vectors. This works, but only for positive shorts.
__m128i m128iLow = _mm_unpacklo_epi16(m128iInput[0], _mm_setzero_si128());
__m128i m128iHigh = _mm_unpackhi_epi16(m128iInput[0], _mm_setzero_si128());
__m128 m128Low = _mm_cvtepi32_ps(m128iLow);
__m128 m128High = _mm_cvtepi32_ps(m128iHigh);

// Puts the 2 __m128 vectors into 1 __m256.
__m256 singleComplete = _mm256_castps128_ps256(m128Low);
singleComplete = _mm256_insertf128_ps(singleComplete, m128High, 1);

// Finally do the math
__m256 scaledVect = _mm256_div_ps(singleComplete, floatScale);

// and puts the result where needed.
_mm256_storeu_ps(floatsOutput[0], scaledVect);

Как я могуконвертировать мои подписанные шорты в поплавки?Или, может быть, есть более эффективный способ решения этой проблемы?

---

РЕДАКТИРОВАТЬ: я пробовал разные ответы по сравнению с алгоритмом без SIMD, делая это 10M раз по сравнению с массивом 2048, на AMD Ryzen 7 2700 в~ 3.2GHz.Я использую Visual 15.7.3 в основном с конфигурацией по умолчанию:

/permissive- /Yu"stdafx.h" /GS /GL /W3 /Gy /Zc:wchar_t /Zi /Gm- /O2 /sdl 
/Fd"x64\Release\vc141.pdb" /Zc:inline /fp:precise /D "NDEBUG" /D "_CONSOLE"
/D "_UNICODE" /D "UNICODE" /errorReport:prompt /WX- /Zc:forScope
/arch:AVX2 /Gd /Oi /MD /openmp /FC /Fa"x64\Release\" /EHsc /nologo
/Fo"x64\Release\" /Fp"x64\Release\test.pch" /diagnostics:classic 

Обратите внимание, что я очень новичок в SIMD и давно не использовал C ++.Вот что я получаю (я перезапускаю каждый тест отдельно, а не один за другим и получаю лучшие результаты, как это):

• Нет SIMD: 7300 мс
• ответ Вима: 2300 мс
• Ответ chtz на SSE2: 1650 мс
• Ответ chtz на AVX2: 2100 мс

Так что я получаю хорошее ускорение с помощью SIMD и ответа chtz на SSE2, хотя и более многословен и сложен дляпонимаю, быстрее.(По крайней мере, при компиляции с включенным AVX, поэтому он избегает дополнительных инструкций по копированию регистров с использованием 3-операндных VEX-кодированных инструкций. В процессорах Intel версии AVX2 должны быть значительно быстрее, чем 128-битная версия.)

Вот мой тестовый код:

const int size = 2048;
const int loopSize = (int)1e7;

float* noSimd(short* shortsInput) {
    float* floatsOutput = new float[size];

    auto startTime = std::chrono::high_resolution_clock::now();

    for (int i = 0; i < loopSize; i++) {
        for (int j = 0; j < size; j++) {
            floatsOutput[j] = shortsInput[j] / 2048.0f;
        }
    }

    auto stopTime = std::chrono::high_resolution_clock::now();
    long long totalTime = (stopTime - startTime).count();

    printf("%lld noSimd\n", totalTime);

    return floatsOutput;
}

float* wimMethod(short* shortsInput) {
    const auto floatScale = _mm256_set1_ps(1.0f / 2048.0f);
    float* floatsOutput = new float[size];

    auto startTime = std::chrono::high_resolution_clock::now();

    for (int i = 0; i < loopSize; i++) {
        for (int j = 0; j < size; j += 8) {
            __m128i short_vec = _mm_loadu_si128((__m128i*)&shortsInput[j]);
            __m256i int_vec = _mm256_cvtepi16_epi32(short_vec);
            __m256  singleComplete = _mm256_cvtepi32_ps(int_vec);

            // Finally do the math
            __m256 scaledVect = _mm256_mul_ps(singleComplete, floatScale);

            // and puts the result where needed.
            _mm256_storeu_ps(&floatsOutput[j], scaledVect);
        }
    }

    auto stopTime = std::chrono::high_resolution_clock::now();
    long long totalTime = (stopTime - startTime).count();

    printf("%lld wimMethod\n", totalTime);

    return floatsOutput;
}

float* chtzMethodSSE2(short* shortsInput) {
    float* floatsOutput = new float[size];

    auto startTime = std::chrono::high_resolution_clock::now();

    for (int i = 0; i < loopSize; i++) {
        for (int j = 0; j < size; j += 8) {
            // get input:
            __m128i val = _mm_loadu_si128((__m128i*)&shortsInput[j]);
            // add 0x8000 to wrap to unsigned short domain:
            val = _mm_add_epi16(val, const0x8000);
            // interleave with upper part of float(1<<23)/2048.f:
            __m128i lo = _mm_unpacklo_epi16(val, const0x4580);
            __m128i hi = _mm_unpackhi_epi16(val, const0x4580);
            // interpret as float and subtract float((1<<23) + (0x8000))/2048.f
            __m128 lo_f = _mm_sub_ps(_mm_castsi128_ps(lo), constFloat);
            __m128 hi_f = _mm_sub_ps(_mm_castsi128_ps(hi), constFloat);
            // store:
            _mm_storeu_ps(&floatsOutput[j], lo_f);
            _mm_storeu_ps(&floatsOutput[j] + 4, hi_f);
        }
    }

    auto stopTime = std::chrono::high_resolution_clock::now();
    long long totalTime = (stopTime - startTime).count();

    printf("%lld chtzMethod\n", totalTime);

    return floatsOutput;
}

float* chtzMethodAVX2(short* shortsInput) {
    const auto floatScale = _mm256_set1_ps(1.0f / 2048.0f);
    float* floatsOutput = new float[size];

    auto startTime = std::chrono::high_resolution_clock::now();

    for (int i = 0; i < loopSize; i++) {
        for (int j = 0; j < size; j += 8) {

            // get input:
            __m128i val = _mm_loadu_si128((__m128i*)&shortsInput[j]);
            // interleave with 0x0000
            __m256i val_unpacked = _mm256_cvtepu16_epi32(val);

            // 0x4580'8000
            const __m256 magic = _mm256_set1_ps(float((1 << 23) + (1 << 15)) / 2048.f);
            const __m256i magic_i = _mm256_castps_si256(magic);

            /// convert by xor-ing and subtracting magic value:
            // VPXOR avoids port5 bottlenecks on Intel CPUs before SKL
            __m256 val_f = _mm256_castsi256_ps(_mm256_xor_si256(val_unpacked, magic_i));
            __m256 converted = _mm256_sub_ps(val_f, magic);
            // store:
            _mm256_storeu_ps(&floatsOutput[j], converted);
        }
    }

    auto stopTime = std::chrono::high_resolution_clock::now();
    long long totalTime = (stopTime - startTime).count();

    printf("%lld chtzMethod2\n", totalTime);

    return floatsOutput;
}

Answer 1

Вы можете заменить стандартный способ преобразования epi16-> epi32-> float и умножения на 1.f/2048.f, вручную составив float.

Это работает, потому что делитель имеет степень 2, поэтомуручное создание числа с плавающей запятой означает просто другой показатель степени.

Благодаря @PeterCordes, вот оптимизированная версия этой идеи для AVX2, использующая XOR для установки старших байтов 32-разрядного числа с плавающей запятой в то же время, что и отражениезнаковый бит целочисленного значения.FP SUB превращает эти младшие биты мантиссы в правильное значение FP:

// get input:
__m128i val = _mm_loadu_si128((__m128i*)input);
// interleave with 0x0000
__m256i val_unpacked = _mm256_cvtepu16_epi32(val);

// 0x4580'8000
const __m256 magic = _mm256_set1_ps(float((1<<23) + (1<<15))/2048.f);
const __m256i magic_i = _mm256_castps_si256(magic);

/// convert by xor-ing and subtracting magic value:
// VPXOR avoids port5 bottlenecks on Intel CPUs before SKL
__m256 val_f = _mm256_castsi256_ps(_mm256_xor_si256(val_unpacked, magic_i));
__m256 converted = _mm256_sub_ps(val_f, magic);
// store:
_mm256_storeu_ps(output, converted);

Посмотрите это в проводнике компилятора Godbolt с помощью gcc и clang ;на Skylake i7-6700k цикл элементов 2048, который горячий в кеше, занимает ~ 360 тактов, с той же скоростью (с точностью до погрешности измерения), что и версия @ wim, которая выполняет стандартное расширение / преобразование / умножение (с аналогичным количествомцикл раскатывания).Протестировано @PeterCordes с Linux perf.Но на Райзене это может быть значительно быстрее, потому что мы избегаем _mm256_cvtepi32_ps (Райзен имеет пропускную способность 1 на 2 такта для vcvtdq2ps ymm: http://agner.org/optimize/.)

Xor 0x8000 с нижней половиной эквивалентнок добавлению / вычитанию 0x8000, поскольку переполнение / перенос игнорируется. И по совпадению, это позволяет использовать одну и ту же магическую константу для XOR-ввода и вычитания.

Как ни странно, gcc и clang предпочитают заменять вычитание надобавление -magic, которое не будет повторно использовать константу ... Они предпочитают использовать add, потому что оно коммутативное, но в этом случае нет никакой выгоды, потому что они не используют его с операндом памяти.

---

Вот версия SSE2, которая выполняет переворачивание со знаком / без знака отдельно от установки старших 2 байтов 32-разрядного битового шаблона FP.

Мы используем один _mm_add_epi16, два_mm_unpackXX_epi16 и два _mm_sub_ps для 8 значений (_mm_castsi128_ps не используются, и _mm_set будет кэшироваться в регистрах):

// get input:
__m128i val = _mm_loadu_si128((__m128i*)input);
// add 0x8000 to wrap to unsigned short domain:
// val = _mm_add_epi16(val, _mm_set1_epi16(0x8000));
val = _mm_xor_si128(val, _mm_set1_epi16(0x8000));  // PXOR runs on more ports, avoids competing with FP add/sub or unpack on Sandybridge/Haswell.

// interleave with upper part of float(1<<23)/2048.f:
__m128i lo = _mm_unpacklo_epi16(val, _mm_set1_epi16(0x4580));
__m128i hi = _mm_unpackhi_epi16(val, _mm_set1_epi16(0x4580));
// interpret as float and subtract float((1<<23) + (0x8000))/2048.f
__m128 lo_f = _mm_sub_ps(_mm_castsi128_ps(lo), _mm_set_ps1(float((1<<23) + (1<<15))/2048.f));
__m128 hi_f = _mm_sub_ps(_mm_castsi128_ps(hi), _mm_set_ps1(float((1<<23) + (1<<15))/2048.f));
// store:
_mm_storeu_ps(output, lo_f);
_mm_storeu_ps(output+4, hi_f);

Демонстрация использования: https://ideone.com/b8BfJd

Если ваш вклад будетесли бы не было подписано short , то _mm_add_epi16 не понадобилось бы (и, конечно, 1<<15 в _mm_sub_ps необходимо было бы удалить).Тогда у вас будет ответ Марата на SSE: преобразовать короткое целое число в число с плавающей точкой .

Это может легко быть перенесено в AVX2 с вдвое большим числом преобразований на итерацию, нонеобходимо позаботиться о порядке элементов вывода (спасибо @wim за указание на это).

---

Кроме того, для чистого решения SSE можно просто использовать _mm_cvtpi16_ps, но это Intelбиблиотечная функция.Нет единой инструкции, которая делает это.

// cast input pointer:
__m64* input64 = (__m64*)input;
// convert and scale:
__m128 lo_f = _mm_mul_ps(_mm_cvtpi16_ps(input64[0]), _mm_set_ps1(1.f/2048.f));
__m128 hi_f = _mm_mul_ps(_mm_cvtpi16_ps(input64[1]), _mm_set_ps1(1.f/2048.f));

Я не тестировал ни одно решение (ни проверял теоретические пропускные способности или задержки)

Answer 2

В AVX2 нет необходимости отдельно преобразовывать верхнюю и нижнюю части:

const auto floatScale = _mm256_set1_ps(1.0f/2048.0f);

short* shortsInput = /* values from somewhere */;
float* floatsOutput = /* initialized */;

__m128i short_vec = _mm_loadu_si128((__m128i*)shortsInput);
__m256i int_vec =  _mm256_cvtepi16_epi32 (short_vec);
__m256  singleComplete = _mm256_cvtepi32_ps (int_vec);

// Finally do the math
__m256 scaledVect = _mm256_mul_ps(singleComplete, floatScale);

// and puts the result where needed.
_mm256_storeu_ps(floatsOutput, scaledVect);

Это прекрасно компилируется в проводнике компилятора Godbolt , и с горячими вводами / выводами в L1dкэш и выровненные массивы ввода / вывода, преобразует массив из 2048 элементов за ~ 360 тактов на Skylake i7-6700k (протестировано в цикле повтора).Это ~ 0,18 цикла на элемент или ~ 5,7 преобразования на такт.Или ~ 1,4 цикла на вектор, включая магазин.Это в основном узкое место по пропускной способности внешнего интерфейса (3,75 мопов в слитых доменах за такт), даже при развертывании цикла clang, потому что преобразование составляет 5 мопов.single uop даже с простым режимом адресации на Haswell / Skylake, поэтому в этом случае хорошо, что последние gcc / clang преобразуют приращения указателя в индексированную адресацию с помощью одного счетчика цикла.С большинством векторных инструкций источника памяти (например, vpmovsxwd xmm, [mem]) это стоило бы дополнительного уопа: Режимы микросинтеза и адресации .

С одной загрузкой и одним хранилищем, это нормально, чтохранилища не могут работать в хранилище данных Haswell / Skylake port7, которое обрабатывает только неиндексированные режимы адресации.

Развертывание цикла необходимо для максимальной пропускной способности ЦП Intel (если нет узких мест в памяти), посколькузагрузить + конвертировать + магазин уже 4 моп.То же, что и в ответе @ chtz.

Идеально использовать векторное значение для дальнейших вычислений сразу, если вам нужно только прочитать значения с плавающей запятой пару раз.Это всего лишь 3 инструкции (но есть некая задержка для скрытого исполнения exec-порядка).Повторное преобразование при необходимости может быть лучше, чем иметь больший объем кэша для хранения вдвое большего float[] результата в памяти;это зависит от вашего варианта использования и аппаратного обеспечения.