Написать

Тёмная тема

Написать

Тёмная тема

ssec++

SSE целочисленное деление?

Существует _mm_div_ps для деления значений с плавающей точкой, есть _mm_mullo_epi16 для целочисленного умножения. Но есть ли что-то для целочисленного деления (16-битное значение)? Как я могу провести такое разделение?

 Marc Glisse29 мая 2013 г., 23:56
Быстрый способ проверить: скомпилироватьtypedef short vec __attribute__((vector_size(16))); vec f(vec x){ return x/7; } vec g(vec x,vec y){ return x/y; } с помощью gcc и посмотрите, как он эмулирует деление на константу, но уменьшает полное деление на скаляры. Если вы можете найти лучший код, пожалуйста, помогите, заполнив отчет об ошибке с помощью gcc. (лязг должен работать так же, но у меня нетпроверил)
 rwong17 апр. 2015 г., 09:52
Смотрите также:libdivide
 Mysticial29 мая 2013 г., 22:49
Нет, это не 'не существует. Это'Вероятно, это сочетание недостаточно людей, нуждающихся в этом, а также сложность (и размер пространства) целочисленной единицы деления.
 Marc Glisse29 мая 2013 г., 22:27
Этого не хватало людям, так что ...
 Paul R29 мая 2013 г., 23:12
Вы хотите разделить на константу или переменную?
 fogbit30 мая 2013 г., 07:19
@PaulR: мне нужно разделить на переменную, int16 / int16
 user208879031 мая 2013 г., 10:11
Я добавил редактирование, которое делает это с SSE. Тем не менее, теперь, когда я думаю об этом, может оказаться более эффективным просто сохранить шорты в массиве, выполнить деление без SSE, а затем загрузить их обратно в регистр.

Ответы на вопрос(4)

это для новичков благодаря оригинальному решению: этоАгнер Фогьбиблиотека подпрограмм сделал со мной волшебство в оптимизации

Вот'Ситуация, когда вы делите на одно и то же значение переменной несколько раз (как в большом цикле)

#include 

unsigned int a, b, d;
unsigned int divisor = any_random_value;
div_u32 OptimumDivision(divisor);
a/OptimumDivision;
b/OptimumDivision;

тот's для unsigned int - если вам нужно использовать отрицательное значениеdiv_i32 вместо этого, который был быстрее в моих тестах, даже если руководство говорит обратное

Я получаю около 3х или больше производительности

Мат говорит, что действительно можно идти быстрее

Агнер Фогьс (http://www.agner.org/optimize/#vectorclass) метод отлично работает, если деление выполняется одним делителем. Кроме того, этот метод имеет еще больше преимуществ, если делитель известен во время компиляции или если он нечасто меняются во время выполнения.

Однако при выполнении SIMD деления на__m128i Элементы, такие, что нет никакой информации как для делителя, так и для дивиденда во время компиляции, у нас нет выбора, кроме как преобразовать в плавающее и выполнить вычисление. С другой стороны, используя_mm_div_ps не даст нам удивительного улучшения скорости, так как эта инструкция имеет переменную задержку от 11 до 14 циклов на большинстве микроархитектур и иногда может доходить до 38 циклов, если мы рассмотрим Knights Landing. Кроме того, эта инструкция не полностью конвейерна и имеет обратную пропускную способность 3-6 циклов в зависимости от микроархитектуры.

Однако мы можем избежать_mm_div_ps и использовать_mm_rcp_ss вместо.

к несчастью__m128 _mm_rcp_ss (__m128 a) быстрый только потому, что он обеспечивает приближение. А именно (взято изIntel Intrnisics Guide):

Вычислить приблизительную обратную величину нижнего 32-битного элемента с плавающей запятой одинарной точности в a, сохранить результат в нижнем элементе dst и скопировать 3 верхних упакованных элемента из a в верхние элементы dst. Максимальная относительная погрешность для этого приближения составляет менее 1,5 * 2 ^ -12.

Поэтому, чтобы извлечь выгоду из_mm_rcp_ssнам нужно компенсировать потери из-за приближения. В этом направлении проделана большая работа, доступная вУлучшено деление по инвариантным целым числам, автор Niels MöИлер и Торбьер-н Гранлунд:

Из-за отсутствия эффективных инструкций деления в современных процессорах деление выполняется как умножение с использованием предварительно вычисленного однословного приближения обратной величины делителя, после чего следует пара шагов настройки.

Чтобы вычислить 16-разрядное целочисленное деление со знаком, нам нужен только один шаг настройки, и мы соответствующим образом основываем наше решение.

SSE2
static inline __m128i _mm_div_epi16(const __m128i &a_epi16, const __m128i &b_epi16) {
    //
    // Setup the constants.
    //
    const __m128  two     = _mm_set1_ps(2.00000051757f);
    const __m128i lo_mask = _mm_set1_epi32(0xFFFF);
    //
    // Convert to two 32-bit integers
    //
    const __m128i a_hi_epi32       = _mm_srai_epi32(a_epi16, 16);
    const __m128i a_lo_epi32_shift = _mm_slli_epi32(a_epi16, 16);
    const __m128i a_lo_epi32       = _mm_srai_epi32(a_lo_epi32_shift, 16);
    const __m128i b_hi_epi32       = _mm_srai_epi32(b_epi16, 16);
    const __m128i b_lo_epi32_shift = _mm_slli_epi32(b_epi16, 16);
    const __m128i b_lo_epi32       = _mm_srai_epi32(b_lo_epi32_shift, 16);
    //
    // Convert to 32-bit floats
    //
    const __m128 a_hi = _mm_cvtepi32_ps(a_hi_epi32);
    const __m128 a_lo = _mm_cvtepi32_ps(a_lo_epi32);
    const __m128 b_hi = _mm_cvtepi32_ps(b_hi_epi32);
    const __m128 b_lo = _mm_cvtepi32_ps(b_lo_epi32);
    //
    // Calculate the reciprocal
    //
    const __m128 b_hi_rcp = _mm_rcp_ps(b_hi);
    const __m128 b_lo_rcp = _mm_rcp_ps(b_lo);
    //
    // Calculate the inverse
    //
    #ifdef __FMA__
        const __m128 b_hi_inv_1 = _mm_fnmadd_ps(b_hi_rcp, b_hi, two);
        const __m128 b_lo_inv_1 = _mm_fnmadd_ps(b_lo_rcp, b_lo, two);
    #else
        const __m128 b_mul_hi   = _mm_mul_ps(b_hi_rcp, b_hi);
        const __m128 b_mul_lo   = _mm_mul_ps(b_lo_rcp, b_lo);
        const __m128 b_hi_inv_1 = _mm_sub_ps(two, b_mul_hi);
        const __m128 b_lo_inv_1 = _mm_sub_ps(two, b_mul_lo);
    #endif
    //
    // Compensate for the loss
    //
    const __m128 b_hi_rcp_1 = _mm_mul_ps(b_hi_rcp, b_hi_inv_1);
    const __m128 b_lo_rcp_1 = _mm_mul_ps(b_lo_rcp, b_lo_inv_1);
    //
    // Perform the division by multiplication
    //
    const __m128 hi = _mm_mul_ps(a_hi, b_hi_rcp_1);
    const __m128 lo = _mm_mul_ps(a_lo, b_lo_rcp_1);
    //
    // Convert back to integers
    //
    const __m128i hi_epi32 = _mm_cvttps_epi32(hi);
    const __m128i lo_epi32 = _mm_cvttps_epi32(lo);
    //
    // Zero-out the unnecessary parts
    //
    const __m128i hi_epi32_shift = _mm_slli_epi32(hi_epi32, 16);
    #ifdef __SSE4_1__
        //
        // Blend the bits, and return
        //
        return _mm_blend_epi16(lo_epi32, hi_epi32_shift, 0xAA);
    #else
        //
        // Blend the bits, and return
        //
        const __m128i lo_epi32_mask = _mm_and_si128(lo_epi32, const_mm_div_epi16_lo_mask);
        return _mm_or_si128(hi_epi32_shift, lo_epi32_mask);
    #endif
}

Это решение может работать с использованиемSSE2 только и будет использоватьFMA если доступно. Однако может быть так, что использование простого деления происходит так же быстро (или даже быстрее), как и использование аппроксимации.

В присутствииAVX это решение может быть улучшено, так как верхняя и нижняя части могут обрабатываться одновременно с использованием одногоAVX регистр.

Проверка

Поскольку мы имеем дело только с 16-битным кодом, мы можем легко проверить правильность решения за несколько секунд, используя тестирование методом грубой силы:

 void print_epi16(__m128i a)
{
    int i; int16_t tmp[8];
    _mm_storeu_si128( (__m128i*) tmp, a);

    for (i = 0; i < 8; i += 1) {
        printf("%8d ", (int) tmp[i]);
    }
    printf("\n");
}

bool run_mm_div_epi16(const int16_t *a, const int16_t *b)
{
    const size_t n = 8;
    int16_t result_expected[n];
    int16_t result_obtained[n];
    //
    // Derive the expected result
    //
    for (size_t i = 0; i < n; i += 1) {
        result_expected[i] = a[i] / b[i];
    }
    //
    // Now perform the computation
    //
    const __m128i va = _mm_loadu_si128((__m128i *) a);
    const __m128i vb = _mm_loadu_si128((__m128i *) b);
    const __m128i vr = _mm_div_epi16(va, vb);
    _mm_storeu_si128((__m128i *) result_obtained, vr);
    //
    // Check for array equality
    //
    bool eq = std::equal(result_obtained, result_obtained + n, result_expected);
    if (!eq) {
        cout < "Testing of _mm_div_epi16 failed" < endl < endl;
        cout < "a: ";
        print_epi16(va);
        cout < "b: ";
        print_epi16(vb);
        cout < endl;
        cout < "results_obtained: ";
        print_epi16(vr);
        cout < "results_expected: ";
        print_epi16(_mm_loadu_si128((__m128i *) result_expected));
        cout < endl;
    }
    return eq;
}

void test_mm_div_epi16()
{
    const int n = 8;
    bool correct = true;
    //
    // Brute-force testing
    //
    int16_t a[n];
    int16_t b[n];

    for (int32_t i = INT16_MIN; correct && i
 Peter Cordes22 июл. 2018 г., 08:54
Забыл сказать раньше: если деление является частью большего алгоритма,divps часто лучше чем итерация rcpps + Newton для общей пропускной способности на современных процессорах, (Кроме КНЛ, где тыВы собираетесь использовать AVX512ER.) Если вы неУзкое место на не полностью конвейерном разделителе, 128 битdivps это просто один моп, против нескольких. rcpps + Newton имеет схожую задержку с divps. OoO exec в любом случае скрывает задержку, еслине переносится петлей.
 Alen Stojanov23 июл. 2018 г., 03:03
@PeterCordes - спасибо за предложения. Я обновил свой ответ, добавив профили производительности, сравнивая простое деление с шагом аппроксимации.
 Alen Stojanov22 июл. 2018 г., 03:26
Согласовано. Я пытался сделать чистыйSSE2 версия иFMA проверка была перенесена из предыдущего теста. Очевидно, что если идти вверх по отношению к ISA, можно работать лучше, чем исходный код. Я должен вероятно отредактировать код и представитьAVX2 версия с тестом, чтобы получить полный ответ.
 Peter Cordes22 июл. 2018 г., 00:04
Если ты'Повторно проверяя поддержку FMA, вы, вероятно, также должны проверить SSE4.1 и использоватьpmovzx расширять знак, по крайней мере, для нижней половины 128-битного вектора. Кстати, сдвиг SSE считает насыщенным, так что я думаю, что вы можетеsrai_epi16(a_epi16, 16) а потомpunpckl/h с этим, так что знак расширяется до 32-разрядных целых чисел и настраивается для перепаковки сpacks/usdw, (Если вы хотите усечение вместо насыщенности, сначала замаскируйте. Shift / и / OR или shift /pblendw более эффективен, чем / и / pack, но распаковка с перемешиванием может быть выигрышем.)

Для 8-битного деления это может быть реализовано путем создания таблицы магических чисел.

Увидеть "Хакер»восторгстр. 238

подписали:

__m128i _mm_div_epi8(__m128i a, __m128i b)
{
    __m128i abs_b = _mm_abs_epi8(b);

    static const uint16_t magic_number_table[129] =
    {
        0x0000, 0x0000, 0x8080, 0x5580, 0x4040, 0x3380, 0x2ac0, 0x24c0, 0x2020, 0x1c80, 0x19c0, 0x1760, 0x1560, 0x13c0, 0x1260, 0x1120,
        0x1010, 0x0f20, 0x0e40, 0x0d80, 0x0ce0, 0x0c40, 0x0bb0, 0x0b30, 0x0ab0, 0x0a40, 0x09e0, 0x0980, 0x0930, 0x08e0, 0x0890, 0x0850,
        0x0808, 0x07d0, 0x0790, 0x0758, 0x0720, 0x06f0, 0x06c0, 0x0698, 0x0670, 0x0640, 0x0620, 0x05f8, 0x05d8, 0x05b8, 0x0598, 0x0578,
        0x0558, 0x0540, 0x0520, 0x0508, 0x04f0, 0x04d8, 0x04c0, 0x04b0, 0x0498, 0x0480, 0x0470, 0x0458, 0x0448, 0x0438, 0x0428, 0x0418,
        0x0404, 0x03f8, 0x03e8, 0x03d8, 0x03c8, 0x03b8, 0x03ac, 0x03a0, 0x0390, 0x0388, 0x0378, 0x0370, 0x0360, 0x0358, 0x034c, 0x0340,
        0x0338, 0x032c, 0x0320, 0x0318, 0x0310, 0x0308, 0x02fc, 0x02f4, 0x02ec, 0x02e4, 0x02dc, 0x02d4, 0x02cc, 0x02c4, 0x02bc, 0x02b4,
        0x02ac, 0x02a8, 0x02a0, 0x0298, 0x0290, 0x028c, 0x0284, 0x0280, 0x0278, 0x0274, 0x026c, 0x0268, 0x0260, 0x025c, 0x0258, 0x0250,
        0x024c, 0x0248, 0x0240, 0x023c, 0x0238, 0x0234, 0x022c, 0x0228, 0x0224, 0x0220, 0x021c, 0x0218, 0x0214, 0x0210, 0x020c, 0x0208,
        0x0202
    };

    Uint8 load_den[16];
    _mm_storeu_si128((__m128i*)load_den, abs_b);

    uint16_t mul[16];

    for (size_t i = 0; i < 16; i++)
    {
        uint16_t cur_den = load_den[i];
        mul[i] = magic_number_table[cur_den];
    }
    // for denominator 1, magic number is 0x10080 that 16bit-overflow occurs.
    __m128i one = _mm_set1_epi8(1);
    __m128i is_one = _mm_cmpeq_epi8(abs_b, one);

    // -128/-128 is a special case where magic number does not work.
    __m128i v80 = _mm_set1_epi8(0x80);
    __m128i is_80_a = _mm_cmpeq_epi8(a, v80);
    __m128i is_80_b = _mm_cmpeq_epi8(b, v80);
    __m128i is_80 = _mm_and_si128(is_80_a, is_80_b);

    // __m128i zero = _mm_setzero_si128();
    // __m128i less_a = _mm_cmpgt_epi8(zero, a);
    // __m128i less_b = _mm_cmpgt_epi8(zero, b);
    // __m128i  sign = _mm_xor_si128(less_a, less_b);
    __m128i abs_a = _mm_abs_epi8(a);
#if 0
    __m128i p = _mm_unpacklo_epi8(abs_a, zero);
    __m128i q = _mm_unpackhi_epi8(abs_a, zero);
    __m256i c = _mm256_castsi128_si256(p);
    c = _mm256_insertf128_si256(c, q, 1);
#else
    // Thanks to Peter Cordes
    __m256i c = _mm256_cvtepu8_epi16(abs_a);
#endif
    __m256i magic = _mm256_loadu_si256((const __m256i*)mul);
    __m256i high = _mm256_mulhi_epu16(magic, c);
    __m128i v0h = _mm256_extractf128_si256(high, 0);
    __m128i v0l = _mm256_extractf128_si256(high, 1);
    __m128i res = _mm_packus_epi16(v0h, v0l);
    __m128i div = _mm_blendv_epi8(res, abs_a, is_one);
    // __m128i neg = _mm_sub_epi8(zero, div);
    // __m128i select = _mm_blendv_epi8(div, neg, sign);
    __m128i select = _mm_sign_epi8(div, _mm_or_si128(_mm_xor_si128(a, b), one));
    return _mm_blendv_epi8(select, one, is_80);
}

без знака:

__m128i _mm_div_epu8(__m128i n, __m128i den)
{
    static const uint16_t magic_number_table[256] =
    {
        0x0001, 0x0000, 0x8000, 0x5580, 0x4000, 0x3340, 0x2ac0, 0x04a0, 0x2000, 0x1c80, 0x19a0, 0x0750, 0x1560, 0x13c0, 0x0250, 0x1120,
        0x1000, 0x0f10, 0x0e40, 0x0d80, 0x0cd0, 0x0438, 0x03a8, 0x0328, 0x0ab0, 0x0a40, 0x09e0, 0x0980, 0x0128, 0x00d8, 0x0890, 0x0048,
        0x0800, 0x07c8, 0x0788, 0x0758, 0x0720, 0x06f0, 0x06c0, 0x0294, 0x0668, 0x0640, 0x021c, 0x05f8, 0x05d8, 0x01b4, 0x0194, 0x0578,
        0x0558, 0x013c, 0x0520, 0x0508, 0x04f0, 0x04d8, 0x04c0, 0x04a8, 0x0094, 0x0480, 0x006c, 0x0458, 0x0448, 0x0034, 0x0024, 0x0014,
        0x0400, 0x03f4, 0x03e4, 0x03d4, 0x03c8, 0x03b8, 0x03ac, 0x039c, 0x0390, 0x0384, 0x0378, 0x036c, 0x0360, 0x0354, 0x014a, 0x0340,
        0x0334, 0x032c, 0x0320, 0x0318, 0x010e, 0x0304, 0x02fc, 0x02f4, 0x02ec, 0x02e4, 0x02dc, 0x02d4, 0x02cc, 0x02c4, 0x02bc, 0x02b4,
        0x02ac, 0x02a4, 0x02a0, 0x0298, 0x0290, 0x028c, 0x0284, 0x007e, 0x0278, 0x0072, 0x026c, 0x0066, 0x0260, 0x025c, 0x0254, 0x0250,
        0x004a, 0x0244, 0x0240, 0x023c, 0x0036, 0x0032, 0x022c, 0x0228, 0x0224, 0x001e, 0x001a, 0x0016, 0x0012, 0x000e, 0x000a, 0x0006,
        0x0200, 0x00fd, 0x01fc, 0x01f8, 0x01f4, 0x01f0, 0x01ec, 0x01e8, 0x01e4, 0x01e0, 0x01dc, 0x01d8, 0x01d6, 0x01d4, 0x01d0, 0x01cc,
        0x01c8, 0x01c4, 0x01c2, 0x01c0, 0x01bc, 0x01b8, 0x01b6, 0x01b4, 0x01b0, 0x01ae, 0x01ac, 0x01a8, 0x01a6, 0x01a4, 0x01a0, 0x019e,
        0x019c, 0x0198, 0x0196, 0x0194, 0x0190, 0x018e, 0x018c, 0x018a, 0x0188, 0x0184, 0x0182, 0x0180, 0x017e, 0x017c, 0x017a, 0x0178,
        0x0176, 0x0174, 0x0172, 0x0170, 0x016e, 0x016c, 0x016a, 0x0168, 0x0166, 0x0164, 0x0162, 0x0160, 0x015e, 0x015c, 0x015a, 0x0158,
        0x0156, 0x0154, 0x0152, 0x0051, 0x0150, 0x014e, 0x014c, 0x014a, 0x0148, 0x0047, 0x0146, 0x0144, 0x0142, 0x0140, 0x003f, 0x013e,
        0x013c, 0x013a, 0x0039, 0x0138, 0x0136, 0x0134, 0x0033, 0x0132, 0x0130, 0x002f, 0x012e, 0x012c, 0x012a, 0x0029, 0x0128, 0x0126,
        0x0025, 0x0124, 0x0122, 0x0021, 0x0120, 0x001f, 0x011e, 0x011c, 0x001b, 0x011a, 0x0019, 0x0118, 0x0116, 0x0015, 0x0114, 0x0013,
        0x0112, 0x0110, 0x000f, 0x010e, 0x000d, 0x010c, 0x000b, 0x010a, 0x0009, 0x0108, 0x0007, 0x0106, 0x0005, 0x0104, 0x0003, 0x0102
    };

    static const uint16_t shift_table[256] =
    {
        0x0001, 0x0100, 0x0100, 0x0080, 0x0100, 0x0040, 0x0040, 0x0020, 0x0100, 0x0080, 0x0020, 0x0010, 0x0020, 0x0040, 0x0010, 0x0020,
        0x0100, 0x0010, 0x0040, 0x0080, 0x0010, 0x0008, 0x0008, 0x0008, 0x0010, 0x0040, 0x0020, 0x0080, 0x0008, 0x0008, 0x0010, 0x0008,
        0x0100, 0x0008, 0x0008, 0x0008, 0x0020, 0x0010, 0x0040, 0x0004, 0x0008, 0x0040, 0x0004, 0x0008, 0x0008, 0x0004, 0x0004, 0x0008,
        0x0008, 0x0004, 0x0020, 0x0008, 0x0010, 0x0008, 0x0040, 0x0008, 0x0004, 0x0080, 0x0004, 0x0008, 0x0008, 0x0004, 0x0004, 0x0004,
        0x0100, 0x0004, 0x0004, 0x0004, 0x0008, 0x0008, 0x0004, 0x0004, 0x0010, 0x0004, 0x0008, 0x0004, 0x0020, 0x0004, 0x0002, 0x0040,
        0x0004, 0x0004, 0x0020, 0x0008, 0x0002, 0x0004, 0x0004, 0x0004, 0x0004, 0x0004, 0x0004, 0x0004, 0x0004, 0x0004, 0x0004, 0x0004,
        0x0004, 0x0004, 0x0020, 0x0008, 0x0010, 0x0004, 0x0004, 0x0002, 0x0008, 0x0002, 0x0004, 0x0002, 0x0020, 0x0004, 0x0004, 0x0010,
        0x0002, 0x0004, 0x0040, 0x0004, 0x0002, 0x0002, 0x0004, 0x0008, 0x0004, 0x0002, 0x0002, 0x0002, 0x0002, 0x0002, 0x0002, 0x0002,
        0x0100, 0x0001, 0x0004, 0x0008, 0x0004, 0x0010, 0x0004, 0x0008, 0x0004, 0x0020, 0x0004, 0x0008, 0x0002, 0x0004, 0x0010, 0x0004,
        0x0008, 0x0004, 0x0002, 0x0040, 0x0004, 0x0008, 0x0002, 0x0004, 0x0010, 0x0002, 0x0004, 0x0008, 0x0002, 0x0004, 0x0020, 0x0002,
        0x0004, 0x0008, 0x0002, 0x0004, 0x0010, 0x0002, 0x0004, 0x0002, 0x0008, 0x0004, 0x0002, 0x0080, 0x0002, 0x0004, 0x0002, 0x0008,
        0x0002, 0x0004, 0x0002, 0x0010, 0x0002, 0x0004, 0x0002, 0x0008, 0x0002, 0x0004, 0x0002, 0x0020, 0x0002, 0x0004, 0x0002, 0x0008,
        0x0002, 0x0004, 0x0002, 0x0001, 0x0010, 0x0002, 0x0004, 0x0002, 0x0008, 0x0001, 0x0002, 0x0004, 0x0002, 0x0040, 0x0001, 0x0002,
        0x0004, 0x0002, 0x0001, 0x0008, 0x0002, 0x0004, 0x0001, 0x0002, 0x0010, 0x0001, 0x0002, 0x0004, 0x0002, 0x0001, 0x0008, 0x0002,
        0x0001, 0x0004, 0x0002, 0x0001, 0x0020, 0x0001, 0x0002, 0x0004, 0x0001, 0x0002, 0x0001, 0x0008, 0x0002, 0x0001, 0x0004, 0x0001,
        0x0002, 0x0010, 0x0001, 0x0002, 0x0001, 0x0004, 0x0001, 0x0002, 0x0001, 0x0008, 0x0001, 0x0002, 0x0001, 0x0004, 0x0001, 0x0002
    };

    static const uint16_t mask_table[256] =
    {
        0x0000, 0xffff, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0xffff, 0x0000, 0x0000, 0x0000, 0xffff, 0x0000, 0x0000, 0xffff, 0x0000,
        0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0xffff, 0xffff, 0xffff, 0x0000, 0x0000, 0x0000, 0x0000, 0xffff, 0xffff, 0x0000, 0xffff,
        0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0xffff, 0x0000, 0x0000, 0xffff, 0x0000, 0x0000, 0xffff, 0xffff, 0x0000,
        0x0000, 0xffff, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0xffff, 0x0000, 0xffff, 0x0000, 0x0000, 0xffff, 0xffff, 0xffff,
        0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0xffff, 0x0000,
        0x0000, 0x0000, 0x0000, 0x0000, 0xffff, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,
        0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0xffff, 0x0000, 0xffff, 0x0000, 0xffff, 0x0000, 0x0000, 0x0000, 0x0000,
        0xffff, 0x0000, 0x0000, 0x0000, 0xffff, 0xffff, 0x0000, 0x0000, 0x0000, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff,
        0x0000, 0xffff, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,
        0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,
        0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,
        0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,
        0x0000, 0x0000, 0x0000, 0xffff, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0xffff, 0x0000, 0x0000, 0x0000, 0x0000, 0xffff, 0x0000,
        0x0000, 0x0000, 0xffff, 0x0000, 0x0000, 0x0000, 0xffff, 0x0000, 0x0000, 0xffff, 0x0000, 0x0000, 0x0000, 0xffff, 0x0000, 0x0000,
        0xffff, 0x0000, 0x0000, 0xffff, 0x0000, 0xffff, 0x0000, 0x0000, 0xffff, 0x0000, 0xffff, 0x0000, 0x0000, 0xffff, 0x0000, 0xffff,
        0x0000, 0x0000, 0xffff, 0x0000, 0xffff, 0x0000, 0xffff, 0x0000, 0xffff, 0x0000, 0xffff, 0x0000, 0xffff, 0x0000, 0xffff, 0x0000
    };

    uint8_t load_den[16];
    _mm_storeu_si128((__m128i*)load_den, den);

    uint16_t mul[16];
    uint16_t mask[16];
    uint16_t shift[16];

    for (size_t i = 0; i < 16; i++)
    {
        const uint16_t cur_den = load_den[i];
        mul[i] = magic_number_table[cur_den];
        mask[i] = mask_table[cur_den];
        shift[i] = shift_table[cur_den];
    }
#if 0
    __m128i a = _mm_unpacklo_epi8(n, _mm_setzero_si128());
    __m128i b = _mm_unpackhi_epi8(n, _mm_setzero_si128());
    __m256i c = _mm256_castsi128_si256(a);
    c = _mm256_insertf128_si256(c, b, 1);
#else
    // Thanks to Peter Cordes
    __m256i c = _mm256_cvtepu8_epi16(n);
#endif
    __m256i magic = _mm256_loadu_si256((const __m256i*)mul);
    __m256i high = _mm256_mulhi_epu16(magic, c);
    __m256i low = _mm256_mullo_epi16(magic, c);
    __m256i low_down = _mm256_srli_epi16(low, 8);
    __m256i high_up = _mm256_slli_epi16(high, 8);
    __m256i low_high = _mm256_or_si256(low_down, high_up);
    __m256i target_up = _mm256_mullo_epi16(c, _mm256_loadu_si256((const __m256i*)shift));
    __m256i cal1 = _mm256_sub_epi16(target_up, low_high);
    __m256i cal2 = _mm256_srli_epi16(cal1, 1);
    __m256i cal3 = _mm256_add_epi16(cal2, low_high);
    __m256i cal4 = _mm256_srli_epi16(cal3, 7);
    __m256i res = _mm256_blendv_epi8(high, cal4, _mm256_loadu_si256((const __m256i*)mask));

    __m128i v0h = _mm256_extractf128_si256(res, 0);
    __m128i v0l = _mm256_extractf128_si256(res, 1);

    return _mm_packus_epi16(v0h, v0l);
}
 Peter Cordes01 дек. 2018 г., 19:30
Если ты'Вы собираетесь использовать AVX2, вы можете и должны использоватьVPMOVZXBW создаватьc с одной инструкцией вместо распаковки lo / hi отдельно и перемешивания вместе._mm256_cvtepu8_epi16
 sugwan kim01 дек. 2018 г., 21:27
@Peter Спасибо за ваш совет.
 Peter Cordes01 дек. 2018 г., 20:24
Вы также можете оптимизировать вычисление знака: так как вам нужно только старший бит, чтобы исправитьblendv, вы можете__m128i sign = _mm_xor_si128(a,b);, Или даже лучше, условно отрицать использованиеvpsignb, Кроме этогоa^b будет ноль дляa==bгде результат должен быть +1, а не 0. Так(a^b) | 1 будет работать, потому что у вас уже естьset1(1) вектор, который сделает его ненулевым, не влияя на бит знака. Так__m128i select = _mm_sign_epi8(div, a^b | one);, (vblendv это 2 мопа и работает только на случайном порте на Haswell)

Пожалуйста, смотрите Агнер ТуманВ векторном классе он реализовал быстрый алгоритм для целочисленного деления с SSE / AVX для 8-битных, 16-битных и 32-битных слов (но не 64-битных)http://www.agner.org/optimize/#vectorclass

Посмотрите в файле vectori128.h код и описание алгоритма как его хорошо написанное руководство VectorClass.pdf

Вот фрагмент, описывающий алгоритм из его руководства ".

Целочисленное деление В наборе команд x86 и его расширениях нет инструкций, полезных для целочисленного векторного деления, и такие инструкции были бы довольно медленными, если бы они существовали. Поэтому библиотека векторных классов использует алгоритм быстрого целочисленного деления. Основной принцип этого алгоритма может быть выражен в этой формуле: a / b ≈ a * (2n / b) >> n Этот расчет проходит следующие шаги: 1. найти подходящее значение для n 2. вычислить 2n / b 3. вычислить необходимые поправки для ошибок округления 4. выполнить умножение и сдвиг вправо и применить исправления для ошибок округления

Эта формула полезна, если несколько чисел разделены одним и тем же делителем b. Шаги 1, 2 и 3 необходимо выполнить только один раз, а шаг 4 повторяется для каждого значения дивиденда a. Математические детали описаны в файле vectori128.h. (См. Также Т. Гранлунд и П. Л. Монтгомери: Деление по инвариантным целым числам с использованием умножения, материалы SIGPLAN. "...

Редактировать: ближе к концу файла vectori128.h показывает, как сделать короткое деление со скалярной переменной "Для вычисления параметров, используемых для быстрого деления, требуется больше времени, чем для деления. Следовательно, выгодно использовать один и тот же объект делителя несколько раз. Например, чтобы разделить 80 беззнаковых коротких целых на 10: "

short x = 10;
uint16_t dividends[80], quotients[80];         // numbers to work with
Divisor_us div10(x);                          // make divisor object for dividing by 10
Vec8us temp;                                   // temporary vector
for (int i = 0; i < 80; i += 8) {              // loop for 4 elements per iteration
    temp.load(dividends+i);                    // load 4 elements
    temp /= div10;                             // divide each element by 10
    temp.store(quotients+i);                   // store 4 elements
}

Изменить: целочисленное деление на вектор шорт

#include 
#include "vectorclass.h"

int main() {    
    short numa[] = {10, 20, 30, 40, 50, 60, 70, 80};
    short dena[] = {10, 20, 30, 40, 50, 60, 70, 80};

    Vec8s num = Vec8s().load(numa);
    Vec8s den = Vec8s().load(dena);

    Vec4f num_low = to_float(extend_low(num));
    Vec4f num_high = to_float(extend_high(num));
    Vec4f den_low = to_float(extend_low(den));
    Vec4f den_high = to_float(extend_high(den));

    Vec4f qf_low = num_low/den_low;
    Vec4f qf_high = num_high/den_high;
    Vec4i q_low = truncate_to_int(qf_low);
    Vec4i q_high = truncate_to_int(qf_high);

    Vec8s q = compress(q_low, q_high);
    for(int i=0; i
 user208879031 мая 2013 г., 10:05
@PaulR, поэтому я посмотрел на это более внимательно. Быстрый метод в VectorClass работает только для скаляров. Для целочисленного деления на вектор / вектор, я думаю, выисправить. Единственное решение, если вам нужно сделать векторное / векторное целое число, это либо преобразовать в числа с плавающей точкой, выполнить деление и преобразовать обратно в шорты. Я добавил редактирование, делая это. Я думаю, вы могли бы также сохранить шорты в массиве, сделать скалярное деление, а затем загрузить их обратно.
 user208879031 мая 2013 г., 12:33
@PaulR, может быть, у нас есть другое определение константы. Когда я слышу константу, я думаю о константе времени компиляции. Я думаю, что под константой вы подразумеваете скалярное значение. Его метод хорошо работает для переменного скаляра, а также для постоянного скаляра времени компиляции. Деление на константу времени компиляции происходит еще быстрее. Различия описаны на странице 13 его руководства, если вы заботитесь.
 Paul R31 мая 2013 г., 11:53
Да, это в значительной степени подводит итог - вы можете выполнить деление на (одну) константу, не прибегая к плавающей точке, но все остальное требует преобразования в плавающую точку и обратно.
 user208879031 мая 2013 г., 12:49
@PaulR, я думаю, мы сейчас согласны.
 Paul R30 мая 2013 г., 12:40
Это действительно полезно только для деления на константу - ОП сказал, что он хочет делить на переменную, и в этом случае я думаю, что единственный вариант - распаковать, преобразовать в float и сделать деление в float.
 Paul R31 мая 2013 г., 12:42
Хорошо - яМы никогда не пробовали этот метод во время выполнения (только во время компиляции), но кажется, что стоимость установки относительно высока, если вы не собираетесь обрабатывать значительное количество элементов с одинаковым скалярным дивидендом, поэтому я думаю, что это 'по-прежнему скалярная константа в том смысле, что в идеале это нет меняется внутри внутреннего цикла. Если у вас есть дивиденды, которые могут изменяться на протяжении итерации внутреннего цикла, метод становится совершенно неэффективным.
 user208879030 мая 2013 г., 16:07
@PaulR, я добавил правку в ответ на твой комментарий. В конце файла vectori128.h он приводит пример того, как выполнять деление шорт с помощью переменных. Я не изучил детали, потому что целочисленное деление еще не было узким местом в моем коде (когда мне приходилось его использовать, я делал это, сдвигая биты вправо, что очень быстро).

Ваш ответ на вопрос

Популярные вопросы

0 ответов

Реализация системного вызова Linux с использованием LKM

0 ответов

pytz и astimezone () не могут быть применены к наивному datetime

0 ответов

Получить метаданные изображения - Android

0 ответов

Связывание данных с DropDownList MVC Razor

0 ответов

Загрузите файл Excel и извлеките из него данные и поместите эти данные в базу данных, используя MVC asp.net