Большая штука на рисунке -- умножитель в поле GF(2^4) построен на операциях типа И-ИсключающееИЛИ, каскадно, в два каскада. Из структуры умножителя видно, что в основе лежат аффинные преобразования, точнее такое разложение можно представить тремя аффинными преобразованиями, тремя инструкциями gf2p8affine. На первом каскаде - выполняется распределение и преобразования в две матрицы, выполняются две операции типа аффинного преобразования, на втором - выполняется операция аффинного преобразования и сложение аргументов.
На картинке представлена базовая идея построения умножителя [0,1]. Ничего особенно удивительного тут не вижу, потому что ранее представлял уже алгоритм формирования матрицы аффинного преобразования для операции умножения на константу. Тут будет что-то близкое: на первом каскаде формируется матрица умножения на константу А - матрица аффинного преобразования, второй каскад в чистом виде процессорная инструкция gf2p8affine.
Сразу можно отметить, что операция изоморфного преобразования ради умножения - три каскада, а тут мы видим два каскада - так быстрее должно получаться, потому что две операции аффинного преобразования выполняются параллельно.
Для решения предлагается составить матрицу Тёплица. Выполнить LU разложение. Перемножить матрицу L с вектором b, матрицу U с b. Произвести редуцирование с использованием матрицы Q:
d = L(a)·b e = U(a)·b c = L(a)·b + Q·U(a)·b = (L + Q·U)b
Нижняя и верхняя матрица заполняются битами аргумента (a) |a0 0 0 0 0 0 0 0| |a1 a0 0 0 0 0 0 0| |a2 a1 a0 0 0 0 0 0| L =|a3 a2 a1 a0 0 0 0 0| |a4 a3 a2 a1 a0 0 0 0| |a5 a4 a3 a2 a1 a0 0 0| |a6 a5 a4 a3 a2 a1 a0 0| |a7 a6 a5 a4 a3 a2 a1 a0| | 0 a7 a6 a5 a4 a3 a2 a1| | 0 0 a7 a1 a2 a3 a4 a5| | 0 0 0 a7 a1 a2 a3 a4| U =| 0 0 0 0 a7 a1 a2 a3| | 0 0 0 0 0 a7 a1 a2| | 0 0 0 0 0 0 a7 a1| | 0 0 0 0 0 0 0 a7| | 0 0 0 0 0 0 0 0|
Матрица Тёплица разбивается на две треугольные матрицы L U дополняющие друг друга: нижнюю (L) и верхнюю(U). Вместе они могут быть получены методом циклического сдвига строк. В системе команд Intel AVX512_VBMI есть подходящая инструкция для изготовления сдвиговой матрицы, vpmultishiftqb. Далее сдвиговую матрицу можно разделить на две треугольные с использованием маски.
Да, ладно! Вся эта логика, теоремы и доказательства построены вокруг одного тривиального утверждения: матрица умножения складывется из трех пребразований - это матрица умножения без переноса (младшая часть) плюс редуцированная по образующему полиному старшая часть умножения без переносов. Весь вопрос в том, как составить налету матрицу треугольно-диагональную (матрицу Тёплица) из двоичных разрядов числа A. Матрицы редуцирования мы уже вводили! См. предыдущие сообщения.
S = _mm_set1_epi64x(0x0807060504030201);
m = _mm_multishift_epi64_epi8( S, m);
8 7 6 5 4 3 2 1 9 8 7 6 5 4 3 2 10 9 8 7 6 5 4 3 11 10 9 8 7 6 5 4 12 11 10 9 8 7 6 5 13 12 11 10 9 8 7 6 14 13 12 11 10 9 8 7 15 14 13 12 11 10 9 8
// Алгоритм умножения c = a*b в поле GF(2^8)
// циклические сдвиги
S = _mm_set1_epi64x(0x0807060504030201);
// маска выделяет нижнюю матрицу L
ML = _mm_set1_epi64x(0x0103070F1F3F7FFF);
// матрица редуцирования для полинома 0x11B
Q = _mm_set1_epi64x(0xB1D3A6FD4B962C58);
v = _mm_set1_epi8(a); // размножить аргумент
v = _mm_multishift_epi64_epi8(S, v); // циклические сдвиги
// младшая часть умножения без переносов
lo = _mm_gf2p8affine_epi64_epi8(_mm_and_si128(ML, v), b, 0);
// старшая часть умножения без переносов
hi = _mm_gf2p8affine_epi64_epi8(_mm_andnot_si128(ML, v), b, 0);
// Q - матрица редуцирования по образующему полиному поля
hi = _mm_gf2p8affine_epi64_epi8(hi, Q, 0);
return hi ^ lo;
Код не отлаживал. Вероятно надо транспонировать или отразить некторые матрицы.
Почему остановился? Начал составлять алгоритм с другого конца. С этого мне сложно отладить. У нас есть готовые тесты на матрицу умножения на константу. Так вот, один из вариантов получения матрицы редуцирования Q - это вычислить матрицу умножения на специальную константу примерно так, отбросив старший бит:
Poly = 0x11B;
Q = m_mulC(Poly&0xFF, Poly);
Если отбросить старший бит, то это будет корень образующего полинома. Второй вывод: матрицы Тёплица я не хочу использовать, более подходящий термин - циркулянты. Мы вводили функцию генерации циркулянта - сдвиговой матрицы. Таким образом LU разложение дается через матрицу - циркулянт. Я по шагам показываю, как я синтезировал алгоритм..
const uint64_t mask = 0x0103070F1F3F7FFF;
C = m_circulant(a);
L = C & mask;
U = C & ~mask;
Надо заметить, что матрица L - это матрица малдщей части умножения без переносов на константу (a), а матрица U - матрица старшей части. В результате получаем референсный алгоритм умножения:
// Алгоритм умножения через аффинные преобразования
uint8_t gf2p8_mul_affine(uint8_t a, uint8_t b, uint32_t Poly)
{
const uint64_t mask = 0x0103070F1F3F7FFF;
uint64_t Q = m_mulC(Poly&0xFF, Poly);
uint64_t C = m_circulant(a);
uint64_t L = C & mask;
uint64_t U = C & ~mask;
return affine(L, b) ^ affine(Q, affine(U, b));
}
В качестве теста сравниваем вывод функции gf2p8_mul_affine с выводом функции умнжения gf2p8_mul, для любых a и b. Результат совпадает!
Матрицы редуцирования в контексте задачи - константы, для каждого полинома поля - своя.
Таблица матриц редуцирования (фрагмент) 0x11B: Q=0xB1D3A6FD4B962C58 0x11D: Q=0x71E2B51B478E1C38 0x165: Q=0xDDBAA952A495F7EE 0x171: Q=0x8D1A34685D37E3C6 0x177: Q=0x4DD7E3C6C1CFD3A6 0x1C3: Q=0x5BEDDAB468D0FBAD 0x1F5: Q=0x2346AF5E9F1D1911
При переходе к векторным инструкциям GFNI применим алгоритм составления циркулянта
const __m128i S = _mm_set1_epi64x(0x0807060504030201);
const __m128i E = _mm_set1_epi64x(0x0102040810204080);
__m128i m_circulant(uint8_t a) {
v = _mm_set1_epi8(a); // размножить аргумент
v = _mm_multishift_epi64_epi8(S, v); // циклические сдвиги
return _mm_gf2p8affine_epi64_epi8(E, v, 0);// транспонирование
}
// Алгоритм умножения c = a*b в поле GF(2^8)
// маска выделяет нижнюю матрицу L
ML = _mm_set1_epi64x(0x0103070F1F3F7FFF);
v = m_circulant(a);
// младшая часть умножения без переносов
lo = _mm_gf2p8affine_epi64_epi8(b, _mm_and_si128(ML, v), 0);
// старшая часть умножения без переносов
hi = _mm_gf2p8affine_epi64_epi8(b, _mm_andnot_si128(ML, v), 0);
// Q - матрица редуцирования по образующему полиному поля
hi = _mm_gf2p8affine_epi64_epi8(hi, Q, 0);
return hi ^ lo;
Следующий шаг - оптимизация. Мы могли бы объединить матрицы в одну используя алгоритм перемножения матриц.
__m128i m_mul(__m128i A, __m128i B) {
B = _mm_gf2p8affine_epi64_epi8(E, B, 0);
A = _mm_gf2p8affine_epi64_epi8(A, B, 0);
return A;
}
// Алгоритм синтеза Матрицы умножения на константу
__m128i m_mulC(uint8_t a, __m128i Q)
{
const __m128i ML = _mm_set1_epi64x(0x0103070F1F3F7FFF);
__m128i v = m_circulant(a);
return _mm_and_si128(ML, v) ^ m_mul(Q, _mm_andnot_si128(ML, v));
}
Количество аффинных преобразований можно сократить до двух, выполняемых параллельно.
// Матрица умножения на константу (a) в поле GF(2^8)
__m128i m_mulC(uint8_t a, __m128i Q)
{
const __m128i MLT = _mm_set1_epi64x(0xFFFEFCF8F0E0C080);
__m128i v = _mm_set1_epi8(a); // размножить аргумент
v = _mm_multishift_epi64_epi8(S, v); // циклические сдвиги
return _mm_gf2p8affine_epi64_epi8(E, _mm_and_si128(MLT, v), 0)
^ _mm_gf2p8affine_epi64_epi8(Q, _mm_andnot_si128(MLT, v), 0);
}
MLT - транспонированная маска для нижней матрицы L. Надо протестировать. В качестве теста, сверяю все три реализации функции m_mulC. Результат совпадает!
Цель достигнута, мы научились эффективно синтезировать матрицу умножения на константу в произвольном поле GF(2^8) без необходимости выполнения изоморфных преобразований. Алгоритм использует каскад из двух инструкций аффинного преобразования. Для каждого поля используется своя матрица редуцирования (Q).
Есть еще несколько тождественных преобразований с матрицами LU, которые могут ускорить алгоритм при условии, что матрица редуцирования Q-константа.
// вариант 2 завершения алгоритма m_mulC
return _mm_gf2p8affine_epi64_epi8(E^Q, _mm_and_si128(MLT, v), 0)
^ _mm_gf2p8affine_epi64_epi8(Q, v, 0);
// вариант 3 завершения алгоритма m_mulC
return _mm_gf2p8affine_epi64_epi8(E, v, 0)
^ _mm_gf2p8affine_epi64_epi8(Q^E, _mm_andnot_si128(MLT, v), 0);
Эти преобразования можно выразить математическими формулами:
M = L(a)+Q·U(a) = (L(a) + U(a)) + (Q−E)·U(a) = C(a) + (Q−E)·U(a)
M = L(a)+Q·U(a) = (L(a) − Q·L(a)) + Q·(L(a)+U(a)) = (E−Q)·L(a) + (Q)·C(a)
Сформулировал вариант алгоритма, при котором возможно параллельное исполнение инструкций аффинного преобразования. Чтобы избавиться от масок, ввожу две параллельные инструкции циклического сдвига "multishift", матрицы LU расчитываются параллельно.
// Матрица умножения на константу (a) в поле GF(2^8)
__m128i m_mulC_bp(uint8_t a, __m128i Q)
{
const __m128i SU = _mm_set1_epi64x(0x0807060504030201);
const __m128i SL = _mm_set1_epi64x(0x00FFFEFDFCFBFAF9);
__m128i v = _mm_set1_epi64x(a);
__m128i L = _mm_multishift_epi64_epi8(SL, v); // циклические сдвиги
__m128i U = _mm_multishift_epi64_epi8(SU, v); // циклические сдвиги
return _mm_gf2p8affine_epi64_epi8(E, L, 0)
^ _mm_gf2p8affine_epi64_epi8(Q, U, 0);
}
Доказательсвто оптимальности кода. Хочу показать, как выглядит оптимизация кода. Для этого вывожу ассемблерный код и симулирую загрузки ядра процессора полученного кода с использованием анализатора производительности машинных кодов (llvm-mca).
$ gcc -O3 -S -o mul_c.s mul_c.c $ llvm-mca.exe --mcpu=icelake-server -timeline mul_c.s
/* Ассемблерный код */
vpmultishiftqb %xmm2, %xmm4, %xmm3
vpmultishiftqb %xmm2, %xmm5, %xmm2
vgf2p8affineqb $0, %xmm3, %xmm0, %xmm0
vgf2p8affineqb $0, %xmm2, %xmm1, %xmm1
Временная диаграмма загрузки ядра процессора [0,0] DeER . . vpmultishiftqb %xmm2, %xmm4, %xmm3 [0,1] DeER . . vpmultishiftqb %xmm2, %xmm5, %xmm2 [0,2] D=eER. . vgf2p8affineqb $0, %xmm3, %xmm0, %xmm0 [0,3] D=eER. . vgf2p8affineqb $0, %xmm2, %xmm1, %xmm1 [1,0] D==eER . vpmultishiftqb %xmm2, %xmm4, %xmm3 [1,1] D==eER . vpmultishiftqb %xmm2, %xmm5, %xmm2 [1,2] .D==eER . vgf2p8affineqb $0, %xmm3, %xmm0, %xmm0 [1,3] .D==eER . vgf2p8affineqb $0, %xmm2, %xmm1, %xmm1
Во временной диаграмме показана загрузка и стадии обработки инструкций (D-декодирование, e-исполнение, R-запись в регистр). Из диаграммы видно, что инструкции могут выполняться параллельно. Производительность достигает два такта на цикл алгоритма, второй цикл выполняется за два такта. Алгоритм m_mulC_bp "bit parallel" признается готовым к использованию.
-- не читал, но все ссылаются, как на первоисточник идеи.
[1] Low complexity bit parallel architectures for polynomial basis multiplication over GF(2m).
A. Reyhani-Masoleh, M. Hasan: Mathematics, Computer Science IEEE Transactions on Computers, 2004.
