пятница, 27 декабря 2019 г.

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) • xs ⊕ M0(x) = M1 • (xs mod P) ⊕ M0
Эту операцию называют "схлопывание" (folding), слов не хватает, слова "редуцирование" и "сложение" уже заняты и наполнены смыслом. Сразу предлагаем воспользоваться правилом: редуцирование от редуцирования на результат не влияет. Поэтому мы можем сначала выполнить редуцирование сомножителей, а потом применять операцию умножения.
Константы (xs 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 = x128 mod P
K3 = x192 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•(x32 mod P) mod P.
Потом, отлаживаю алгоритм редуцирования и два сдвига,
M2•(x64 mod P)) ⊕ M1•(x32 mod P)) mod P.
Затем - четыре сдвига,
M4•(x128 mod P)) ⊕ M3•(x96 mod P) ⊕ M2•(x64 mod P) ⊕ M1•(x32 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') пересчитать константу. Если константа расчиыватеся, как xs mod P(x) Надо рассчитать константу xs-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=X191-32 и K4=X127-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

четверг, 19 декабря 2019 г.

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

Вышли новые рекомендации к Российской национальной криптографии, режим шифрования MGM, который обещают использовать как основной и пожалуй единственный, если не считать рекомендации к протоколу CRISP. Режим MGM использует полиномиальное умножение в поле Галуа, по полиному x128 + x7 + x2 + x1 + 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 бит. Редуцирование выполняется по полиному x64 + x4 + x3 + x1 + 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 бит. Редуцирование выполняется по полиному x128 + x7 + x2 + x1 + 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);
}
Из сатьи взял алгоритм со сдвигами для редуцирования по полиному x128 + x7 + x2 + 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, хвост от полинома (x7 + x2 + 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(264) и GF(2128). Умножение состоит из двух частей - умножение без переноса и редуцирование. В цикле достаточно выполнять только умножение без переноса, а этап редуцирования можно вынести из цикла.
Несколько комментариев к статье
Операция 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, на новых процессорах этот вариант будет быстрее.

суббота, 30 ноября 2019 г.

Блочный шифр ГОСТ "Кузнечик" -- оптимизация

Речь идет об оптимизации алгоритма блочного шифра ГОСТ "Кузнечик". Я давно отложил задачу на сладкое, а тут время появилось ее доделать.
В основе шифра Кузнечик лежат две операции. 1) нелинейное биективное преобразование - табличная подстановка. 2) Что-то вполне линейное, в смысле алгебра линейная в конечном поле, на основе арифметики Галуа, с редуцированием по полиному 0x1С3.

Я увидел возможность преобразовать алгоритм, сделать его быстрее.
Итак структура алгоритма такова [ГОСТ 34.12-2015]:
E[K] = X[K10]LSX[K9]... LSX[K2]LSX[K1](a)
Я пропущу описание того, что сделано в первой серии. Я умею выражаться длинными векторами, для этого использую векторные расширения языка Си. Блочный шифр -- это преобразование вектора состоящего из 16 элементов разрядностью 8 бит. Для этого использую описание типа:

typedef uint8_t uint8x16_t __attribute__((__vector_size__(16)));

Алгоритм защифровывания выглядит так
uint8x16_t kuzn_encrypt(KuznCtx* ctx, const uint8x16_t a)
{
    uint8x16_t S = a ^ ctx->K[0];
    int i;
    for (i=0; i<9; i++){
        S = LS(S) ^ ctx->K[i+1];
    }
    return S;
}

Шифрование -- это вот такой алгоритм, внутри которого есть функция LS, всего девять циклов по развернутому ключу {K1,K2, ...K10}.
Зашифровывание -- это последовательное применение преобразований:
X[K1], S, L, X[K2], S, L, ... X[K9], S, L, X[K10]
Преобразование X[K](a) = K ⊕ a
Преобразование S -- это таблица подстановок -- "нелинейное биективное преобразование".

uint8x16_t LS(uint8x16_t a)
{
    int i;
    for(i=0; i< 16;i++) a[i] = sbox[a[i]]; // -- подстановка, 
    uint8x16_t v = R16(a); // -- то самое преобразование L
    return v;
}

Вот добрались до сути. КАК оптимизировать преобразование R16?
Преобразование R выполняется над векторами 16 элементов по 8 бит и содержит 16 циклов.
Кроме всяких пермутаций это преобразование содержит много операций умножения в конечном поле GF(2) с полиномом p(x) = x8 + x7 + x6 + x + 1, или Px = 0x1C3.
Путем долгих мучительных преобразований можно записать, что каждый элемент входного вектора зависит от каждого элемента выходного вектора. Всего 16 наборов констант для каждого элемента. vi = SUM(aj * LMTij) . LMT - это матрица коэффициентов, которую по ходу пришлось пересчитывать и транспонировать, чтобы преобразование стало выглядеть линейно.

uint8x16_t GMULv16(uint8x16_t L, uint8x16_t a);
uint8x16_t R16(uint8x16_t a)
{
    int i;
    uint8x16_t v = {0};
    for(i=0; i< 16;i++) {
        v ^= GMULv16(LMT[i], a[i]); // -- векторное умножение в поле Px = 0x1C3
    }
    return v;
}
Операйця GMULv16 выполняет умножение каждого элемента вектора с редуцированием по полиному 0x1C3.
На этом этапе я остановился в прошлый раз. Большим достижением казалость то, что алгоритм удалось сделать без ветвления и с ипользованием векторных инструкций. GMULv16 можно заменить на 16 таблиц по 256 значений по 16 байт каждое, 64кБайт. Но это не наш путь.

В цикле можно использовать умножение без переноса, а финальное редуцирование выполнять за пределами цикла.
vi = SUM((aj * LMTij) mod P(x)) =  SUM(aj * LMTij) mod P(x)

В результате получаем такой алгоритм:

uint8x16_t R16(uint8x16_t a)
{
    int i;
    uint16x16_t vh = {0};// -- 16 элементов разрядностью 16 бит. 256 бит.
    for(i=0; i< 16;i++) {
        vh ^= CL_MULv16(LMT[i], a[i]); // -- векторное умножение без переносов
    }
    uint8x16_t v = RR (vh, Px);// -- финальное редуцирование по полиному P(x),
    return v;
}
Размер таблицы - 16 элементов по 16 байт, 256 байт. В цикле получаем умножение 128 бит на 8 бит. Это умножение можно представить 4я инструкциями умножения полиномов в поле разрядностью 64 бит.

Может быть есть компромис между размером таблицы и безумным умножением?
Я так понимаю, что большие таблицы медленно обрабатываются, потому что в кеш память не влезают. Если мы делаем алгоритм для контроллера, то таблицы 65кБайт не лезут во флеш память, не лезут в память оперативную, сильно увеличивают объем кода. Обработка таблицы будет идти со скоростью работы памяти, а не со скоростью ядра процессора. Для встроенных приложений нужен компактный код. На нашей целевой архитектуре есть инструкция полиномиального умножения (ARM Cortex-A8+ NEON, ARMv8, Intel Core+AVX), но в контроллерах ARMv7e-m ничего такого нет, вынужденно используем таблицы. Нам понадобятся оба варианта алгоритма. Таблица не лезит в память, поэтому я сделал две таблицы, одна - это таблица умножения на числа 0,1,2,3...15, младшие биты числа. Вторая тоже таблица умножения на числа 0,16,32,64...256, - умножение на старшие биты числа.

uint8x16_t R16(uint8x16_t a)
{
    int i; 
    uint8x16_t v = {0}; 
    uint8x16_t a0 = a & 0xF,  a1 = (a >>4) & 0xF; 
    for(i=0; i< 16;i++) { 
        v ^= GMULv16_L[i][a0[i]] ^ GMULv16_H[i][a1[i]]; 
    }
    return v; 
}

В этом варианте умножения используются 32 таблицы по 16 значений по 16 байт, 8кБайт. Надо проверить, какой из двух вариантов быстрее: с подстановкой 8кБайт или с безумным умножением...
Проверил. Сильно зависит от процессора. На 3-м поколении Intel Core-i7 умножение без переносов работает крайне медленно, в пять раз медленнее чем табличный вариант. На 7-м поколении Core-i7 умножение работает в два раза медленнее чем табличный вариант. Разница есть. Табличный вариант позволяет получить скорость зашифровывания около 8 Мбайт/сек. Это все равно медленно.

[11.12.2019] Дальше я исследую возможность переноса алгоритма с умножением полиномов на векторные инструкции ARM NEON. У меня под рукой есть кросс-компилятор GNU ARM Embedded Toolchain, GCC. И есть под рукой плата BeagleBone, на которой могу запустить полученный код. Платформа содержит процессор ARM Cortex-A8 с архитектурой ARMv7-A и поддерживает инструкции NEON.
В системе команд NEON имеется инструкция полиномиального умножения VMULL.P8, которая оперирует с векторами poly8x8_t -- 8 элементов по 8бит. Инструкция производит результат - вектор 128 бит, 8 элементов по 16 бит. Подробнее см. [NEON™ Version: 1.0 Programmer’s Guide] на сайте ARM.

Мой код выглядит так:

#include <arm_neon.h>
uint8x16_t R16(uint8x16_t a)
{
    poly16x8_t v0 = {0};
    poly16x8_t v1 = {0};
    int i;
    for(i=0;i<16;i++){
        poly8x8_t  a8 = vdup_n_u8(sbox[a[i]]);// размножаем значение на все элементы
 poly8x16_t p8 = LMT[i];
        v0 ^= vmull_p8(a8, vget_low_p8 (p8));  // младшая часть вектора poly8x8
        v1 ^= vmull_p8(a8, vget_high_p8(p8));// старшая часть вектора poly8x8
    }
    /// редуцирование вынесли из цикла
    . . .
}
В коде использую псевдофункции, которые описаны в заголовке и в результате генерируют инструкции NEON. Исследую возможные комбинации ключей компилятора, чтобы задействовать инструкции NEON:
$ arm-eabi-gcc -print-multi-lib
Нахожу подходящую комбинацию ключей:
$ arm-eabi-gcc -march=armv7-a+fp -mcpu=cortex-a8 -mfpu=neon-vfpv4 \
  -mfloat-abi=hard -O3 -S -o - kuzn.c
Содержимое цикла преобразцется в набор инструкций

        vdup.8   d18, r3        -- размножить значение 8x8
        vld1.64 {d24-d25}, [r1:64]!  -- загрузить LMT 8x8x2
        vmull.p8 q13, d18, d24  -- умножить вектора 8x8 
        vmull.p8 q9,  d18, d25  -- умножить вектора 8x8
        veor     q8,  q8,  q13  -- исключающее ИЛИ 16x8
        veor     q11, q11, q9   -- исключающее ИЛИ 16x8
Ничего лишнего, одно умножение - одна инструкция процессора. Все инструкции используют вектора 8x16, умножение проивзодится над векторами 8x8, результат умножения имеет размерность 16x8. Не путать - умножение проиводится без переносов.
Можно ли считать этот результат хорошим. Не знаю. Объективно надо замерить скорость и сравнить результат. Код получился компактный, но не факт что инструкция умножения будет выполняться быстрее, чем обращение к памяти. Делаю 1Миллион циклов зашифровывания и замеряю время -- получаю результат = 2.5 секунды, на платформе beaglebone Texas Instrumens Sitara AM335x. Это получается скорость шифрования 6.4Мбайта/сек. Замеряю результат табличного варианта алгоритма зашифровывания = 4.2 секунды. Мы выиграли время и место. Дальше за счет разворачивания цила удалось повысить скорость зашифровывания при использовании умножения до 7.4Мбайт/сек. Это ХОРОШИЙ результат!

вторник, 19 ноября 2019 г.

Резюме 2019/2020

Ищу работу!

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

Люблю программировать. Умею писать операционные системы, умею писать криптографию, умею на время создавать безумное количество электроники для автоматизации, чтобы сразу работало.
Умею вести производство военной техники, умею разрабатывать военную технику с гарантией 7 лет и сроком службы до 40 лет.
Умею заниматься снабжением военных и гражданских проектов практически любых. Умею работать со спец. счетами и гос контрактами. Умею отбиваться от налоговых и прокурорских проверок, хотя это просто везение, научился выживать в условиях травли. Умею сохранять данные в облаках и на выделенных серваках, шифровать и нарезать мелко-мелко большие данные. Умею жонглировать серверами, когда их выносит ФСБ или МВД. Умею писать обработку протоколов ФТС. Умею писать блок чейны и архивы. Умею разбирать архивы документов XML со скоростью ... с любой скоростью. Имею опыт разработки корпоративных информационных систем на этих технологиях.

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

Почему-то не умею извлечь прибыль из этого.
Скучаю, ищу работу...Безумно устал от всей суеты. В отпуск хочу, в увольнительную. Хочу программы писать со скуки и воспитывать программистов. Достаточно одного.


Резюме 2011

Ищу работу!
Чего умею?
  • Умею быть преподавателем. Лет 5 отработал преподавателем электроники. Создал свой собственный курс по микросхемотехнике в Политехническом университете.
  • Могу управлять небольшим коллективом разработчиков.
  • Умею проектировать электронику: измерительную, цифровую и силовую. 
  • Умею писать встроенное программное обеспечение для ARM, 8051.
  • Умею производить спроектированную электронику. Умею организовать контрактное производство электроники штучное, мелко серийное и серийное.
  • Умею создавать сети для организаций, на 100-1000 рыл. Вообще-то умею создавать виртуальные транс-национальные корпорации не слезая со стула, вернее проектировать и настраивать аппаратуру и ПО для них. Умею в одиночку обслуживать 40 Кисок (CISCO) и десяток серверов. Хорошо знаком с GNU/Linux, FreeBSD, Oracle Solaris. Могу дописывать и исправлять код каждой из них, если понадобится.
  • Знаю толк в информационной безопасности. Могу настроить инфраструктуру организации, включая сервисы e-mail, web, vpn, ip телефонию и многое другое. 
  • Программировать люблю, на С. Умею проектировать все, что касается ядра операционной системы. Могу написать свою операционную систему с нуля. Так и сделал. Файловую систему для SSD. Люблю описывать сетевые протоколы, люблю писать интерпретаторы разных языков. Умею проектировать алгоритмы управления роботами. Умею расчитывать траектории движения без использования вещественных чисел с компенсацие ошибок округления. Умею писать реализации криптографических алгоритмов. Умею все это реализовать в аппаратуре программируеых матриц. Могу процессор спроектировать, правда не за один месяц.
  • Могу писать пркиладное программное обеспечение, имею богатый опыт разработки ПО для систем реального времени.
  • Могу эксперименты ставить по физике полупроводников, знаю толк в террагерцовой оптоэлектронике, правда лет 10 этим не занимался.
Почему-то не умею извлечь прибыль из этого, странно.
Скучаю, ищу работу...

Упражнения с BACnet - журнал операций

Я когда начинал описывать протокол BACnet, сделал такой сетевой "драйвер", который загружает формат Wireshark PCAP и подает на вход службы разбора пакетов APDU. Это позволило написать и отладить разбор форматов ASN кодирования BACnet.
С тех пор прошло много времени. Много кода... Сейчас Мне понадобилось сохранять конфигурацию контроллера в ответ на команду Re-инициализации. Вообще, как себе представляешь можно работать в контроллере с базой данных объектов? И при этом требуется сохранять конфигурацию между сессиями, при перезагрузке. Я придумал, что надо писать во флеш журнал операций модификации данных. Просто так, линейно вперед, пока место не закончится. Когда контроллер просыпается после перезагрузки он должен накатить журнал операций. Журнал, чтобы не плодить разные сущности я решил паковать прямо командами протокола, транзакциям. Например, транзакция понятная - создать объект, дописать свойства к объекту. По сути мы можем сохранять все операции записи типа WriteProperty, прямо из буфера обмена без изменения писать в журнал - во флеш. Вот такая замечательная идея.
Как отладить замечательную идею?
Журнал во флеш, флеш в контроллере. Контроллер перегружается... Не удобно.
Кроме контроллера есть сервер сбора данных, роутер BACnet. Предлагаю сохранить его конфигурацию в файл и поизучать каким-нибудь инструментом. Нужна отладка для журнала файловой системы.
Берем приложение, пригодное для отладки, для рассматривания журналов. Да вот оно - Wireshark. В итоге я пишу функцию, функционал, для записи данных не во флеш, а в файл со структурой PCAP и изучаю, что получилось в итоге. А в итоге имею последовательность команд протокола CreateObject, CreateObject, CreateObject... и на этом выводе примерно за день отлаживаю функцию экспорта объектов из базы данных в журнал.

База данных - это структура типа дерево объектов.
Foreach (DeviceInfo->ObjectTree, (Callback) save_object_cb, to_file);
Т.е. описываем некоторую процедуру обхода дерева с записью каждого элемента в файл/журнал. Так выглядит сериализация дерева.
Каждая запись состоит из шапки::{длина данных, тип записи}, пакета данных, контрольной_суммы. После чего можно записать следующую запись. В шапку можно записать и тип записи, в нашей концепции - это тип сервиса - CreateObject, DeleteObject, WritePropertyMultiple.
По сути нам не нужен этот идентификатор сервиса, потому что CreateObject без параметров можно расценивать как Delete, А операцию повторного создания объекта расценивать, как WritePropertyMultiple. Писать в журнал поток команд протокола подкупает своей простотой.
При ре-инициализации устройства через перезагрузку можно переписать журнал заново. При загрузке просто накатить этот журнал.

Есть проблема: как минимизировать объем данных в журнале, как делать инкрементную запись -- писать только те параметры, которые изменились с прошлой версии базы.
[07.01.2020]
Написал вывод в журнал, с учетом файлов. Получилось две операции: CreateObject, AtomicWriteFile.
Нашел свой старый проект, журналируемая файловая система для флеш. Тогда мы использовали четыре вида полей: APPEND - дописывание файлов, ATTRIB - дописывание/изменение атрибутов, CREATE - создание, DELETE - удаление, SKIP - технологическая операция - заполнение пробелов, для выравнивания на блок записи.
Сейчас думаю использовать старые идеи только расширить понятие файловой системы до объектов. CREATE -- создание объектов CreateObject, ATTRIB - запись атрибутов в формате WritePropertyMultiple, APPEND - дописывание файлов в форме AtomicWriteFile, DELETE - пометить на удаление из базы DeleteObject.

среда, 6 ноября 2019 г.

ARMv7m -- операции BFI и BFX и битовые строки

В приложении к криптографии я задался вопросом как "правильно" писать исходный код на языке Си, чтобы компилятор использовал инструкции BFI (Bit Field Insert) и BFX(Bit Field Extract).

Мой алгоритм берет 4 бита и использует их в качестве индекса массива. В примере я описываю "нелинейное биективное преобразование" по таблице подстановок (ГОСТ 34.12-2015 п.5.1.1).
Для отладки алгоритма я использую компиляцию в ассемблер и рассматриваю, какими командами компилятор выражается.
> gcc -march=native -O3 -o - -S magma.c
-- на экран выводится ассемблерный код оптимизированный под мой процессор
Или под целевую платформу ARMv7e-m
$ arm-none-eabi-gcc  -mthumb -mcpu=cortex-m4 -march=armv7e-m -mfloat-abi=hard -mfpu=fpv4-sp-d16 -o - -O3 magma.c -S

Исходно пишу на Си, рассматриваю результат.
Пишу определение, которое соответствует операции BFX

#define BEXTR(x, n, len) (((x) >> (n)) & ((1 << (len))-1))
В исходнике пишу так:

   for (i=0; i<8; i++){
       s = sbox[i][BEXTR(a,(i*4),4)];
       r |= (s & 0xF) <<(i*4);
   }
sbox - это таблица подстановок.
В ассемблерном коде возникает команда
ubfx...
Но мне никак не удается подобрать обратную операцию.
Пишу определение
#define BFI(x, y, n, len) x = ((x) & ~(((1 << (len))-1)<<(n))) | ((y & ((1 << (len))-1))<<(n))
Но в результате компилятор НЕ использует инструкцию bfi, определение НЕ работает.
Смотрю в документацию (Arm C Language Extensions, ACLE Q2 2019) нахожу пояснение, инструкция BFI описывается средствами языка Си, т.е. должно быть соответствующее ей выражение.
Через некоторое время, дошло, что в расширении языка С бывают свои битовые строки. Такое вот описание битовых полей (см код ниже) позволило однозначно задействовать инструкцию извлечения битовой строки и последующее помещение битовой строки обратно в регистр.

uint32_t T(uint32_t a)
{
    register union {
      struct {
        uint32_t u0:4;
        uint32_t u1:4;
        uint32_t u2:4;
        uint32_t u3:4;
        uint32_t u4:4;
        uint32_t u5:4;
        uint32_t u6:4;
        uint32_t u7:4;
      };
      uint32_t x;
    } r;
    r.x  = a;
    r.u0 = sbox[0][r.u0];
    r.u1 = sbox[1][r.u1];
    r.u2 = sbox[2][r.u2];
    r.u3 = sbox[3][r.u3];
    r.u4 = sbox[4][r.u4];
    r.u5 = sbox[5][r.u5];
    r.u6 = sbox[6][r.u6];
    r.u7 = sbox[7][r.u7];
    return r.x;
}
Это и есть биективное преобразование согласно ГОСТ, именно в таком виде оно попало в мою реализацию. -- Магия!!

В результате компилятор создает такой код:

   ubfx r2, r0, #16, #4   -- загрузить битовое поле
   add  r2, r2, r3
   ldrb r2, [r2,#64] @ zero_extendqisi2 -- загрузить байт из таблицы
   bfi  r0, r2, #16, #4   -- вставить битовое поле

CRISP 1.0 -- Протокол защищенного обмена для индустриальных систем


Протокол не важен. Важно что это первый нормальный набор векторов для Magma (ГОСТ 34.12-2015)в режиме генерации имитовставки (IMIT) и в режим гаммирования (CTR) для данных не выровненных. Сами режимы описаны в документе (ГОСТ 34.13-2015).

Протокол защищенного обмена выпущен ТС26 (Технический комитет по стандартизации российской криптографии) в виде методических рекомендаций.
Основная проблема в реализации режимов блочного шифрования - Это, не глядя в чужой код, получить точное соответствие с векторами. Сразу никогда не получается.
Где-то это мои проблемы - читать не умею. Где-то проблемы разработчика - описано плохо. Но чаще всего - это неточности в описании алгоритма, который надо как-то отладить, и ошибок в реализации может быть больше одной на пути рабочей версии.
Почему путь тернист. Берем сам алгоритм блочного шифра, "Магма" (ГОСТ 34.12-2015): описан, есть примеры - тестовые вектора. Взял алгоритм, сделал реализацию, но - не работает. Нигде в стандарте не сказано, в каком представлении даны числа в тестовых векторах. Порядок следования байт может быть Little-Endian или Network (Big-Endian). С ключами шифрования - тоже самое, они могут быть в нормальном порядке, задом наперед по 64 бита, задом наперед/совсем задом наперед - все число от старшего бита к младшему вывернуто. В процессе отладки я подбираю, как правильно записать, в какой последовательности, входные и выходные данные, чтобы числа сошлись. Сходятся.
Потом беру описание режима блочного шифра (ГОСТ 34.13-2015). В тестовых векторах та же путаница не ясно, как должны сходится эти вектора с изменение порядка следования байт или без изменения. Примеры почему то даны только для случая выровненных данных. В случае с усечением или без выравнивания на 64 бита, алгоритм остается не отлаженным.


 Теоретически мы можем предположить, что алгоритм работает так:
(подготовка ключа- надо решить как его правильно выкрутить, чтобы из сетевого представления получить числовое) -- (подготовка данных из последовательного представления в числовое) -- (подать на вход ключ и данные) -- (результат выкрутить в "сетевое" представление). При этом существуют некоторые неточности изложения, в результате которых такие же неоднозначные преобразования могут содержаться внутри самого алгоритма. Парадокс в том, что при этом вектора сходятся, потому что дают те же числа только в ином порядке.

Например, смотрим определения, чем отличается "бинарное представление", "байтовое представление" и числовое. Как разгадать этот ребус. По-русски ведь написано.
x||y -- конкатенация двоичных строк x и y ... при которой левая часть - слева, а правая - справа.
Потом пишут 1||0 - вот это чего за число. или это 0x02? Или это число 0x01? Оказывается это число 0x80.Оказывается оно таким после серии проб и ошибок.
binary() -- это представление символьной строки в виде байтовой строки. -- отпад.
byte() -- это представление числа в виде байтовой строки, при котором соответствующая итоговой байтовой строке двоичная строка .... -- вынос мозга.
вот у меня есть число 0xBAD как его представить в виде байтовой строки: "0B AD" или "AD 0B". Отсюда и берутся такие неоднозначности. Числа можно представить в "сетевом" порядке или в порядке свойственном аппаратуре LE|BE.
На самом деле, если четко договориться по терминам, что такое "двоичная строка", что такое "байтовая строка", что такое "число", что такое "символьная строка", как одно в другое переводится, то проблем нет. Авторы пытаются прикрыть неточность выражения мыслей обилием математических выражений типа Zn||...||Z2||Z1||Z0 - такие выражения никак не облегчаются понимание вопроса по выравниванию данных. Биты которые однозначно слева в "байтовом представлении" могут быть как слева так и справа. Но после прочтения не всегда это ясно.

Вот например каким количеством вариантов можно описать действие
IV = LSB_32(byte(SeqNum, 6)) -- по сути это означает что мы должны взять младшие 32 бита от байтового представления числа, при генерации байтового представления нужно взять 6 байт от числа, т.е 48 бит. Пишу варианты:
uint64_t SeqNum;
uint32_t IV = htonll(SeqNum); (1)
uint32_t IV = htonll(SeqNum<<16); (2)
uint32_t IV = (SeqNum); (3)
uint32_t IV = htonll(SeqNum)>>32; (3)
...
Или вот другая формула SN = 0^5||MSB_35(byte(SeqNum,6))
Если под функцией byte() автор понимает преставление байт из числового представления в сетевое, то это запишется так:
uint64_t SN = (swap64(SeqNum<<16)) & ((1ULL<<35)-1)); (3)
uint64_t SN = ((SeqNum>>(48-35)) & ((1ULL<<35)-1)); (4)
-- я реально перебираю все эти варианты, ориентируясь всего на один признак, результат равен или нет тестовому вектору. После подбора я нахожу, что автор имел ввиду:
SN = старшие_биты(0^5)||число(MSB_35(число(SeqNum,6)). -- сшивание битовых строк/чисел.
При укладке данных SN в "байтовую строку" происходит преобразование "в сетевом" порядке следования байт.
Строка = ...|| byte(SN,5) ||... -- число SN (40бит, 5 байт) разворачивается в порядке Big-Endian, в сетевом представлении. Откуда это следует? -- Я догадался.
Можно вчитываться в документацию, пытаться понять, все ли точно выражено в описании или все таки бинарное представление итоговой байтовой строки от числового представления сдвинутого на пять бит делается как-то иначе. А если в изложении несколько таких неточностей, или я не так понял и в моей реализации неточность или ошибка.
Может автор прав, есть однозначная запись битовых строк, операция MSB применима к битовым строкам, операция byte() преобразует число в битовую строку. Проверяем определения... Я понял, все дело в том, что битовые строки почему-то неявным образом переводятся в байтовые с выравниванием по старшему биту.
При использовании блочного шифра, я почему-то упускаю: надо данные из сетевого представления переводить в числовое перед использованием функции шифрования или надо результат функции шифрования переводить из локального в сетевое представление перед использованием. Такие вот проблемы. 

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

четверг, 31 января 2019 г.

Криптография на эллипсе

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

Самое сильное из моих впечатлений - это часики. Это понятная аналогия, проясняет сознание.
Возьмем к примеру окружность, часики с круглым циферблатом, уравнение окружности:
x^2+y^2 = 1

Примем за начало отсчета точку O = (0,1) . 12часов 00 минут.
К чему мы клоним, чтобы сразу стало понятно. Мы вводим группу вращения стрелок на циферблате. Время складывается: минуты в часы...
Стараемся думать про одну стрелку. за минуту стрелка отклонилась на (x1, y1). Эту точку можно выразить через синусы и косинусы единичного перемещения 360/60 = 6 градусов. Координаты единичного перемещения обозначим точкой G=(x1,y1)

Утверждаем что это у нас Группа точек на кривой, Группа в математическом смысле.
Свойство Группы:
1) Существование нейтрального элемента, такого что P+O = P, O = (0,1)
2) Существование обратного элемента для каждого члена группы, -P= (-x,y). Перевели стрелки назад. P+(-P)=O
Вводим операцию удвоения точки, с ней можно будет ввести операцию умножения на скаляр через удвоение и сложение.
3) 2P = ... осторожно можно споткнуться... = (cos(2ф), sin(2ф)) = (x^2-y^2, yx+xy)
4) Закон сложения точек (x1, y1)+(x2, y2) = (cos(a+b), sin(a+b)) = (x1x2-y1y2, y1x2+x1y2)

А теперь можно заставить часики ходить...
Алгоритм №1 умножения на скаляр Q = kP, k-раз по минуте.

Q:=O;
for i=.. downto 0 begin
  Q := 2Q;
  if (k_i !=0) Q:= Q+P;
end

Минуты считаются по модулю 60. Число 60 не является простым, его можно на множители раскладывать. Число Р назовем генератором группы, обозначим буквой G чтобы всех запутать.

Алгоритм №2 умножение лесенкой Монтгомери.

Q:=O; P=G;
for i=.. downto 0 begin
  if (k_i !=0){
      Q:= Q+P, P=2P;
  } else {
      P: =P+Q, Q=2Q;
  }
end
return Q
Эти алгоритмы не зависят от того как выглядит операция удвоения и сложения. Алгоритмов умножения можно придумать великое множество: справа налево, слева направо, комбинированные, с окном, со сложением и вычитанием, с разложениями и окнами.
Оба алогоритма можно свести к одной или двум операциям: удвоение точки Q=2Q и Q=2Q+G
Или иными словами мы на каждом шаге алгоритма вычисляем либо удвоение Q_{2n} зная Q_{n}, или Q_{2n+1} зная Q_{n}, Q_{n+1} и Q_1


Я знаю сколько было времени на часах, когда я их запустил - это мой секрет, могу выразить его в минутах d (число минут). Могу рассказать всем, что если умножить Q = dG получится некоторая точка с координатами (Q.x, Q,y) - которая однозначно связана с моим секретом - это будет точка для проверки подписи. Я хочу подписать сообщение. Мне нужно представить сообщение в виде числа m. Тогда подписанное сообщение - это показание часиков R = (m*d)G. Которое можно проверить с использованием открытой точки: R = mQ.

Цифровая подпись, ее неподдельная сущность, держится на том, что никто не может вычислить обратное число d, зная R, m и Q. Или плохо старается.

Все известные алгоритмы нахождения обратного числа держаться на Алгоритм № 3 НОД наибольший общий делитель. На базе него можно получить алгоритм деления или нахождения обратного числа по отношению к операции умножения. Для изготовления понадобится число типа скаляр и операция над точками - уполовинивание. Уполовинивание связано с неопределенностью при операциях с нечетными числами, которую надо как-то разрешать на каждом шаге алгоритма.
...

И тут пришел Монтгомери со своими кривыми алгоритмами и решил все "упросить": проекция x в операции удвоения не зависит от координаты y!
2P = (2x^2-1, 2xy)
Это значит, что мы можем считать удвоение без использования второй координаты. После этого берем паузу и думаем, а как теперь считать сложение точек без использования y- координаты.
x = x_2 x_3 - y_2 y_3 =... надо выразить через X координаты точек P Q и G.
x = 2 x_2 x_3 - x_1
Утверждение такое:
x_{2n} = 2x_{n}^2 -1
x_{2n+1} = 2x_{n} x_{n+1} - x_1
Начальное состояние для вычисления умножения при n=0 (x_{n}, x_{n+1}) = (1, x1).

По сути венец творения Монтгомери - это утверждение, что операцию вычисления x координаты при сложении точек на эллиптической кривой, можно представить в общем виде, как
x_{m+n} = f(x_m, x_n, x_{m-n}) вот и думай теперь над своим алгоритмом.


Откуда взяты идеи с часами и Алгоритмы Монтгомери:
https://eprint.iacr.org/2017/293.pdf -- оттуда