Введение

Деление — достаточно затратная операция. Например, на CPU Cannon Lake задержки 32-битного деления находятся в интервале 10-15 тактов, а на Zen4 — 9-14 тактов. Задержки 32-битного умножения на обоих CPU составляют 3-4 такта.

Ни в одном из популярных ISA SIMD (SSE, AVX, AVX-512, ARM Neon, ARM SVE) нет целочисленного деления, оно есть только в RISC-V Vector Extension. Однако во всех этих ISA есть деление с плавающей запятой.

В этой статье мы представим два подхода к реализации SIMD-деления 8-битных беззнаковых чисел:

1. с помощью деления с плавающей запятой,

2. с помощью алгоритма деления столбиком.

Мы попробуем векторизировать показанную ниже процедуру C++. Процедура не может делать никаких допущений о делимых, особенно если все они равны. Таким образом, мы не можем использовать деление на константу.

void scalar_div_u8(const uint8_t* a, const uint8_t* b, uint8_t* out, size_t n) {
    for (size_t i=0; i < n; i++) {
        out[i] = a[i] / b[i];
    }
}

Компиляторы не могут векторизировать эту процедуру. Например, GCC 14.1.0 генерирует такой ассемблерный код (вырезан весь мусор, связанный с выравниванием кода):

2b40:       48 85 c9                test   %rcx,%rcx
2b43:       74 30                   je     2b75 <_Z13scalar_div_u8PKhS0_Phm+0x35>
2b45:       45 31 c0                xor    %r8d,%r8d
2b60:       42 0f b6 04 07          movzbl (%rdi,%r8,1),%eax
2b65:       42 f6 34 06             divb   (%rsi,%r8,1)
2b69:       42 88 04 02             mov    %al,(%rdx,%r8,1)
2b6d:       49 ff c0                inc    %r8
2b70:       4c 39 c1                cmp    %r8,%rcx
2b73:       75 eb                   jne    2b60 <_Z13scalar_div_u8PKhS0_Phm+0x20>
2b75:       c3                      ret

Операции с плавающей запятой

8-битное число можно без потери точности преобразовать в число с плавающей запятой одинарной точности.

Общий процесс деления состоит из следующих этапов:

1. преобразуем тип 8-битных делимого и делителя в 32-битные integer,

2. преобразуем беззнаковые integer в числа с плавающей запятой,

3. выполняем деление с плавающей запятой,

4. преобразуем результат с плавающей запятой в 32-битные integer,

5. выполняем обратное преобразование типа 32-битного integer в 8-битный окончательный результат.

Деление с округлением

Вот настоящая реализация процедуры SSE. Стоит отметить, что перед обратным преобразованием в integer нужно явным образом усечь число с плавающей запятой. По умолчанию это преобразование округляет аргумент, так что мы можем получить ошибочные результаты (смещённые на 1).

1. Загружаем четыре 8-битных делимых.

   uint32_t buf_a;
   memcpy(&buf_a, &a[i], 4);

2. И четыре 8-битных делителя.

   uint32_t buf_b;
   memcpy(&buf_b, &b[i], 4);

3. Переносим их в регистр SSE и преобразуем тип в 32-битный.

   const __m128i a_u8  = _mm_cvtsi32_si128(buf_a);
   const __m128i a_u32 = _mm_cvtepu8_epi32(a_u8);
   
   const __m128i b_u8  = _mm_cvtsi32_si128(buf_b);
   const __m128i b_u32 = _mm_cvtepu8_epi32(b_u8);

4. Преобразуем тип 32-битных integer во float.

   const __m128  a_f32 = _mm_cvtepi32_ps(a_u32);
   const __m128  b_f32 = _mm_cvtepi32_ps(b_u32);

5. Выполняем деление, а затем усечение.

   const __m128  c_f32   = _mm_div_ps(a_f32, b_f32);
   const __m128  c_f32_2 = _mm_round_ps(c_f32, _MM_FROUND_TO_ZERO | _MM_FROUND_NO_EXC);

6. Преобразуем float обратно в integer.

   const __m128i c_i32 = _mm_cvtps_epi32(c_f32_2);

7. Преобразуем 32-битный тип в 8-битный: собираем младшие 8-битные числа в одно 32-битное слово и сохраняем его в выходной массив.

   const __m128i c_u8  = _mm_shuffle_epi8(c_i32, _mm_setr_epi8(
       0, 4, 8, 12,
       -1, -1, -1, -1,
       -1, -1, -1, -1,
       -1, -1, -1, -1
   ));
   
   const uint32_t buf = _mm_cvtsi128_si32(c_u8);
   memcpy(&out[i], &buf, 4);

Деление без округления

Команда округления ROUNDPS имеет довольно большую задержку, по крайней мере, на CPU Intel. На IceLake — 8 тактов, а на Zen4 — всего 3 такта.

Мы можем избежать округления чисел плавающей запятой, умножив делимое на 256 (выполнив сдвиг влево на 8 битов), а затем сдвинув вправо на 8 окончательный результат. Сдвиг вправо можно выполнить без лишних затрат, потому что мы и так используем перетасовку для получения отдельных битов, так что вопрос только в константе. Сдвиг влево на 8 подходит только для SSE-кода — мы можем использовать перетасовку байтов (byte shuffle), чтобы выполнить сдвиг и расширение целых чисел. В случае AVX2-кода перетасовка байтов выполняется на 128-битных трактах, поэтому нам потребуется больше работы, чтобы подготовить ввод для этой операции.

Эта SSE-процедура практически идентична процедуре из предыдущего раздела:

1. Загружаем четыре делимых.

   uint32_t buf_a;
   memcpy(&buf_a, &a[i], 4);
   const __m128i a_u8  = _mm_cvtsi32_si128(buf_a);

2. Преобразуем делимое << 8 в 32-битные числа.

   const __m128i a_u32 = _mm_shuffle_epi8(a_u8, _mm_setr_epi8(
       -1, 0, -1, -1,
       -1, 1, -1, -1,
       -1, 2, -1, -1,
       -1, 3, -1, -1
   ));

3. Загружаем четыре делителя и преобразуем их в 32-битные числа.

   uint32_t buf_b;
   memcpy(&buf_b, &b[i], 4);
   const __m128i b_u8  = _mm_cvtsi32_si128(buf_b);
   const __m128i b_u32 = _mm_cvtepu8_epi32(b_u8);

4. Преобразуем тип всех 32-битных integer во float.

   const __m128  a_f32 = _mm_cvtepi32_ps(a_u32);
   const __m128  b_f32 = _mm_cvtepi32_ps(b_u32);

5. Выполняем деление.

   const __m128  c_f32   = _mm_div_ps(a_f32, b_f32);

6. Преобразуем результат деления в 32-битные integer.

   const __m128i c_i32 = _mm_cvtps_epi32(c_f32);

7. Преобразуем результат деления >> 8 в 8-битные числа: получаем бит 1 каждого 32-битного слова.

   const __m128i c_u8  = _mm_shuffle_epi8(c_i32, _mm_setr_epi8(
       0 + 1, 4 + 1, 8 + 1, 12 + 1,
       -1, -1, -1, -1,
       -1, -1, -1, -1,
       -1, -1, -1, -1
   ));
   
   const uint32_t buf = _mm_cvtsi128_si32(c_u8);
   memcpy(&out[i], &buf, 4);

Использование приближенной обратной величины

В SSE есть команда RCPPS, вычисляющая приближенную инверсию её аргумента: 1/x. Она позволит нам использовать выражение делимое ⋅ approx(1/делитель).

В спецификации говорится, что относительная погрешность не превышает 1.5 ⋅ 2 ^{− 12}. Но в нашей задаче абсолютная погрешность оказывается слишком большой, чтобы использовать результат команды напрямую. В таблице ниже показаны результаты для первых x значений, для которых погрешность значима.

Однако путём проб и ошибок мы выяснили, что после умножения делимого на значение 1.00025 результат RCPPS можно использовать. Если точнее, подойдёт любой коэффициент от 1.00024 до 1.00199.

Деление столбиком

Алгоритм деления столбиком известен нам со школы.

1. Мы начинаем с остатка и результата деления, равного нулю.

2. Затем на каждом этапе мы расширяем остаток, добавляя в конец следующий разряд делимого, начиная с самых старших разрядов.

3. Если остаток оказывается больше делителя, то мы вычисляем q := остаток/делитель. Это число, разряд, находится в интервале [0, 9] для десятичного деления и равно 0 или 1 для двоичного деления.

4. Выполняем декремент остатка на q ⋅ divisor.

5. Разряд q мы добавляем в начало результата деления.

6. Если в делимом ещё остаются разряды, переходим к пункту 2.

Удобное свойство этого алгоритма заключается в вычислении и результата деления, и остатка. Его недостатки:

• Каждый этап зависит от предыдущего, из-за чего сложно использовать внеочередное выполнение в CPU.

• Количество итераций равно количество значимых разрядов в делимом минус количество разрядов в делителе плюс 1. В наихудшем случае оно равно количеству разрядов в делимом.

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

uint8_t long_div_u8(uint8_t dividend, uint8_t divisor) {
    uint8_t reminder = 0;
    uint8_t quotient = 0;

    for (int i=7; i >= 0; i--) {
        // Освобождаем место для i-того бита делимого в 0-й позиции
        reminder <<= 1;

        // вставляем этот бит
        reminder |= (dividend >> i) & 1;

        if (reminder >= divisor) {
            // задаём i-тый бит в результате деления
            quotient |= 1 << i;

            // изменяем остаток
            reminder -= divisor;
        }
    }

    return quotient;
}

Деление столбиком можно применять и к целым числам со знаком. Нужно вычислить абсолютные значения делимого и делителя, а затем выполнить алгоритм. Если операнды имели разные знаки, то результат нужно сделать отрицательным.

Таккак 8-битные целые числа со знаком находятся в интервале от − 128 до 127, интервал их абсолютных значений равен 0…128. То есть операнды деления в столбик всё равно будут иметь восемь битов.

Векторизация

Алгоритм состоит из следующих операций:

1. Извлекаем i-тый бит из делителя и сдвигаем его в остаток:

   reminder = (reminder << 1) | ((dividend >> i) & 1);

2. Сравниваем остаток и делитель;

3. Условным образом задаём i-тый бит в результате деления и изменяем остаток:

   quotient |= 1 << i;
   reminder -= divisor;

SSE и AVX2

Этап 1: обновление остатка

В SSE и AVX2 можно очень легко скопировать самый старший бит в самый младший бит 0. Мы сравниваем с нулём число, интерпретируемое как единица со знаком. В результате получаем 0x00 или 0xff ( − 1).

Можно переписать основной цикл так, чтобы делимое сдвигалось влево на 1, благодаря чему мы сможем скопировать этим способом все его биты.

Сдвиг влево на 1 — это простое сложение, то есть очень быстрая операция.

Этап 2: сравнение

В SSE и AVX2 нет беззнакового сравнения, только со знаком. Можно сравнить два беззнаковых числа сравнением со знаком: нужно обратить их самые старшие биты. Это выполняется операции XOR с 0x80.

Стоит отметить, что нам придётся выполнять xor один раз для делителя и восемь раз для остатка.

Этап 3: условные операции

Получив после сравнения маску, мы можем легко вычислить операции с масками.

На самом деле, нам нужен безусловный сдвиг результата деления на 1 с последующим условным заданием i-того бита.

Реализация (SSE)

В представленном ниже коде на C показаны все подробности реализации.

uint8_t long_div_u8(uint8_t dividend, uint8_t divisor) {
    uint8_t reminder = 0;
    uint8_t quotient = 0;
    uint8_t bit = 0x80;

    uint8_t divisor_xored = divisor ^ 0x80;

    for (int i=7; i >= 0; i--) {

        // msb => 0 или -1
        const uint8_t msb = (int8_t)dividend < 0 ? 0xff : 0x00;

        // вставляем бит
        reminder -= msb;

        if ((int8_t)(reminder ^ 0x80) >= (int8_t)(divisor_xored) {
            // задаём  i-тый бит в результате деления
            quotient |= bit;

            // изменяем остаток
            reminder -= divisor;
        }

        bit >>= 1;
        dividend <<= 1;
        // освобождаем место для i-того бита делимого в 0-й позиции
        reminder <<= 1;
    }

    return quotient;
}

Настоящая реализация на SSE.

void sse_long_div_u8(const uint8_t* A, const uint8_t* B, uint8_t* C, size_t n) {
    const __m128i msb  = _mm_set1_epi8(int8_t(0x80));
    const __m128i zero = _mm_set1_epi8(0x00);

    for (size_t i=0; i < n; i += 16) {
        __m128i dividend = _mm_loadu_si128((__m128i*)(&A[i]));

        const __m128i divisor = _mm_loadu_si128((__m128i*)(&B[i]));
        __m128i divisor_xored = _mm_xor_si128(divisor, msb);

        __m128i bit      = msb;
        __m128i reminder = _mm_set1_epi16(0);
        __m128i quotient = _mm_set1_epi16(0);
        for (int j=0; j < 8; j++) {
            // копируем msb делимого в остаток
            const __m128i t0 = _mm_cmplt_epi8(dividend, zero);
            reminder = _mm_sub_epi8(reminder, t0);

            // беззнаковое сравнение делителя и остатка
            const __m128i reminder_xored = _mm_xor_si128(reminder, msb);
            const __m128i gt = _mm_cmpgt_epi8(divisor_xored, reminder_xored);

            // получаем условное вычитание и бит результата деления
            const __m128i cond_divisor      = _mm_andnot_si128(gt, divisor);
            const __m128i cond_quotient_bit = _mm_andnot_si128(gt, bit);

            // условным образом обновляем остаток и результат деления
            reminder = _mm_sub_epi16(reminder, cond_divisor);
            quotient = _mm_or_si128(quotient, cond_quotient_bit);

            // следующий бит для результата деления
            bit = _mm_srli_epi32(bit, 1);

            // помещаем следующий бит из делимого в MSB
            dividend = _mm_add_epi8(dividend, dividend);

            // освобождаем место для бита из делимого
            reminder = _mm_add_epi32(reminder, reminder);
        }

        _mm_storeu_si128((__m128i*)(&C[i]), quotient);
    }
}

AVX-512

Этап 1: обновляем остаток

В отличие от кода для SSE/AVX2, здесь проще выполнять сдвиг вправо, чтобы поместить i-тый бит в нулевую позицию. Значит, изолирование самого младшего бита и объединение его с результатом деления можно выразить одной тернарной операцией.

reminder <<= 1;
t0         = divisor >> i;
reminder   = reminder | (t0 & 1); // тернарная операция

Этап 2: сравнение

AVX512 поддерживает беззнаковое сравнение байтов, в результате которого мы получаем маску.

Этап 3: условные операции

Простое использование операций с масками.

Реализация

Показана настоящая реализация с AVX512. В отличие от SSE-кода, внутренний цикл развёрнут вручную. Кроме того, здесь нет явного использования встроенной функции с тернарной логикой, однако взглянув на ассемблерный код, можно увидеть, что компилятор выполнил слияние двоичной операции.

void avx512_long_div_u8(const uint8_t* A, const uint8_t* B, uint8_t* C, size_t n) {
    const __m512i one = _mm512_set1_epi8(1);

    for (size_t i=0; i < n; i += 64) {
        const __m512i dividend = _mm512_loadu_si512((const __m512*)(&A[i]));
        const __m512i divisor  = _mm512_loadu_si512((const __m512*)(&B[i]));

        const __m512i dividend_bit7 = _mm512_and_epi32(_mm512_srli_epi32(dividend, 7), one);
        const __m512i dividend_bit6 = _mm512_and_epi32(_mm512_srli_epi32(dividend, 6), one);
        const __m512i dividend_bit5 = _mm512_and_epi32(_mm512_srli_epi32(dividend, 5), one);
        const __m512i dividend_bit4 = _mm512_and_epi32(_mm512_srli_epi32(dividend, 4), one);
        const __m512i dividend_bit3 = _mm512_and_epi32(_mm512_srli_epi32(dividend, 3), one);
        const __m512i dividend_bit2 = _mm512_and_epi32(_mm512_srli_epi32(dividend, 2), one);
        const __m512i dividend_bit1 = _mm512_and_epi32(_mm512_srli_epi32(dividend, 1), one);
        const __m512i dividend_bit0 = _mm512_and_epi32(_mm512_srli_epi32(dividend, 0), one);

        __m512i quotient = _mm512_set1_epi32(0);
        __m512i reminder = dividend_bit7;

        {
            const __mmask64 ge = _mm512_cmpge_epu8_mask(reminder, divisor);
            reminder = _mm512_mask_sub_epi8(reminder, ge, reminder, divisor);
            quotient = _mm512_mask_add_epi8(quotient, ge, quotient, one);
        }

        reminder = _mm512_add_epi32(reminder, reminder);
        reminder = _mm512_or_epi32(reminder, dividend_bit6);

        {
            const __mmask64 ge = _mm512_cmpge_epu8_mask(reminder, divisor);
            reminder = _mm512_mask_sub_epi8(reminder, ge, reminder, divisor);
            quotient = _mm512_add_epi32(quotient, quotient);
            quotient = _mm512_mask_add_epi8(quotient, ge, quotient, one);
        }

        reminder = _mm512_add_epi32(reminder, reminder);
        reminder = _mm512_or_epi32(reminder, dividend_bit5);

        {
            const __mmask64 ge = _mm512_cmpge_epu8_mask(reminder, divisor);
            reminder = _mm512_mask_sub_epi8(reminder, ge, reminder, divisor);
            quotient = _mm512_add_epi32(quotient, quotient);
            quotient = _mm512_mask_add_epi8(quotient, ge, quotient, one);
        }

        reminder = _mm512_add_epi32(reminder, reminder);
        reminder = _mm512_or_epi32(reminder, dividend_bit4);

        {
            const __mmask64 ge = _mm512_cmpge_epu8_mask(reminder, divisor);
            reminder = _mm512_mask_sub_epi8(reminder, ge, reminder, divisor);
            quotient = _mm512_add_epi32(quotient, quotient);
            quotient = _mm512_mask_add_epi8(quotient, ge, quotient, one);
        }

        reminder = _mm512_add_epi32(reminder, reminder);
        reminder = _mm512_or_epi32(reminder, dividend_bit3);

        {
            const __mmask64 ge = _mm512_cmpge_epu8_mask(reminder, divisor);
            reminder = _mm512_mask_sub_epi8(reminder, ge, reminder, divisor);
            quotient = _mm512_add_epi32(quotient, quotient);
            quotient = _mm512_mask_add_epi8(quotient, ge, quotient, one);
        }

        reminder = _mm512_add_epi32(reminder, reminder);
        reminder = _mm512_or_epi32(reminder, dividend_bit2);

        {
            const __mmask64 ge = _mm512_cmpge_epu8_mask(reminder, divisor);
            reminder = _mm512_mask_sub_epi8(reminder, ge, reminder, divisor);
            quotient = _mm512_add_epi32(quotient, quotient);
            quotient = _mm512_mask_add_epi8(quotient, ge, quotient, one);
        }

        reminder = _mm512_add_epi32(reminder, reminder);
        reminder = _mm512_or_epi32(reminder, dividend_bit1);

        {
            const __mmask64 ge = _mm512_cmpge_epu8_mask(reminder, divisor);
            reminder = _mm512_mask_sub_epi8(reminder, ge, reminder, divisor);
            quotient = _mm512_add_epi32(quotient, quotient);
            quotient = _mm512_mask_add_epi8(quotient, ge, quotient, one);
        }

        reminder = _mm512_add_epi32(reminder, reminder);
        reminder = _mm512_or_epi32(reminder, dividend_bit0);


        {
            const __mmask64 ge = _mm512_cmpge_epu8_mask(reminder, divisor);
            reminder = _mm512_mask_sub_epi8(reminder, ge, reminder, divisor);
            quotient = _mm512_add_epi8(quotient, quotient);
            quotient = _mm512_mask_add_epi8(quotient, ge, quotient, one);
        }

        _mm512_storeu_si512((__m512*)&C[i], quotient);
    }
}

Результаты эксперимента

Краткие итоги:

• На CPU Intel быстрее всех оказалась реализация деления в столбик с AVX-512.

• На Ryzen быстрее всех был способ с приближенным обратным значением AVX2.

• Стоит отметить, что автоматическая векторизация GCC создавала более качественный по сравнению с написанными вручную вариантами AVX2, или сравнимый с ними.

Все программы бенчмарков компилировались с опциями -O3 -march=native отдельно на каждой машине.

Ryzen 7

• Компилятор: gcc (Debian 14.1.0-5) 14.1.0

• CPU: AMD Ryzen 7 7730U с Radeon Graphics

Skylake-X

• CPU: Intel(R) Xeon(R) W-2104, 3.20GHz

• Компилятор: gcc (GCC) 11.2.0

IceLake

• Компилятор: gcc (GCC) 13.3.1 20240611 (Red Hat 13.3.1-2)

• CPU: Intel(R) Xeon(R) Gold 6338, 2.00GHz

Исходный код

Пример реализации выложен на GitHub.