Как реализовать спад к нулю в математике со знаком с фиксированной точкой, в SSE? - PullRequest
Новая
       57

Как реализовать спад к нулю в математике со знаком с фиксированной точкой, в SSE?

0 голосов
/ 03 октября 2018

Существует много физических событий, подобных распаду (например, трение тела или утечка заряда ), которые обычно моделируются в итераторах, таких как x' = x * 0.99, которые обычно очень легкописать в арифметике с плавающей запятой.

Однако у меня есть требование сделать это 16-битным способом с фиксированной запятой "8,8", в sse.Для эффективной реализации на типичном ALU упомянутая формула может быть переписана как x = x - x/128; или x = x - (x>>7), где >> - «арифметическое», расширяющее знак смещение вправо.

И я застрял здесь, потому что _mm_sra_epi16() производитабсолютно нелогичное поведение, которое легко проверить с помощью следующего примера: 

#include <cstdint>
#include <iostream>
#include <emmintrin.h>

using namespace std;

int main(int argc, char** argv) {
    cout << "required: ";
    for (int i = -1; i < 7; ++i) {
        cout << hex << (0x7fff >> i) << ", ";
    }
    cout << endl;
    cout << "produced: ";
    __m128i a = _mm_set1_epi16(0x7fff);
    __m128i b = _mm_set_epi16(-1, 0, 1, 2, 3, 4, 5, 6);
    auto c = _mm_sra_epi16(a, b);
    for (auto i = 0; i < 8; ++i) {
        cout << hex << c.m128i_i16[i] << ", ";
    }
    cout << endl;
    return 0;
}

Вывод будет следующим: 

required: 0, 7fff, 3fff, 1fff, fff, 7ff, 3ff, 1ff,
produced: 0, 0, 0, 0, 0, 0, 0, 0,

Это относится только ко всем первым изменениям, как на самом деле _mm_sra1_epi16 функция, случайно названная sra и получившая __m128i второй аргумент, но забавное предложение без причины.Так что это не может быть использовано в SSE.

С другой стороны, я слышал, что алгоритм деления чрезвычайно сложен, поэтому _mm_div_epi16 отсутствует в SSE и также не может быть использован.Что делать и как реализовать / векторизовать эту популярную технику "распада"?

1 Ответ

0 голосов
/ 03 октября 2018

x -= x>>7 тривиально реализовать с помощью SSE2, используя постоянный счетчик сдвигов для повышения эффективности.Это компилируется в 2 инструкции, если AVX доступен, в противном случае movdqa требуется для копирования v перед разрушительным сдвигом вправо. 

__m128i downscale(__m128i v){
    __m128i dec = _mm_srai_epi16(v, 7);
    return _mm_sub_epi16(v, dec);
}

GCC даже автоматически векторизирует его ( Godbolt ). 

void foo(short *__restrict a) {
    for (int i=0 ; i<10240 ; i++) {
        a[i] -= a[i]>>7;  // inner loop uses the same psraw / psubw
    }
}

В отличие от float, фиксированная точка имеет постоянную абсолютную точность во всем диапазоне, а не постоянную относительную точность.Таким образом, для небольших положительных чисел v>>7 будет равен нулю, и ваш декремент остановится.(Отрицательные входные значения уменьшаются до -1, потому что арифметическое смещение вправо округляется в сторону -infinity.)

Если небольшие входные значения, где смещение может опуститься до 0, вы можете захотеть выполнить ИЛИ с помощью _mm_set1_epi16(1), чтобы убедиться в уменьшенииненулевойНезначительное влияние на входы большого размера.Тем не менее, это в конечном итоге приведет к уменьшению цепочки от 0 до -1.(И затем вернемся к 0, потому что -1 | 1 == -1 в дополнении 2). 

__m128i downscale_nonzero(__m128i v){
    __m128i dec = _mm_srai_epi16(v, 7);
    dec = _mm_or_si128(dec, _mm_set1_epi16(1));
    return _mm_sub_epi16(v, dec);
}

Если начать с отрицательного значения, последовательность будет -большой, логарифмической до -128, линейной до -4, -3,-2, -1, 0, -1, 0, -1, ...

Ваш код получил все нули, потому что _mm_sra_epi16 использует младшие 64 бита второго источникавектор как 64-битный счетчик сдвига, который применяется ко всем элементам. Прочтите руководство .Таким образом, вы сдвинули все биты из каждого 16-битного элемента.

Это не идиотизм, но для подсчета сдвига на элемент требуется AVX2 (для 32/64-битных элементов) или AVX512BW для _mm_srav_epi16 или 64-битные арифметические сдвиги вправо, что будет иметь смысл для того, как вы пытаетесь его использовать.(Но счетчик сдвига не имеет знака, поэтому -1 также собирается сдвинуть все биты.)

Действительно, эта инструкция должна называться _mm_sra1_epi16()

Да, это имело бы смысл.Но помните, что когда они были названы, AVX2 _mm_srav_* еще не существовало.Кроме того, это конкретное имя не было бы идеальным, потому что 1 и i не являются наиболее визуально различимыми.(i для непосредственного использования, для формы psraw xmm1, imm16 вместо формы psraw xmm1, xmm2/m128 инструкции asm: http://felixcloutier.com/x86/PSRAW:PSRAD:PSRAQ.html).

Другой способ заключается в том, что в инструкции asm MMX / SSE2 есть дваформы: немедленный (с одинаковым количеством всех элементов) и вектор. Вместо того, чтобы заставлять вас передавать счет всем элементам, векторная версия берет скалярное число в нижней части векторного регистра. Я думаю, что предполагаемое использование-case после movd xmm0, eax или чего-то еще.

Если вам необходимо подсчет смещения для каждого элемента без AVX512 , см. различные вопросы и ответы об эмуляции, например, Сдвиг 4 целых чисел вправо на разные значения SIMD .

Некоторые из обходных путей используют умножения на степени 2 для переменного смещения влево, а затем сдвиг вправо, чтобы поместить данные в нужное место. (Но вам нужнокаким-то образом подготовить вектор 1<<n SIMD, так что это работает, если один и тот же набор счетчиков повторно используется для многих векторов, или, особенно, если это константа времени компиляции).

WВ 16-битных элементах вы можете использовать только один _mm_mulhi_epi16 для подсчета правого смещения во время выполнения без потери точности или пределов диапазона.mulhi(x*y) точно так же, как (x*(int)y) >> 16, поэтому вы можете использовать y=1<<14 для сдвига вправо на 16-14 = 2 в этом элементе.

...