CRC -- реализация через умножение без переносов

CRC - (Циклический избыточный код; англ. Cyclic redundancy check) используется для контроля целостности данных при передаче по сети или записи данных на носитель. Разные протоколы используют различные варианты, которые различаются разрядностью и образующими полиномами. Для рассчета CRC используется арифметика Галуа, в частности операция умножения. Умножение в конечном поле Галуа можно разделить на две операции: умножение без переносов и редуцирование. В большинстве случаев программисты не сильно задумываются над математикой и используют готовые таблицы и алгоритмы.
Допустим в нашем распоряжении есть процессорная инструкция умножения без переносов. Мы хотим оптимизировать и тем самым ускорить в несколько раз расчет контрольных сумм CRC, в десятки раз. Целью статьи является разработка общего подхода к синтезу алгоритмов содержащих умножение в конечном поле. Процесс создания алгоритмов для расчета CRC оказался весьма трудоемким. Хотя все эти алгоритмы подобны друг другу, но требуют отладки и расчета коэффициентов.
{Статья большая, я ее буду медленно дописывать и исправлять... прошу прощения за возможные неточности.}
На платформе Intel операции умножения без переносов соответствует инструкция pclmul. На платформе ARMv8-А также есть аналогичные инструкции, на базе которых можно эффективно реализовать операцию умножения без переносов, так называемая инструкция полиномиального умножения (PMULL) разрядностью операнда 64бит. На платформе ARMv7-A+NEON есть поддержка в виде инструкций vmull.p8, векторное полиномиальное умножение 8бит. На платформе ARMv8-М +MVE Helium ожидаем поддержку в виде векторных инструкций полиномиального умножения p8 бит и p16 бит

Все варианты операции CRC отличаются разрядностью, полиномами и порядком следования бит в числе.

В нашем программном обеспечении используются операции CRC:
CRC8/BACnet
CRC16/BAСnet
CRC16/MODBUS
CRC32K/Koopman
Кроме того, мы используем
CRC32С/Castagnoli
CRC32
CRC32B

Для описания алгоритмов я буду использовать инлайн-функции, которые призваны генерировать векторные инструкции. Базовые операции:
REFLECT128 -- Разворачивание последовательности бит в векторе
BSWAP128 -- Переворачивание порядка следования байт в векторе
LOAD128U -- Загрузка данных в векторный регистр без выравнивания
SLL128U -- логический сдвиг влево числа 128 бит.
Умножение без переноса:
CL_MUL32 () -- умножение без переносов двух чисел 32 бита, результат 64 бита
CL_MUL32H() -- старшие 32 бита произведения двух 32 битных чисел
CL_MUL32L() -- младшие 32 бита произведения
CL_MUL128() -- умножение младшей или старшей части вектора 64 бит на старшую или младшую часть второго вектора. Результат - число 128 бит.

Что такое CRC

По сути вот что это такое, с математической точки зрения,

CRC = (M(x) • xt) mod P(x), 
где t = deg(P(x)) - степень полинома, M(x)- сообщение(число)

-- сообщение M(x) представляют в виде длинного двоичного числа и подвергают операции редуцирования в конечном поле.
Идея оптимизации сводится к накоплению буфера с отложенной операцией редуцирования, операциям сдвига в конечном поле и последующему финальному редуцированию на буфере. Буфером является векторный регистр.
Традиционно обработка CRC ведется на потоке данных, по одному байту. В реальных задачах возникает необходимость расчета CRC при передаче данных по сети или записи на носитель, с расчетом по шапке и телу пакета. Расчет может производится по частям: сначала от шапки пакета, затем от тела пакета или сразу от всего пакета. Когда чтение данных осуществляется блоками, расчет может производится так же по блокам.
Мы хотим разработать скоростной алгоритм, обрабатывать до 16 байт на один цикл, ориентируясь на максимальную разрядность векторного регистра. Планируем использовать операции умножения без переносов. Мы хотим рассчитывать параллельно до 16 контрольных сумм пользуясь тем, что преобразования линейные. Результат параллельного расчета фрагментов сообщения можно складывать. Операция сложения представлена битовой операцией исключающего ИЛИ (⊕). Сложение может выполняться до или после редуцирования.
Редуцирование - это уменьшение разрядности числа после операции умножения без переноса.
Метод ускорения CRC с использованием инструкций умножения без переноса (pclmul) был предложен до нас[2,3], Intel. Однако в литературе представлен только метод для CRC-32. Кроме того, в предложенной методике операция редуцирования выполняется без использования оперций умножения. В данной статье я привожу варианты алгоритмов другой разрядности 8,16,24,32, 64, и нескольких образующих полиномов, привожу методику рассчета сдвиговых коэффициентов, и констант для редуцирования полиномов. Предложен свой вариант алгоритмов для отраженных полиномов, отраженного порядка бит.
Список литературы, рекомендую предварительно ознакомиться:
[1] A painless guide to CRC error detection algorithms
[2] Intel/ Fast CRC Computation Using PCLMULQDQ Instruction
и несколько реализаций на GitHub и в исходниках Apple, в частности BZIP2.
[3] Intel® Carry-Less Multiplication Instruction and its Usage for Computing the GCM Mode

Операция схлопывания

Сообщение М можно разложить на два фрагмента со смещением:
M1(x) • x^s ⊕ M0(x) = M1 • (x^s mod P) ⊕ M0
Эту операцию называют "схлопывание" (folding), слов не хватает, слова "редуцирование" и "сложение" уже заняты и наполнены смыслом. Сразу предлагаем воспользоваться правилом: редуцирование от редуцирования на результат не влияет. Поэтому мы можем сначала выполнить редуцирование сомножителей, а потом применять операцию умножения.
Константы (x^s mod P) для разных показателей степени 32 64 96 128 .. можно рассчитать заранее и использовать для ускорения алгоритма.

Допустим, максимальное число, над которым удобно производить операцию "схлопывания" = 128 бит. Это прежде всего определяется разрядностью операции CL_MUL. В нашем случае операция может выполняться над 64 битными числами и результат имеет разрядность 128бит. Тогда определим операцию схлопывания так, преобразование вектора 256 бит [M4:M3:M2:M1] в 128 бит:
[C2:C1] = (M4 • K3) ⊕ (M3 • K4) ⊕ [M2:M1]
-- результат безрассудного умножения имеет разрядность 128 бит
K4 = x¹²⁸ mod P
K3 = x¹⁹² mod P
K4 и K3 - это сдвиговые константы, возведение в степень - это операция сдвига на соответствующее число разрядов. Таким образом мы заменяем операцию сдвига в поле на операцию умножения без переноса. Соответственно, для последующего преобразования в CRC-32 нам понадобится операция редуцирования из 128-в-96 и 96-в-64 бит и финальное редуцирование 64-в-32 бита. Надо научиться рассчитывать константы для сдвигов и схлопывания.

Редукция Barrett'a

Редукция основана на такой идее, которая изначально родилась для целых чисел. Хотим подобрать такое число Ux, чтобы результат умножения Ux*Px - давал единицу с высокой точностью. Например, хотим поделить целое число в конечном поле 32 бита на 10, а операция деления напрочь отсутствует в нашем арсенале. Тогда мы можем найти такое целое число, чтобы результат операции (A*Ux)>>32 в точности давал бы результат A/10 для всех целых чисел. По сути мы ищем такое целое число, чтобы оно работало как multiplicative inverse (обратное число по отношению к образующему полиному).
Редукцию Барретта можно было бы определить так:

            T1 = floor(R(x)/x^32) * [1/P(x)];   R/P
            T2 = floor(T1/x^32) * P(x);         int(R/P)*P;
            CRC= (R+int(R/P)*P) mod x^32;       R-int(R/P)*P

только под делением мы понимаем операцию уполовинивания в поле, под сложением - битовую операцию XOR и под произведением - умножение без переносов, под делением на x^32 понимаем сдивиг на 32 разряда вправо. Короче, математическое определение редуцирования мало что дает на практике. Можно сказать, зная алгоритм редуцирования для целых чисел я угадал, как он будет выглядеть в конечном поле.
Операция  floor(R(x)/x^32) в конечном поле это просто сдвиг на 32 разряда вправо, а для отраженного полинома - сдвиг влево.

Функция barrett_calc рассчитывает магическое число Баррета, которое в дальнейшем позволяет редуцировать полином.

/*! Алгоритм №1. рассчитывает константу Барретт'а */
uint64_t barrett_calc(uint64_t poly, int bits)
{
    int n = bits;
    uint64_t r = poly;
    uint64_t v = 0;
    while (--n){
        if (r & (1ULL<<n)) {
            r ^= poly>>(bits-n);
            v |= 1ULL<<n;
        }
    }
    if (r) v|=1;// для деления с остатком округляем в большую сторону
    return v;
}

{Я совсем не помню откуда я взял этот алгоритм.}
Операция редуцирования методом Барретта 64-в-32

/*! Алгоритм №2. Редукция Барретт'а*/
uint32_t R (uint32_t v1, uint32_t v0)
{
    uint32_t t = CL_MUL32H(v1, U) ^ v1;// U имеет разрядность 33 бита
    return CL_MUL32L(t, P) ^ v0;
}

Константа U рассчитывается по алгоритму Барретта. Константа U имеет разрядность 33 бита. Поскольку старший бит в число 32 бита не укладывается, к результату умножения CL_MUL32H может добавляться значение операнда, - та самая единица в старшем разряде - это следует не потерять.

Такой же алгоритм можно использовать для обратного порядка обработки CRC, обратного порядка следования бит.

/*! Алгоритм №3. Редукция Барретт'а для отраженных чисел */
uint32_t RR (uint32_t v1, uint32_t v0)
{
    uint32_t t = CL_MUL32L(v1, UR);
    return CL_MUL32H(t, PR) ^ t ^ v0;// PR имеет разрядность 33 бита
}

Во втором случае используется отраженный полином Pr, с отраженным порядком следования бит. И отраженная константа Ur.

Мы только что определили, как выглядит операция поиска обратного (деления с остатком?) не-пойми-чего [1/P] с помощью алгоритма уполовинивания. см barrett_calc
Надо проверить результат barrett_calc. Если все верно, то произведение [1/P]•P должно давать что-то похожее на единицу с достаточной точностью. Проверяем...
(0x104D101DF • 0x104C11DB7) = 1.00000000.490D678D (65 бит);
Результат интересный, на единицу похож смутно, отдаленно, но похож. Если сдвинуть на 64 разряда, то единица получается. А вот в младших 32 битах лежит магическое число, для следующей операции пригодится - обратите внимание остаток равен сдвиговой константе.
Можно проверить ту же операцию на отраженном полиноме.
(0x1F7011641 • 0x1DB710641) = 1.63CD6124.00000001 (65 бит)
Чтобы поверить в магию умножения над отразить все биты в числе, все 65 бит... Короче, это тоже можно считать единицей, с заданной точностью 32 бита. Т.е. полиномы можно умножать задом наперед, но результат будет вывернут.

Правило для умножения без переноса для чисел с отраженным порядком бит выглядит так:

reflect(A) • reflect(B) == reflect(A • B)<<1

Результат нужно сдвинуть на один разряд. Таким образом, мы можем предложить алгоритм для отраженных полиномов, с использованием той же операции умножения без переносов. Алгоритм будет выглядеть очень похожим, с той лишь разницей, что при работе с отраженными полиномами всегда нужно сдвигать результат умножения на разряд влево.
Мой личный опыт строится на историческом осознании роли CRC и развитии CRC в нашем программном обеспечении, А также из чтения статей и изучения чужого опыта. В контроллерах мы применяем табличные методы расчета CRC. Причем, всегда стоит вопрос минимизации исходного кода прошивки, поэтому мы используем маленькие таблицы. Привожу пример табличного алгоритма, 4бит на шаг алгоритма.

// Алгоритм №4 -- табличный метод 4 бит.
uint32_t CRC32_update   (uint32_t crc, uint8_t val){
    crc = (crc<<4) ^ CRC32_Lookup4[((crc>>28) ^ (val>>4)) & 0xF ];
    crc = (crc<<4) ^ CRC32_Lookup4[((crc>>28) ^ (val   )) & 0xF ];
    return crc;
}

// Алгоритм №5 -- табличный метод 8 бит.
uint32_t CRC32_update   (uint32_t crc, uint8_t val){
    crc = (crc<<8) ^ CRC32_Lookup8[((crc>>24) ^ (val   )) & 0xFF ];
    return crc;
}

Математическое обоснование этих действий.

 [M3:M2:M1]•x8 ==  ((M3•x8 ⊕ M2)•x8 ⊕ M1)•x8

В примере сообщение нарезано по байтам. Видно, что алгоритм будет подобен себе для любой нарезки сообщения: на 1 бит, на 2, на4, на 8, на 16...
Для синтеза таблиц умножения использую метод...

// Алгоритмом № Расчет таблицы умножения, прямой порядок
/*! \brief генерация таблицы подстановки CRC
 \param poly полином
 \param bits число бит в полиноме
 \param size число элементов в таблице 16 или 256
 */
void crc_gen_table(uint64_t poly, int bits, int size)
{
    uint64_t table[size];// = {0};
    uint64_t p =poly;
    int i,j;
    table[0] = 0;
    table[1] = p;
    for (i=1;(1<<i)<size;i++)
    {
        if (p&(1ULL<<(bits-1))) {
            p &= ~((~0ULL)<<(bits-1));
            p = (p<<1) ^ poly;
        } else
            p = (p<<1);
        table[1<<i] = p;
        for(j=1; j<(1<<i); j++) {
            table[(1<<i)+j] = p ^ table[j];
        }
    }
// Распечатать таблицу
    printf("POLY=0x%0*"PRIX64"\n", bits/4, poly);
    for(i=0;i<size;i++){
        printf("0x%0*"PRIX64", ", bits/4, table[i]);
        if ((i&0x3)==0x3) printf("\n");
    }
}

Пользуясь Алгоритмом № умею синтезировать любые таблицы от CRC-8 до CRC-64. Для синтеза таблиц отраженного CRC использую аналогичный алгоритм. По сути таблица - это таблица умножения полинома на число 0,1..15, для 4битной таблицы и на число 0..255 для 8 битной таблицы. Таким образом, можно таблицу заменить на операцию полиномиального умножения с редуцированием по образующему полиному. В общем виде эти алгоритмы запишутся так:

Алгоритм № -- метод c умножением 8 бит.
[c4:c3:c2:c1] ::= [c3:c2:c1:0]⊕((c4⊕ Mi)•x8 mod P)

Как пользоваться редуцированием.
Один шаг над восьмибитными числами, обновляет значение CRC по байтам:

uint32_t CRC32_update_8 (uint32_t crc, uint8_t data){
    uint32_t val=data;
    val = BSWAP32(val); // перестановка байт
    crc = crc ^ val;
// Barrett's reduction
    uint32_t t = CL_MUL32H(crc>>24, 0x04D101DF)^(crc>>24);
    uint32_t v = CL_MUL32L(t      , 0x04C11DB7)^(crc<< 8);
    return v;
}

На самом деле многие алгоритмы и теоремы в конечном поле можно вовсе не доказывать. Доказательством может служить простой перебор чисел. Можно проверить утверждение для всех чисел в поле.
Я пишу серию алгоритмов для обработки по 8 бит, по 16 бит, по 32 по 64. И подбираю операции, отлаживаю для каждого шага. Потом из четырех полученных алгоримов создают пятый, который может работать со всеми длинами. Алгоритм можно вывести из математики, но на самом деле это две паралельные реальности, я смотрю на математику и пишу алгоритм, и наоборот проверяю какие-то действия в алгоритме и переписываю математику, которая этому соответствует. Что-то видно из математики что-то видно из алгоритма.

uint32_t CRC32_update_32 (uint32_t crc, uint32_t val){
 val = BSWAP32(val); // перестановка байт
 crc = crc ^ val;
// Barrett's reduction
 uint32_t t = CL_MUL32H(crc, 0x04D101DF)^(crc);
 uint32_t v = CL_MUL32L(t  , 0x04C11DB7);
 return v;
}

Вот этот вариант алгоритма проходит 32 бита (4 байта) за шаг. Ускорились в 4 раза!

// Алгоритм №. Расчет сдвиговых констант (xt mod P):

uint64_t xt_mod_P(uint64_t poly, int t, int bits)
{
    uint64_t v = poly;
    do {
        if (v & MSB){ // старший бит проверяем
          v = (v<<1) ^ poly;
        } else {
          v = (v<<1);
        }
    } while (--t);
    return v;
}

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

/* Алгоритм №. Расчет отраженных сдвиговых констант (xt mod P):
  poly - отраженный полином 33 бита.
*/
uint64_t xt_mod_P_reflect(uint64_t poly, int t, int bits)
{
    uint64_t v = poly;
    do {
 if (v & 1){
     v = (v>>1) ^ poly;
 } else {
     v = (v>>1);
 }
    } while (--t);
    return v;
}

Кроме того, хочу предложить алгоритм расчета сдвиговой константы для сдвига вправо, x^-t mod P(x). Этот вариант привожу для полноты картины.

// Алгоритм №. расчет сдвиговых констант, сдвиг вправо
uint64_t xt_mod_P_neg(uint64_t poly, int t, int bits){
    uint64_t v=1;
    do {
        if (v&1)
            v = (v>>1) ^ (1ULL<<(bits-1)|poly>>1);
        else
            v = (v>>1);
    } while (--t);
    return v;
}

После этого все операции сдвига я начал оценивать именно с этой позиции, каждое умножение в поле это сдвиги и набор сдвигов можно представить через операции умножения. Каждый байт сообщения (слово, двойное слово) мы подвергаем сдвигу в поле и в конце должны просуммировать все эти значения. Мы должны решить, как именно разбивать сообщение на группы, на байты и слова. В конечном итоге каждый байт исходного сообщения подвергается своему сдвигу и суммируется с результатом.
Еще одно интересное свойство, которе хотелось бы добавить для полноты изложения - можно предложить обратную операцию, отмена действия UNDO, операция CRC обратима. Можно предложить операцию дописывания шапки сообщения, независимо от тела. Сначала расчитывается CRC от тела сообщения, потом досчитывается CRC от шапки, и в итоге можно вычислить CRC от всего сообщения. И наконец, можно предложить обратное преобразование, которое находит ошибочный бит в сообщении -- восстановление сообщения. Ну и совсем хакерская мысль, можно предложить математическое преобразование, которое меняет сообщение, но сохраняет CRC неизменным. Можно переписать сообщение и подогнать значение CRC, добавив магию в соль.

Итак, начнем с голой идеи. Мне сложно сразу написать правильный, работающий алгоритм, практически невозможно отладить, отладить можно только каждое отдельное действие. Поэтому я создаю алгоритм по шагам. На первом шаге я делаю всего одну опрерацию - редуцирование, - отлаживаю ее.
На втором шаге делаю две операции: редуцирование и сдвиг в поле
M1•(x³² mod P) mod P.
Потом, отлаживаю алгоритм редуцирования и два сдвига,
M2•(x⁶⁴ mod P)) ⊕ M1•(x³² mod P)) mod P.
Затем - четыре сдвига,
M4•(x¹²⁸ mod P)) ⊕ M3•(x⁹⁶ mod P) ⊕ M2•(x⁶⁴ mod P) ⊕ M1•(x³² mod P) mod P.
Потом делаю алгоритм, который обрабатывает все сообщение целиком с циклическим схлопыванием и финальным редуцированием.

 Алгоритм №. Векторный алгоритм рассчета CRC-32
Шаг 1. загрузить вектор с переворотом 
 [c4:c3:c2:c1] := [M4:M3:M2:M1] ⊕ [crc:0:0:0]
Шаг 2. Цикл для каждого фрагмента 128-бит
    [A3:A2:A1] := [c4:c3]•(x192 mod P) ⊕ [c2:c1]•(x128 mod P)
 [c4:c3:c2:c1] := [ 0:A3:A2:A1] ⊕ [M4:M3:M2:M1]
Шаг 3. Сдвиги,
 [c3:c2:c1] := [с4:c3]•(x96 mod P) ⊕ [с2:c1]•(x32 mod P) 
    [c2:c1] := [c2:c1] ⊕ с3•(x64 mod P)
Шаг 4. Редуцирование 
    [T2:T1] := c2•Ux
    [c2:c1] := T2•Px ⊕ [c2:c1]
Возвращает: crc = c1

1. Вычисляем, например, CRC32 по словам 32 бита. Алгоритм №.
2. Пишу преобразование двойным словам 64 биата

uint32_t CRC32_update_64 (uint32_t crc, uint8_t* data){
 c = BSWAP64(data); // перестановка байт
 c = c ^ crc;       
// схлопывание [c2:c1]•x32 mod P(x) (96 бит в 64 бита)
 с = CL_MUL32(c[1], xt_mod_P(Px, 64))
   ^ CL_MUL32(c[0], xt_mod_P(Px, 32));
// Barrett's reduction
 uint32_t t = CL_MUL32H(crc, Ux)^(crc);
 uint32_t v = CL_MUL32L(t  , Px);
 return v;
}

// Алгоритм №. 
uint32_t    CRC32_update_128(uint32_t crc, uint8_t *data, int len) {
   uint64x2_t c;
   c = BSWAP128(data); data+=16;
   c[1] ^= (uint64_t)crc<<32;
   int blocks = (len+15 >> 4);
   while (--blocks>0) {
      c = CL_MUL64(c[1], K3) 
        ^ CL_MUL64(c[0], K4)
        ^ BSWAP128(data); data+=16;
   }
// схлопывание [c4:c3:c2:c1]•x32 mod P(x) (128бит в 64 бита)
   c = CL_MUL32(c[1]>>32, K4) 
     ^ CL_MUL32(c[1]>> 0, K5)
     ^ CL_MUL32(c[0]>>32, K6)
     ^ CL_MUL32(c[0]>> 0, K7);
// Редуцирование 64 бит в 32 бит
   uint32_t t = CL_MUL32H(c[0]>>32, Ux);
   crc = CL_MUL32L(t, Px) ^ c[0];
   return crc;
}

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

Алгоритм №. Векторный алгоритм рассчета CRC-32, 256 бит
Шаг 1. загрузить вектор с переворотом 
 [c8:c7:c6:c5:c4:c3:c2:c1] := [M8:M7:M6:M5:M4:M3:M2:M1] ⊕ [crc:0:0:0:0:0:0:0]
Шаг 2. Цикл для каждого фрагмента 256-бит
 [A3:A2:A1] := [c4:c3]•(x256+64 mod P) ⊕ [c2:c1]•(x256 mod P)
 [A7:A6:A5] := [c8:c7]•(x256+64 mod P) ⊕ [c6:c5]•(x256 mod P)
 [c8:c7:c6:c5:c4:c3:c2:c1] := [ 0:A7:A6:A5: 0:A3:A2:A1] ⊕ [M8:M7:M6:M5:M4:M3:M2:M1]
Шаг 3. Редуцирование 256 в 128
 [c4:c3:c2:c1] := [c4:c3:c2:c1]  
    ⊕ [c8:c7]•(x128+64 mod P) ⊕ [c6:c5]•(x128 mod P)
Шаг 4. если есть фрагмент 128-бит,
    [A3:A2:A1] := [c4:c3]•(x192 mod P) ⊕ [c2:c1]•(x128 mod P)
 [c4:c3:c2:c1] := [ 0:A3:A2:A1] ⊕ [M4:M3:M2:M1]
Шаг 5. Сдвиги,
 [c2:c1] := c4•(x128 mod P) ⊕ c3•(x96 mod P)
    ⊕ c2•(x64 mod P) ⊕ c1•(x32 mod P) 
Шаг 6. Редуцирование 
 [T2:T1] := c2•Ux
 [c2:c1] := T2•Px ⊕ [c2:c1]
Возвращает: crc = c1

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

// Алгоритм №. 
uint32_t CRC32K_update_8 (uint32_t crc, uint8_t *data){
 uint32_t val;
 val = data[0];
 crc = crc ^ val;
 poly64x2_t c = {crc};
 c = SLL128U(c,32+24);// сдвиг влево числа 128 бит
// Barrett's reduction
    poly64x2_t t;
    t  = CL_MUL128(c, (poly64x2_t){Ux,Px}, 0x00);
    c ^= CL_MUL128(t, (poly64x2_t){Ux,Px}, 0x10);
    return c[1];
}

// Алгоритм №. 
uint32_t CRC32K_update_16(uint32_t crc, uint8_t *data){
 uint32_t val;
 val = *(uint16_t*)data;
 crc = crc ^ val;
// подбираю связь между величиной сдвига и шагом алгоритма 
 poly64x2_t c = {crc};
 c = SLL128U(c,32+16);
// Barrett's reduction
    poly64x2_t t;
    t  = CL_MUL128(c, (poly64x2_t){Ux, Px}, 0x00);
    c ^= CL_MUL128(t, (poly64x2_t){Ux, Px}, 0x10);
    return c[1];
}

// Алгоритм №. 
uint32_t CRC32K_update_24(uint32_t crc, uint8_t *data){
 uint32_t val=0;
 __builtin_memcpy(&val, data,3);
 crc = crc ^ val;
// подбираю связь между величиной сдвига и шагом алгоритма 
 poly64x2_t c = {crc};
 c = SLL128U(c,32+8);
// Barrett's reduction
    poly64x2_t t;
    t  = CL_MUL128(c, (poly64x2_t){Ux, Px}, 0x00);
    c ^= CL_MUL128(t, (poly64x2_t){Ux, Px}, 0x10);
    return c[1];
}

// Алгоритм №. 
uint32_t CRC32K_update_32(uint32_t crc, uint8_t *data){
 uint32_t val;
 val = *(uint32_t*)data;
 crc = crc ^ val;
//  сдвиги заменяем на умножение 
 poly64x2_t c = {0,(uint64_t)crc<<32};
 c = CL_MUL128(c& 0xFFFFFFFF00000000ULL, (poly64x2_t){K7<<1}, 0x01);
// Barrett's reduction
 poly64x2_t t;
 t  = CL_MUL128(c, (poly64x2_t){Ux, Px}, 0x00);// {Ux, Px}
 c ^= CL_MUL128(t, (poly64x2_t){Ux, Px}, 0x10);
 return c[1];
}

// Алгоритм №. 
uint32_t CRC32K_update_128(uint32_t crc, uint8_t *data){
 poly64x2_t c = (poly64x2_t)LOAD128U(data);
 c[0]^=crc;
// сдвиги
 c = CL_MUL128(c<<32, (poly64x2_t){K4<<1}, 0x00)
   ^ CL_MUL128(c& 0xFFFFFFFF00000000ULL, (poly64x2_t){K5<<1}, 0x00)
   ^ CL_MUL128(c<<32, (poly64x2_t){K6<<1}, 0x01)
   ^ CL_MUL128(c& 0xFFFFFFFF00000000ULL, (poly64x2_t){K7<<1}, 0x01);
// Barrett's reduction
 poly64x2_t t;
 t  = CL_MUL128(c, (poly64x2_t){Ux, Px}, 0x00);// {Ux, Px}
 c ^= CL_MUL128(t, (poly64x2_t){Ux, Px}, 0x10);
 return c[1];
}

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

// Алгоритм №. 
uint32_t CRC32K_update_N(uint32_t crc, uint8_t *data, int len){
 poly64x2_t c = (poly64x2_t)LOAD128U(data);data+=16;
 c[0]^=crc;
// Циклическое схлопывание folding
 int blocks = (len+15 >> 4);
 while (--blocks>0){
  c = CL_MUL128(c, (poly64x2_t){K3<<1}, 0x00)
    ^ CL_MUL128(c, (poly64x2_t){K4<<1}, 0x01);
  c = SLL128U(c, 32) ^ (poly64x2_t)LOAD128U(data); data+=16;
 }
// Сдвиги редуцирование
 c = CL_MUL128(c<<32, (poly64x2_t){K4<<1}, 0x00)
   ^ CL_MUL128(c& 0xFFFFFFFF00000000ULL, (poly64x2_t){K5<<1}, 0x00)
   ^ CL_MUL128(c<<32, (poly64x2_t){K6<<1}, 0x01)
   ^ CL_MUL128(c& 0xFFFFFFFF00000000ULL, (poly64x2_t){K7<<1}, 0x01);
// Barrett's reduction
 poly64x2_t t;
 t  = CL_MUL128(c, (poly64x2_t){Ux, Px}, 0x00);// {Ux, Px}
 c ^= CL_MUL128(t, (poly64x2_t){Ux, Px}, 0x10);
 return c[1];
}

Что хотелось бы доработать в алгоритме. Убрать все сдвиги. Например, если у меня получается что умножается число со сдвигом, то из двух операций можно оставить одну - умножение. Везде где мне вздумалось при выводе алгоритма написать CL_MUL(a, K<<n), это можно решить иначае CL_MUL(a, K') пересчитать константу. Если константа расчиыватеся, как x^s mod P(x) Надо рассчитать константу x^s-1 mod P(x)

 Алгоритм №. Векторный алгоритм CRC-32B отраженный порядок 
Шаг 1. загрузить вектор M, 32x4 =128 бит
 [c4:c3:c2:c1] := [M4:M3:M2:M1] ⊕ [0:0:0:crc]
Шаг 2. Цикл для каждого фрагмента 128-бит
    [A3:A2:A1] := [c2:c1]•(x191 mod P) ⊕ [c4:c3]•(x127 mod P)
 [c4:c3:c2:c1] := [A3:A2:A1:0] ⊕ [M4:M3:M2:M1]
Шаг 3. Сдвиги,
 [c2:c1] := 
      с1•(x127 mod P)
    ⊕ с2•(x95 mod P)
    ⊕ с3•(x63 mod P)
    ⊕ с4•(x31 mod P)
Шаг 4. Редуцирование 
    [T2:T1] := c1•Ux
    [c2:c1] := T1•Px ⊕ [c2:c1]
Возвращает: crc = c2

В данном Алгоритме № есть еще один сдвиг, лишняя операция, которую можно заменить на умножение. Число A покомпонентно в цикле сдвигается на 32 бита. Этот сдвиг тоже можно заменить путем коррекции коэффициентов K3=X^191-32 и K4=X^127-32.

 Алгоритм №. Векторный алгоритм CRC-32B отраженный порядок 
Шаг 1. загрузить вектор M, 32x4 =128 бит
 [c4:c3:c2:c1] := [M4:M3:M2:M1] ⊕ [0:0:0:crc]
Шаг 2. Цикл для каждого фрагмента 128-бит
    [A3:A2:A1] := [c2:c1]•(x159 mod P) ⊕ [c4:c3]•(x95 mod P)
 [c4:c3:c2:c1] := [0:A3:A2:A1] ⊕ [M4:M3:M2:M1]
Шаг 3. Сдвиги,
 [c2:c1] := 
      с1•(x127 mod P)
    ⊕ с2•(x95 mod P)
    ⊕ с3•(x63 mod P)
    ⊕ с4•(x31 mod P)
Шаг 4. Редуцирование 
    [T2:T1] := c1•Ux
    [c2:c1] := T1•Px ⊕ [c2:c1]
Возвращает: crc = c2

Привожу результат практический:

CRC-32/BZIP2 poly=0x04c11db7
init=0xffffffff refin=false refout=false xorout=0xffffffff check=0xfc891918
Px = 0x104C11DB7
Ux = 0x104D101DF
K7 = xt_mod_P(Px,  32) = 0x04C11DB7
K6 = xt_mod_P(Px,  64) = 0x490D678D
K5 = xt_mod_P(Px,  96) = 0xF200AA66
K4 = xt_mod_P(Px, 128) = 0xE8A45605
K3 = xt_mod_P(Px, 192) = 0xC5B9CD4C
K2 = xt_mod_P(Px, 512) = 0xE6228B11
K1 = xt_mod_P(Px, 576) = 0x8833794C

Таблица умножения и константы:

CRC-32B poly=0x04c11db7
init=0xffffffff refin=true refout=true xorout=0xffffffff check=0xCBF43926

POLY=0xEDB88320
Таблица умножения:
0x00000000, 0x1DB71064, 0x3B6E20C8, 0x26D930AC,
0x76DC4190, 0x6B6B51F4, 0x4DB26158, 0x5005713C,
0xEDB88320, 0xF00F9344, 0xD6D6A3E8, 0xCB61B38C,
0x9B64C2B0, 0x86D3D2D4, 0xA00AE278, 0xBDBDF21C,

BarrettR u = x^64/P(x) Ur=0x1F7011641 Pr=0x1DB710641
[1/p]*p(x) = 163CD612400000001
Test =CBF43926 ..ok

K7( 32-1) = 00000001
K6( 64-1) = B8BC6765
K5( 96-1) = CCAA009E
K4(128-1) = 9BA54C6F
K3(192-1) = 65673B46

CRC-32C (Castagnoli) poly=0x1edc6f41
init=0xffffffff refin=true refout=true xorout=0xffffffff check=0xE3069283

POLY=0x82F63B78
Таблица умножения:
0x00000000, 0x105EC76F, 0x20BD8EDE, 0x30E349B1,
0x417B1DBC, 0x5125DAD3, 0x61C69362, 0x7198540D,
0x82F63B78, 0x92A8FC17, 0xA24BB5A6, 0xB21572C9,
0xC38D26C4, 0xD3D3E1AB, 0xE330A81A, 0xF36E6F75,

BarrettR u = x^64/P(x) Ur=0xDEA713F1 Pr=0x105EC76F1
[1/p]*p(x) = 0DD45AAB800000001
Test =E3069283 ..ok

K7( 32-1) = 00000001
K6( 64-1) = DD45AAB8
K5( 96-1) = 493C7D27
K4(128-1) = 3171D430
K3(192-1) = 3743F7BD

CRC-32K/BACnet (Koopman) poly= 0x741b8cd7 (0xEB31D82E)
init=0xffffffff refin=true refout=true xorout=0xffffffff check=0x2D3DD0AE

POLY=0xEB31D82E
Таблица умножения:
0x00000000, 0x83CF0F3C, 0xD1FDAE25, 0x5232A119,
0x7598EC17, 0xF657E32B, 0xA4654232, 0x27AA4D0E,
0xEB31D82E, 0x68FED712, 0x3ACC760B, 0xB9037937,
0x9EA93439, 0x1D663B05, 0x4F549A1C, 0xCC9B9520,

BarrettR u = x^64/P(x) Ur=0x17D232CD Pr=0x1D663B05D
[1/p]*p(x) = 018C5564C00000001
Test =2D3DD0AE ..ok

K7( 32-1) = 00000001
K6( 64-1) = 18C5564C
K5( 96-1) = 9D65B2A5
K4(128-1) = 69F48E4D
K3(192-1) = CBB06D55

[31.12.2019]

Написал векторный алогоритм для CRC-16 с отраженным порядком.

Алгоритм №. Векторный алгоритм CRC-16B отраженный порядок
Шаг 1. загрузить вектор M, 16x8 =128 бит
 [c8:c7:c6:c5:c4:c3:c2:c1] := [M8:M7:M6:M5:M4:M3:M2:M1] ⊕ [0:0:0:0:0:0:0:crc]
Шаг 2. Цикл для каждого фрагмента 128-бит
          [A5:A4:A3:A2:A1] := [c4:c3:c2:c1]•(x191 mod P) ⊕ [c8:c7:c6:c5]•(x127 mod P)
 [c8:c7:c6:c5:c4:c3:c2:c1] := [A5:A4:A3:A2:A1:0:0:0] ⊕ [M8:M7:M6:M5:M4:M3:M2:M1]
Шаг 3. Сдвиги,
 [c2:c1] := 
      с1•(x127 mod P)
    ⊕ с2•(x111 mod P)
    ⊕ с3•(x95 mod P)
    ⊕ с4•(x79 mod P)
    ⊕ с5•(x63 mod P)
    ⊕ с6•(x47 mod P)
    ⊕ с7•(x31 mod P)
    ⊕ с8•(x15 mod P)
Шаг 4. Редуцирование 
    [T2:T1] := c1•Ur
    [c2:c1] := T1•Pr ⊕ [c2:c1]
Возвращает: crc = c2

Тут сообщение и вектор (с) разбиты на элементы по 16 бит. Общая идея ясна алгоритм подобен Алгоритму № для CRC-32B.
Привожу наборы констант.

CRC-16/X25 poly=0x1021
init=0xffff refin=true refout=true xorout=0xffff check=0x906e
POLY=0x8408
Таблица умножения:
0x0000, 0x1081, 0x2102, 0x3183,
0x4204, 0x5285, 0x6306, 0x7387,
0x8408, 0x9489, 0xA50A, 0xB58B,
0xC60C, 0xD68D, 0xE70E, 0xF78F,

BarrettR u = x^16/P(x) Ur=0x1911 Pr=0x10811
[1/p]*p = 00000000019D80001
Test =906E ..ok

K9( 15) = 0001
K8( 31) = 19D8
K7( 47) = 1CBB
K6( 63) = 042B
K5( 79) = 81BF
K4( 95) = 2C27
K3(111) = 8555
K2(127) = 7EEA
K1(191) = A95D

CRC-16/MODBUS poly=0x8005
init=0xffff refin=true refout=true xorout=0x0000 check=0x4b37

POLY=0xA001
Таблица умножения:
0x0000, 0xCC01, 0xD801, 0x1400,
0xF001, 0x3C00, 0x2800, 0xE401,
0xA001, 0x6C00, 0x7800, 0xB401,
0x5000, 0x9C01, 0x8801, 0x4400,

BarrettR u = x^16/P(x) Ur=0x1BFFF Pr=0x14003
[1/p]*p = 000000001D0020001
Test =4B37 ..ok

K9( 15) = 0001
K8( 31) = 9001
K7( 47) = FC01
K6( 63) = D101
K5( 79) = CCC1
K4( 95) = C551
K3(111) = C3FD
K2(127) = C100
K1(191) = CCD0

[01.01.2020]

Зачем публиковать, если нет ни одного ценителя. Ау! Ценители отзовитесь. Нужно публиковать алгоритм?

Алгоритм №. Векторный алгоритм CRC-64/XZ отраженный порядок
Шаг 1. загрузить вектор M, 64x2 =128 бит
 [c2:c1] := [M2:M1] ⊕ [crc:0]
Шаг 2. Цикл для каждого фрагмента 128-бит
 [A2:A1] := c1•(x191 mod P) ⊕ c2•(x127 mod P)
 [c2:c1] := [A2:A1] ⊕ [M2:M1]
Шаг 3. Сдвиги,
 [c2:c1] := с1•(x127 mod P) ⊕ с2•(x63 mod P)
Шаг 4. Редуцирование 
    [T2:T1] := c1•Ur
    [c2:c1] := T1•Pr ⊕ [c2:c1]
Возвращает: crc = c2

Таблица параметров для алгоритма:

CRC-64/XZ poly=0x42F0E1EBA9EA3693 отраженный порядок
init=0xFFFFFFFFFFFFFFFF refin=true refout=true xorout=0xFFFFFFFFFFFFFFFF check=0x995DC9BBDF1939FA

POLY=0xC96C5795D7870F42
Таблица умножения:
0x0000000000000000, 0x7D9BA13851336649, 0xFB374270A266CC92, 0x86ACE348F355AADB,
0x64B62BCAEBC387A1, 0x192D8AF2BAF0E1E8, 0x9F8169BA49A54B33, 0xE21AC88218962D7A,
0xC96C5795D7870F42, 0xB4F7F6AD86B4690B, 0x325B15E575E1C3D0, 0x4FC0B4DD24D2A599,
0xADDA7C5F3C4488E3, 0xD041DD676D77EEAA, 0x56ED3E2F9E224471, 0x2B769F17CF112238,

BarrettR u = x^128/P(x) Ur=0x9C3E466C172963D5 Pr=0x192D8AF2BAF0E1E85 (65-бит)
[1/p]*p(x) = 48663A84688941C50000000000000001
Test =995DC9BBDF1939FA ..ok

K5( 63) = x63  mod P = 0000000000000001
K4(127) = x127 mod P = DABE95AFC7875F40
K3(191) = x191 mod P = E05DD497CA393AE4

Алгоритм №. Векторный алгоритм CRC-64/WE прямой порядок
Шаг 1. загрузить вектор M с разворотом байт, 64x2 =128 бит
 [c2:c1] := [M2:M1] ⊕ [0:crc]
Шаг 2. Цикл для каждого фрагмента 128-бит
 [A2:A1] := c2•(x192 mod P) ⊕ c1•(x128 mod P)
 [c2:c1] := [A2:A1] ⊕ [M2:M1]
Шаг 3. Сдвиги,
 [c2:c1] := с2•(x128 mod P) ⊕ с1•(x64 mod P)
Шаг 4. Редуцирование 
    [T2:T1] := c2•Ux
    [c2:c1] := T2•Px ⊕ [c2:c1]
Возвращает: crc = c1

Параметры алгоритма:

CRC-64/WE poly=0x42F0E1EBA9EA3693
init=0xFFFFFFFFFFFFFFFF refin=false refout=false xorout=0xFFFFFFFFFFFFFFFF check=0x62EC59E3F1A4F00A

POLY=0x42F0E1EBA9EA3693
Таблица умножения:
0x0000000000000000, 0x42F0E1EBA9EA3693, 0x85E1C3D753D46D26, 0xC711223CFA3E5BB5,
0x493366450E42ECDF, 0x0BC387AEA7A8DA4C, 0xCCD2A5925D9681F9, 0x8E224479F47CB76A,
0x9266CC8A1C85D9BE, 0xD0962D61B56FEF2D, 0x17870F5D4F51B498, 0x5577EEB6E6BB820B,
0xDB55AACF12C73561, 0x99A54B24BB2D03F2, 0x5EB4691841135847, 0x1C4488F3E8F96ED4,

Barrett u = x^128/P(x) Ux = 0x1578D29D06CC4F872 Px = 0x142F0E1EBA9EA3693
[1/p]*p(x) = 0x1000000000000000005F5C3C7EB52FAB6
Test =62EC59E3F1A4F00A ..ok

K5( 64) = 42F0E1EBA9EA3693
K4(128) = 05F5C3C7EB52FAB6
K3(192) = 4EB938A7D257740E

Вычисление CRC по обному биту

Не знаю какой в этом смылсл, если мы только что представили вычисление CRC по 128 бит за цикл. НО хочется вернуться к истокам. Рассмотрим под микроскопом каждый отдельный бит. Как выглядит операция CRC над битами. Прежде всего можно ввести операцию сдвиг-редуцирование, как основа для построения более сложных выражений:

XT(с) := MSB(c)==1? (с<<1)^Px: (c<<1);

Такая же операция для отраженного полинома выглядит так:

XT(с) := LSB(c)==1? (с>>1)^Pr: (c>>1);

У операции есть свойство группы, которое позволяет выполнять оптимизации:

XT(a ⊕ b) == XT(a) ⊕ XT(b)

Опрерация обновления CRC может быть представлена в виде:

Cn := XT(Cn-1 ⊕ dn<<(N-1)), где N -разрядность CRC

Можно упростить, раскрыв скобки

Cn := XT(Cn-1) ⊕ dn•Px, где N -разрядность CRC

Для двух бит получаем выражение:

Cn+1 := X2•Cn-1 ⊕ [dn+1,dn]•Px

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

Cn+3 := X4•Cn-1 ⊕ [dn+3..dn]•Px

ИЛИ

CRC = (CRC<<4) ^ PMUL(CRC>>(N-4) ^ Data, Px)

где PMUL -- операция полиномиального умножения и редуцирования. Умножение на полином Px - это по сути - сдвиг.

Cn+3 := X4•Cn-1 ⊕ [dn+3..dn]•XN

[05.01.2020]

Я увидел возможность для дальнейшей оптимизации алгоритма расчета CRC за счет увеличения разрядности редуцирования до 64бит. При этом можно сократить число сдвиговых операций(операций умножения) при редуцировании. Для этого необходимо пересчитать константы Барретта с увеличением разрядности до 64.

Алгоритм №. Векторный алгоритм CRC-32 отраженный порядок
Шаг 1. начальное значение
 [с4:c3:c2:c1] := [0:0:0:crc]
Шаг 2. Цикл для каждого фрагмента 128-бит, кроме последнего
 [с4:c3:c2:c1] := [с4:c3:c2:c1] ⊕ [M4:M3:M2:M1]
 [с4:c3:c2:c1] := [c2:c1]•(x31+128 mod P) ⊕ [c4:c3]•(x31+64 mod P)
Шаг 3. Загрузить последний фрагмент
 [с4:c3:c2:c1] := [с4:c3:c2:c1] ⊕ [M4:M3:M2:M1]
    [c3:c2:c1] := [c2:с1]•(x31+64 mod P) ⊕ [с4:c3]•(x31 mod P)
Шаг 4. Редуцирование 
 [T4:T3:T2:T1] := [c2:c1]•Ur
    [c3:c2:c1] := [T2:T1]•Pr ⊕ [c3:c2:c1]
Возвращает: crc = c3

В данном алгоритме используется модифицированная константа Барретт'а (Ur), рассчитанная для разрядности 64 бита. Привожу константы для алгоритмов с отраженным полиномом.

CRC-32B Zlib CRC-32/ISO-HDLC
poly=0x04C11DB7 refin=true refout=true check=0xCBF43926
Ur = 0xB4E5B025F7011641 Pr = 0x1DB710641

CRC-32C (Castagnoli)
poly=0x1EDC6F41 refin=true refout=true check=0xE3069283
Ur = 0x4869EC38DEA713F1 Pr = 0x105EC76F1

CRC-32K (Koopman) BACnet
poly=0x741B8CD7 init=0xffffffff refin=true refout=true xorout=0xffffffff check=0x2D3DD0AE
Ur = 0xC25DD01C17D232CD Pr = 0x1D663B05D

Те же идеи удалось перенести на CRC-16 и CRC-8 с отраженным порядком бит.

Алгоритм №. Векторный алгоритм CRC-16 отраженный порядок
Шаг 1. начальное значение
 [с8...c1] := [0:0:0:0:0:0:0:crc]
Шаг 2. Цикл для каждого фрагмента 128-бит, кроме последнего
 [с8...c1] := [с8...c1] ⊕ [M8...M1]
 [с8...c1] := [c4...c1]•(x15+128 mod P) ⊕ [c8...c5]•(x15+64 mod P)
Шаг 3. Загрузить последний фрагмент
 [с8...c1] := [с8...c1] ⊕ [M8...M1]
 [c5...c1] := [c4...с1]•(x15+64 mod P) ⊕ [с8...c5]•(x15 mod P)
Шаг 4. Редуцирование 
 [T8...T1] := [c4...c1]•Ur
 [c5...c1] := [T4...T1]•Pr ⊕ [c5...c1]
Возвращает: crc = c5

Алгоритм №. Векторный алгоритм CRC-8 отраженный порядок
Шаг 1. начальное значение
 [с16...c1] := [0...0:crc]
Шаг 2. Цикл для каждого фрагмента 128-бит, кроме последнего
 [с16...c1] := [с16...c1] ⊕ [M16...M1]
 [с16...c1] := [c8...c1]•(x7+128 mod P) ⊕ [c16...c9]•(x7+64 mod P)
Шаг 3. Загрузить последний фрагмент
 [с16...c1] := [с16...c1] ⊕ [M16...M1]
  [c9...c1] := [c8...с1]•(x7+64 mod P) ⊕ [с16...c9]•(x7 mod P)
Шаг 4. Редуцирование 
 [T16...T1] := [c8...c1]•Ur
  [c9...c1] := [T8...T1]•Pr ⊕ [c9...c1]
Возвращает: crc = c9

Параметры алгоритма CRC-8 с отражением

CRC-8/BAC POLY=0x81 refin=true refout=true
Ur = 0x808182878899AAFF(64 бита) Pr=0x103 (9 бит)
X7 mod P = 0x01
X71 mod P = 0x81
X135 mod P = 0xC1

Алгоритм №. Векторный алгоритм CRC-64 отраженный порядок
Шаг 1. начальное значение
 [с2:c1] := [0:crc]
Шаг 2. Цикл для каждого фрагмента 128-бит, кроме последнего
 [с2:c1] := [с2:c1] ⊕ [M2:M1]
 [с2:c1] := c1•K1 ⊕ c2•K2
Шаг 3. Загрузить последний фрагмент
 [с2:c1] := [с2:c1] ⊕ [M2:M1]
 [c2:c1] := с1•K3 ⊕ с2•K4
Шаг 4. Редуцирование 
 [T2:T1] := c1•Ur
 [c2:c1] := T1•Pr ⊕ [c2:c1]
Возвращает: crc = c2

Алгоритм CRC-8, CRC-16, CRC-32 в точности такие же как и для разрядности CRC-64, отличие только в значениях сдвиговых констант Kn и {Ur,Pr}, и число разрядов выходного значения (N). Таким образом алгоритм получается универсальный для всех вариантов CRC с отраженным порядком. НО это еще не финал. Дело в том, что представленный алгоритм работает только с входными данными, длиной кратной 128бит. А мы хотим сделать универсальный алгоритм, чтобы он работал с любыми длинами. Идея такая: можно на шаге 3. загружать последний фрагмент не кратный 128 бит и выполнять соответствующий его длине сдвиг. Для этого вносим соответствующее изменение в алгоритм и просчитываем соответствующие сдвиги, получаем таблицы K3 и K4 из 16 сдвиговых констант.

Алгоритм №. Векторный алгоритм CRC-64 отраженный порядок
Шаг 1. начальное значение
 [с2:c1] := [0:crc]
Шаг 2. Цикл для каждого фрагмента 128-бит, кроме последнего
 [с2:c1] := [с2:c1] ⊕ [M2:M1]
 [с2:c1] := c1•K1 ⊕ c2•K2
Шаг 3. Загрузить последний фрагмент k = (length & 15)
 [с2:c1] := [с2:c1] ⊕ [M2:M1]
 [c2:c1] := с1•K3(k) ⊕ с2•K4(k)
Шаг 4. Редуцирование 
 [T2:T1] := c1•Ur
 [c2:c1] := T1•Pr ⊕ [c2:c1]
Возвращает: crc = c2

Вот теперь моя работа закончена. Я изобрел свой собственный способ расчета CRC! Самый быстрый, самый крутой!

Наборы констант для CRC-16/X-25 BACnet с отраженным порядком
poly=0x1021 init=0xffff refin=true refout=true xorout=0xffff check=0x906e

K1, K2:{0x8E10, 0x81BF} // x^{143}, x^{79}
K3, K4:
[ 1] = {0x4EA8, 0x97B7},// x^{-41},x^{-105}
[ 2] = {0x290C, 0xC1A3},// x^{-33},x^{-97}
[ 3] = {0xCA45, 0x9750},// x^{-25},x^{-89}
[ 4] = {0x1563, 0x5212},// x^{-17},x^{-81}
[ 5] = {0x5188, 0x33C1},// x^{-9}, x^{-73}
[ 6] = {0x0811, 0xD7B6},// x^{-1}, x^{-65}
[ 7] = {0x0100, 0xD06A},// x^{ 7}, x^{-57}
[ 8] = {0x0001, 0xCC8C},// x^{15}, x^{-49}
[ 9] = {0x1189, 0x4EA8},// x^{23}, x^{-41}
[10] = {0x19D8, 0x290C},// x^{31}, x^{-33}
[11] = {0x5ADC, 0xCA45},// x^{39}, x^{-25}
[12] = {0x1CBB, 0x1563},// x^{47}, x^{-17}
[13] = {0x0B44, 0x5188},// x^{55}, x^{-9}
[14] = {0x042B, 0x0811},// x^{63}, x^{-1}
[15] = {0x9FD5, 0x0100},// x^{71}, x^{ 7}
[ 0] = {0x81BF, 0x0001},// x^{79}, x^{15}
Ur=0x859B040B1C581911 Pr=0x10811

[24.04.2021]

Нашел способ использовать один алгоритм с разрядностью 64 бита для расчета CRC-16/CRC-24/CRC-32/CRC-64. Фактически можно использовать алгоритм CRC-64 без изменения со своим набором констант для расчета контрольных сумм меньшей разрядности, например CRC-16. При этом используется правило:

 C = A mod P => C•K = A•K mod P•K

В нашем случае К - это свдиг влево на (64-16=48) бит.

Алгоритм №. Векторный алгоритм CRC-64 прямой порядок
Шаг 1. начальное значение
 [с2:c1] := [crc<<(64-N):0]
Шаг 2. Цикл для каждого фрагмента 128-бит, кроме последнего
 [с2:c1] := [с2:c1] ⊕ [M2:M1]
 [с2:c1] := c2•K2 ⊕ c1•K1
Шаг 3. Загрузить последний фрагмент k = (length & 15)
 [с2:c1] := [с2:c1] ⊕ [M2:M1]
 [c2:c1] := с2•K4(k) ⊕ с1•K3(k)
Шаг 4. Редуцирование 
 [T2:T1] := c2•Ux ⊕ [c2:c1]
 [c2:c1] := T2•Px ⊕ [c2:c1]
Возвращает: crc = c1>>(64-N)

Наборы констант для CRC-16
poly=0x1021 init=0xffff refin=false refout=false xorout=0xffff check=0xd64e

K1, K2:{0xAEFC, 0x650B},// x^{128}, x^{192}
K3, K4:
[ 1] = {0xCBF3000000000000, 0x2AE4000000000000},// x^{-104}, x^{-40}
[ 2] = {0x9B27000000000000, 0x6128000000000000},// x^{-96}, x^{-32}
[ 3] = {0x15D2000000000000, 0x5487000000000000},// x^{-88}, x^{-24}
[ 4] = {0x9094000000000000, 0x9D71000000000000},// x^{-80}, x^{-16}
[ 5] = {0x17B9000000000000, 0x2314000000000000},// x^{-72}, x^{-8}
[ 6] = {0xDBD6000000000000, 0x0001000000000000},// x^{-64}, x^{ 0}
[ 7] = {0xAC16000000000000, 0x0100000000000000},// x^{-56}, x^{ 8}
[ 8] = {0x6266000000000000, 0x1021000000000000},// x^{-48}, x^{16}
[ 9] = {0x2AE4000000000000, 0x3331000000000000},// x^{-40}, x^{24}
[10] = {0x6128000000000000, 0x3730000000000000},// x^{-32}, x^{32}
[11] = {0x5487000000000000, 0x76B4000000000000},// x^{-24}, x^{40}
[12] = {0x9D71000000000000, 0xAA51000000000000},// x^{-16}, x^{48}
[13] = {0x2314000000000000, 0x45A0000000000000},// x^{-8}, x^{56}
[14] = {0x0001000000000000, 0xB861000000000000},// x^{ 0}, x^{64}
[15] = {0x0100000000000000, 0x47D3000000000000},// x^{ 8}, x^{72}
[ 0] = {0x1021000000000000, 0xEB23000000000000},// x^{16}, x^{80}
Ur=0x11303471A041B343 Pr=0x1021000000000000

Наборы констант для CRC-24/OpenPGP
poly=0x864cfb init=0xb704ce refin=false refout=false xorout=0x000000 check=0x21cf02

K1, K2:{0x6243DA, 0xB22B31},// x^{128}, x^{192}
K3, K4:
[ 1] = {0x2471670000000000, 0x7190920000000000},// x^{-96}, x^{-32}
[ 2] = {0xEE50080000000000, 0xC6E2490000000000},// x^{-88}, x^{-24}
[ 3] = {0xCE8F4A0000000000, 0xD19E9A0000000000},// x^{-80}, x^{-16}
[ 4] = {0x1D1CA30000000000, 0xF77C040000000000},// x^{-72}, x^{-8}
[ 5] = {0x6DC6AA0000000000, 0x0000010000000000},// x^{-64}, x^{ 0}
[ 6] = {0x67F3180000000000, 0x0001000000000000},// x^{-56}, x^{ 8}
[ 7] = {0x79152C0000000000, 0x0100000000000000},// x^{-48}, x^{16}
[ 8] = {0xE2DD700000000000, 0x864CFB0000000000},// x^{-40}, x^{24}
[ 9] = {0x7190920000000000, 0x668F480000000000},// x^{-32}, x^{32}
[10] = {0xC6E2490000000000, 0x8309D70000000000},// x^{-24}, x^{40}
[11] = {0xD19E9A0000000000, 0x3609520000000000},// x^{-16}, x^{48}
[12] = {0xF77C040000000000, 0xD9FE8C0000000000},// x^{-8}, x^{56}
[13] = {0x0000010000000000, 0x36EB3D0000000000},// x^{ 0}, x^{64}
[14] = {0x0001000000000000, 0x3B918C0000000000},// x^{ 8}, x^{72}
[15] = {0x0100000000000000, 0xF50BAF0000000000},// x^{16}, x^{80}
[ 0] = {0x864CFB0000000000, 0xFD7E0C0000000000},// x^{24}, x^{88}
Ur=0xF845FE2493242DA4 Pr=0x864CFB0000000000

Наборы констант для CRC-32/BZIP2
poly=0x04c11db7 init=0xffffffff refin=false refout=false xorout=0xffffffff check=0xfc891918

K1, K2:{0xE8A45605, 0xC5B9CD4C},// x^{128}, x^{192}
K3, K4:
[ 1] = {0x052B9A0400000000, 0x876D81F800000000},// x^{-88}, x^{-24}
[ 2] = {0x3C5F6F6B00000000, 0x1ACA48EB00000000},// x^{-80}, x^{-16}
[ 3] = {0xBE519DF400000000, 0xA9D3E6A600000000},// x^{-72}, x^{-8}
[ 4] = {0xD02DD97400000000, 0x0000000100000000},// x^{-64}, x^{ 0}
[ 5] = {0x3C423FE900000000, 0x0000010000000000},// x^{-56}, x^{ 8}
[ 6] = {0xA3011FF400000000, 0x0001000000000000},// x^{-48}, x^{16}
[ 7] = {0xFD7384D700000000, 0x0100000000000000},// x^{-40}, x^{24}
[ 8] = {0xCBF1ACDA00000000, 0x04C11DB700000000},// x^{-32}, x^{32}
[ 9] = {0x876D81F800000000, 0xD219C1DC00000000},// x^{-24}, x^{40}
[10] = {0x1ACA48EB00000000, 0x01D8AC8700000000},// x^{-16}, x^{48}
[11] = {0xA9D3E6A600000000, 0xDC6D9AB700000000},// x^{-8}, x^{56}
[12] = {0x0000000100000000, 0x490D678D00000000},// x^{ 0}, x^{64}
[13] = {0x0000010000000000, 0x1B280D7800000000},// x^{ 8}, x^{72}
[14] = {0x0001000000000000, 0x4F57681100000000},// x^{16}, x^{80}
[15] = {0x0100000000000000, 0x5BA1DCCA00000000},// x^{24}, x^{88}
[ 0] = {0x04C11DB700000000, 0xF200AA6600000000},// x^{32}, x^{96}
Ur=0x04D101DF481B4E5A Pr=0x04C11DB700000000

MGM -- (Multilinear Galois Mode) векторная реализация

Вышли новые рекомендации к Российской национальной криптографии, режим шифрования MGM, который обещают использовать как основной и пожалуй единственный, если не считать рекомендации к протоколу CRISP. Режим MGM использует полиномиальное умножение в поле Галуа, по полиному x¹²⁸ + x⁷ + x² + x¹ + 1

Опишу порядок создания реализации алгоритма MGM. До начала работы надо иметь работающую реализацию алгоритма зашифровывания "Кузнечик" и тестовые вектора от MGM. Тестовые вектора нашел в draft стандарте IETF [Multilinear Galois Mode (MGM)]

Нужны глубокие познания в области умножения полиномов. Мы хотим использовать векторные инструкции полиномиального умножения, умножения без переносов.
На архитектуре Intel обещают операцию pclmul -- умножение без переноса двух чисел 64 бит, результат. Эта инструкция присутствует практически во всех новых процессорах Intel. И ожидаем появления векторной инструкции pclmul, которая позволяет производить одновременно 4 умножения 64бит на векторе 512 бит (AVX512+VPCLMULQDQ). Эта инструкция будет доступна на серверах Ice Lake в следующем году (2020). Но, алгоритмы на ее основе можно отладить уже сегодня с использованием эмулятора платформы Intel SDE. Я скачал SDE и убедился, что на нем можно эмулировать конкретные семейства процессоров Intel. Я использую компилятор GCC 9.2 под MSYS64, который хорошо осведомлен о новых процессорных инструкциях Intel и поддерживает целевую платформу Ice Lake.

Я не думаю, что векторные инструкции сильно ускорят процесс, потому что Intel обещает только один или два блока AVX512, в то время как каждое ядро способно выполнять четыре инструкции одновременно. Есть подозрение, что четыре инструкции 128 бит могут выполняться с такой же скоростью, что и одна инструкция AVX512 -- это пессимистический взгляд на вещи. Я не знаю, как устроен процессор на самом деле, но есть подозрение, что блоки AVX умеют объединяться в два блока AVX2 или в один блок AVX512. Возможно что в новых процессорах AVX2 -- неделимый блок, два блока AVX2 умеют объединяться в один блок AVX512. Если бы я разрабатывал процессоры Intel, я бы так и сделал, для экономии ресурсов. Почитал сообщения на форуме Intel, подтверждают, что быстродействие векторных инструкций AVX512 в два раза ниже чем быстродействие AVX2 на процессорах Xeon Scalable. В спецификации на процессор указывается число блоков AVX512. Подтверждают, что один блок AVX512 имеет такую же производительность на такт, как две инструкции AVX2. Но тактовая частота блоков AVX512 снижается от загрузки, она чуть ниже чем у ядра. Поэтому производительность получается даже ниже, чем на  инструкциях AVX2. Исходя из этих утверждений я хочу разработать алгоритм под векторные инструкции AVX2+VPCLMULQDQ, два умножения на инструкцию.


Начнем разработку с умножения без переносов.

Для начала решаем еще более простую задачу.
Пишем алгоритм умножения без переносов 64 бит на 64бит. с последующим редуцированием 128 бит в 64 бит. Редуцирование выполняется по полиному x⁶⁴ + x⁴ + x³ + x¹ + 1

uint64_t gfmul64(uint64_t a, uint64_t b)
{
 poly64x2_t r = CL_MUL128((poly64x2_t){a}, (poly64x2_t){b}, 0x00);
 uint64_t x1 = r[1];
 uint64_t x0 = r[0];
 x1 = x1 ^ x1>>63 ^ x1>>61 ^ x1>>60;
 return x0 ^ x1 ^ x1<<1 ^ x1<<3 ^ x1<<4;
}

Второй вариант алгоритма пишем с редуцированием на операциях умножения без переносов.

uint64_t gfmul64(uint64_t a, uint64_t b)
{
 poly64x2_t r = CL_MUL128((poly64x2_t){a}, (poly64x2_t){b}, 0x00);
// Редукция Барретт'а:
 poly64x2_t t = CL_MUL128(r, (poly64x2_t){0x1BULL}, 0x01) ^ r;
 r ^= CL_MUL128(t, (poly64x2_t){0x1BULL}, 0x01);
 return r[0];
}

В этом варианте я использую две константы 0x1B для редуцирования по методу Барретта. Константы высчитывал честно по алгоритму, который представлен в предыдущем сообщении. На самом деле первый вариант умножения gfmul64 тоже основан на редуцировании по методу Барретта только умножение разложено на операции сдвига и сложения в поле.
Проверяю лаконичность высказываний компилятора:

        vpclmulqdq  $0,  %xmm1, %xmm0, %xmm1
        vpclmulqdq  $17, %xmm2, %xmm1, %xmm0
        vpxor       %xmm1, %xmm0, %xmm0
        vpclmulqdq  $17, %xmm2, %xmm0, %xmm0
        vpxor       %xmm1, %xmm0, %xmm0

Я боюсь ошибиться. Поэтому сделал еще и третий вариант побитовое умножение.

uint64_t gfmul64_2(uint64_t a, uint64_t b)
{
    uint64_t r = 0;
    int i;
    for (i=0; i< 64; i++){
        if (b & (1ULL<<i)){
            r ^= a;
        }
        a = (a<<1) ^ (((int64_t)a>>63) & 0x1BULL);
    }
    return r;
}

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

Из популярной литературы я знаю алгоритм редуцирования на основе сдвигов [Carry-Less Multiplication and Its Usage for Computing the GCM Mode, алгоритм 4].

Теперь опишем умножение без переносов 128 бит на 128 бит. с последующим редуцированием 256 бит в 128 бит. Редуцирование выполняется по полиному x¹²⁸ + x⁷ + x² + x¹ + 1

Начнем с референсной реализации, образцовой, на которую ссылается Intel, Sun OpenSolaris, кто с кого списывал:

/// (Copyright 2009 Sun Microsystems, Inc.)
void gfmul(__m128i x_in_m, __m128i y_m, uint64_t *res)
{
   uint64_t R = { 0xe100000000000000ULL };
   struct aes_block z = { 0, 0 };
   struct aes_block v;
   uint64_t x;
   int i, j;
   uint64_t *x_in=(uint64_t*)&x_in_m;
   uint64_t *y=(uint64_t*)&y_m;
   v.a=y[1];
   v.b=y[0];
   for (j = 1; j>=0; j--) {
      x = x_in[j];
      for (i = 0; i < 64; i++, x <<= 1) {
         if (x & 0x8000000000000000ULL) {
            z.a ^= v.a;
            z.b ^= v.b;
         }
         if (v.b & 1ULL) {
             v.b = (v.a << 63)|(v.b >> 1);
             v.a = (v.a >> 1) ^ R;
         } else {
             v.b = (v.a << 63)|(v.b >> 1);
             v.a =  v.a >> 1;
         }
      }
   }
   res[0] = z.b;
   res[1] = z.a;
}

Этот алгоритм так себе описан. Он медленный, 128 циклов. Вроде люди пытались описать с использованием векторных типов и псевдо-функций векторных инструкций (wmmintrin), но что-то пошло не так. Я даже знаю, что пошло не так, компилятор Sun не имел никаких векторных расширений и это был черти какой год, GCC тогда тоже не умел их использовать, а Intel умел только через явное использование этих самых псевдо-функций. Пробую описать тоже самое, теми средствами, которые понятны современному компилятору GCC 9.+.
Надо заметить, что в GCM используется обратный порядок следования бит. А нам нужен прямой порядок.

Алгоритм 1 полиномиальное умножение 128x128 бит
Шаг 1: рассчитываем умножение без переноса для двух чисел: 
  A = [A1:A0]; B= [B1:B0];
Покомпонентный результат умножения будет:
  A0•B0 = [C1:C0], A1•B1 = [D1:D0], A0•B1 = [E1:E0], A1•B0 = [F1:F0]
Шаг 2: составляем результат 256-бит:
 [A1:A0]•[B1:B0] = [D1:F1⊕E1⊕D0:F0⊕E0⊕C1:C0]

На шаге 2 можно применить способ умножения Карацубы

Алгоритм 2 полиномиальное умножение 128x128 бит методом Карацубы
Шаг 1: Покомпонентный результат умножения будет:
  A0•B0 = [C1:C0], A1•B1 = [D1:D0], (A0⊕A1)•(B0⊕B1) = [E1:E0]
Шаг 2: составляем результат 256-бит:
 [A1:A0]•[B1:B0] = [D1:D0⊕D1⊕E1⊕C1:D0⊕E0⊕C0⊕C1:C0]

Если считать, что инструкция умножения без переносов выполняется за такт, за то же время, что и операция Исключающее ИЛИ, то экономии не будет. Это можно проверить.
Описываем умножение по Алгоритму 1:

/// Алгоритм 1 полиномиального умножения 128x128
poly64x4_t CL_MUL256(poly64x2_t p, poly64x2_t v)
{
    poly64x2_t q  = CL_MUL128(p, v, 0x10) 
                  ^ CL_MUL128(p, v, 0x01);// Карацуба отдыхает
    poly64x2_t r0 = CL_MUL128(p, v, 0x00) ^ SLL128(q, 64);
    poly64x2_t r1 = CL_MUL128(p, v, 0x11) ^ SRL128(q, 64);
    return (poly64x4_t){r0[0], r0[1], r1[0], r1[1]};
}

Тут я использую встроенные псевдо-функции, которые призваны генерировать векторные инструкции.

/// SRL128 - сдвиг числа 128бит вправо на указанное число бит
static inline poly64x2_t SRL128(v2du a, const int bits) {
    return (poly64x2_t)__builtin_ia32_psrldqi128((v2di)a, bits);
}

/// SRL128 - сдвиг числа 128бит влево на указанное число бит
static inline poly64x2_t SLL128(poly64x2_t a, const int bits) {
    return (poly64x2_t)__builtin_ia32_pslldqi128((v2di)a, bits);
}

/// CL_MUL128 - полиномиальное умножение 64х64 бит
static inline poly64x2_t CL_MUL128(poly64x2_t a, poly64x2_t b, const int c) {
    return (poly64x2_t)__builtin_ia32_pclmulqdq128 ((v2di)a,(v2di)b,c);
}

Из сатьи взял алгоритм со сдвигами для редуцирования по полиному x¹²⁸ + x⁷ + x² + x + 1

Алгоритм 4 редуцирование по полиному x128
Обозначим входные данные [X3:X2:X1:X0] с разбивкой числа 256 бит по 64 бит.
Шаг 1: Сдвигаем X3 на 63, 62 и 57-бит вправо.
  A = X3>>63, B = X3>>62, С = X3>>57
Шаг 2: D = X2 ⊕ A ⊕ B ⊕ C
Шаг 3: Выполняем свдиги над числам 128 бит
 [E1:E0] = [X3:D]<<7
 [F1:F0] = [X3:D]<<2
 [G1:G0] = [X3:D]<<1
Шаг 4: Сложение (XOR) над числми [E1:E0], [F1:F0], [G1:G0] и [X3:D]. 
 [H1:H0] = [X3 ⊕ E1 ⊕ F1 ⊕ G1 : D ⊕ E0 ⊕ F0 ⊕ G0]
Шаг 5: Возвращаем результат [X1 ⊕ H1: X0 ⊕ H0].

Пояснения к алгоритму 4. Мы можем сделать тоже самое с использованием умножения. Константа Барретта при использовании редуцирования будет равна самому полиному, сдвиги {0,1,2,7} соответствуют умножению на число 0x87, хвост от полинома (x⁷ + x² + x + 1).

Алгоритм 5 редуцирование по полиному x128
Обозначим входные данные [X3:X2:X1:X0] с разбивкой числа 256 бит на 64 бит.
Шаг 1: Умножаем старшую часть числа X на константу Барретта (0x87).
  [A3:A2:A1:A0] = [X3:X2]•0x87
Шаг 2: B = A ⊕ X
Шаг 3: Умножаем старшую часть числа B на полином (0x87)
  [C3:C2:C1:C0] = [B3:B2]•0x87
Шаг 4: D = C ⊕ X
Шаг 5: Возвращаем результат [D1:D0].

Благодаря равенству констант алгоритм можно упростить, сократив число умножений с четырех до трех.

Алгоритм 6 редуцирование по полиному x128
Обозначим входные данные [X3:X2:X1:X0] с разбивкой числа 256 бит на 64 бит.
Шаг 1: Умножаем старшую часть числа X на константу Барретта (0x87).
  [0:A2:A1:A0] = [X3:X2]•0x87
Шаг 2: Умножаем старшую часть числа B на полином (0x87)
  [0: 0: 0:C0] = [0 :A2]•0x87
Шаг 3: D = C ⊕ A ⊕ X
Шаг 4: Возвращаем результат [D1:D0].

Совместим алгоритм №6 с алгоритмом №1 умножения полиномов.

Алгоритм 7 редуцирование по полиному x128 + 0x87
Шаг 1:
  X2•0x87 = [A1:A0], X3•0x87 = [B1:B0]
Шаг 2:
  B1•0x87 = [C1:C0]
Шаг 3: Возвращаем [B0⊕A1⊕С1⊕X1:A0⊕C0⊕X0]

Можно заметить, что заменой умножения на сдвиги можно получить алгоритм №4.

// Алгоритм №7: Редуцирование по полиному x128 + 0x87
poly64x2_t gf128_reduction(poly64x2_t r0, poly64x2_t r1)
{
 const poly64x2_t Px ={0x87};
 poly64x2_t b  = CL_MUL128(r1,Px, 0x01);
 poly64x2_t a  = CL_MUL128(r1,Px, 0x00);
 poly64x2_t c  = CL_MUL128( b,Px, 0x01);
 return r0 ^ a ^ c ^ SLL128(b,64);
}

Проверяем результат компиляции:

// $ gcc -O3 -march=native -o - -S test.c
        vpclmulqdq      $1, %xmm2, %xmm0, %xmm3
        vpclmulqdq      $0, %xmm2, %xmm0, %xmm0
        vpclmulqdq      $1, %xmm2, %xmm3, %xmm2
        vpxor   %xmm0, %xmm2, %xmm2
        vpslldq $8, %xmm3, %xmm3
        vpxor   %xmm3, %xmm1, %xmm3
        vpxor   %xmm3, %xmm2, %xmm0

-- ничего лишнего.
Проверяем тестовые вектора, взял из статьи.

Умножение в конечном поле GF(2128), P(x) = x128 + x7 + x2 + x + 1
a = 0x7b5b54657374566563746f725d53475d
b = 0x48692853686179295b477565726f6e5d
gfmul128 (a, b) 
  = 0x040229a09a5ed12e7e4e10da323506d2

Результат. Для алгоритма MGM разработал эффективную реализацию умножения в конечном поле GF(2⁶⁴) и GF(2¹²⁸). Умножение состоит из двух частей - умножение без переноса и редуцирование. В цикле достаточно выполнять только умножение без переноса, а этап редуцирования можно вынести из цикла.

Несколько комментариев к статье

Операция PCLMUL -- тормозная. Она долго грузится, но исполняется за такт. На разных платформах она выполняется по-разному. Хотелось бы сюда вставить табличку, из которой видно насколько эффективно использовать эту инструкцию. Быстродействие складывается из двух параметров: Задержка и Исполнение. Инструкции могут долго загружаться, но выполняться одна за другой. Например, для новых серверов я ожидаю быстродействие 6 тактов задержка на исполнение и 1 такт на исполнение. Это значит, что группа инструкций PCLMUL не связанных между собой могут грузиться 6 тактов, и исполняться по 1-ому такту на команду. В некоторых случаях процессор выполняет инструкции в несколько потоков, тога пиковая производительность получается выше. Следить за тем чтобы команды эффективно грузились и переставлять порядок исполнения - это дело компилятора и процессора. НО зачастую уменьшив взаимосвязь между параметрами алгоритма можно добиться большей производительности.

Для понимания, когда стиот применять инструкции, а когда это совсем не эффективно, привожу таблицу по производительности операции PCLMUL на разных платформах. Таблицу нашел в итнернетах, источник - [Agner Fog. Lists of instruction latencies, throughputs..] Таблицы таймингов (задержки и производительность инструкций) для некторых моделей процессоров Intel представлены в
[#]Intel(c) Intel® 64 and IA-32 Architectures Optimization Reference Manual, приложение С.3 [Appendix C.3, “Latency and Throughput”].

Процессор      Задержка+обработка 
-- AMD -- 
Bulldozer     12 7
Piledriver    12 7
Steamroller   11 7
Excavator      5 5
Ryzen          4 2
Jaguar         3 1
-- Intel --
Silvermont    10 10 -- ужасно медленно!!
Nehalem       12 8
Ivy Bridge    14 8
Goldmont       6 4
Haswell        7 2
KnightsLanding 6 2
Skylake        7 1
SkylakeX       7 1
Broadwell      5 1

Можно утверждать, что в современных процессорах инструкция выполняется достаточно быстро, чтобы вытеснить другие решения. Все познается в сравнении, другие решения могут сроится на таблицах, операциях PXOR и PSHUFB. Например, на платформе Skylake операция PSHUFB выполняется за 1+1 такт, операция PXOR за 1+0.5 такта.

[22.12.2019]

Нашел возможность дальнейшей оптимизации Алгоритма редуцирования (№7), сократил число умножений до двух.

Алгоритм №8 редуцирование по полиному x128 + 0x87
Шаг 1:
  X3•0x87 = [B1:B0]
Шаг 2:
  (B1⊕X2)•0x87 = [D1:D0]
Шаг 3: Возвращаем [B0⊕D1⊕X1:D0⊕X0]

// Алгоритм №8: Редуцирование по полиному x128 + 0x87
poly64x2_t gf128_reduction(poly64x2_t r0, poly64x2_t r1)
{
 const poly64x2_t Px ={0x87};
 poly64x2_t b = CL_MUL128(r1,Px, 0x01) ^ SLL128(r1,64);
 poly64x2_t d = CL_MUL128(b, Px, 0x01) ^ SLL128(b, 64);
 return r0 ^ d;
}

Возвращаясь к режиму GCM, приведу мой вариант редуцирования для отраженного порядка следования бит.

// Алгоритм №9: Редуцирование по полиному x128 + 0x87
poly64x2_t gf128r_reduction(poly64x2_t r0, poly64x2_t r1)
{
 const poly64x2_t Px ={0x87};
 poly64x2_t b = CL_MUL128(r1,Px, 0x01) ^ SLL128(r1,64);
 poly64x2_t d = CL_MUL128(b, Px, 0x01) ^ SLL128(b, 64);
 return r0 ^ d;
}

Алгоритм №9 я сам вывел. Однако, до меня в примерно в такой форме его использовали в OpenSSL () -- этот факт обнаружился уже позже. Алгоритм в OpenSSL сложно распознать, поскольку он дан уже в ассемблерном коде. Я его распознал там по сигнатуре, схожая последовательность инструкций PCLMUL+PALIGNER, как в моем случае.

[23.12.2019]

Сравнил быстродействие варианта умножения методом Карацубы по Алгоритму №2 (три инструкции PCLMUL) и Алгоритма №1. Существенного выиграша не обнаружил. Думаю следует использовать вариант без Карацубы, по Алгоритму №1, на новых процессорах этот вариант будет быстрее.