Вышли новые рекомендации к Российской национальной криптографии, режим шифрования MGM, который обещают использовать как основной и пожалуй единственный, если не считать рекомендации к протоколу CRISP. Режим MGM использует полиномиальное умножение в поле Галуа, по полиному x
128 + x
7 + x
2 + x
1 + 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
64 + x
4 + x
3 + x
1 + 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
128 + x
7 + x
2 + x
1 + 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
128 + x
7 + x
2 + 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
7 + x
2 + 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
64) и GF(2
128). Умножение состоит из двух частей - умножение без переноса и редуцирование. В цикле достаточно выполнять только умножение без переноса, а этап редуцирования можно вынести из цикла.
Несколько комментариев к статье
Операция 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, на новых процессорах этот вариант будет быстрее.
