Показаны сообщения с ярлыком GFNI. Показать все сообщения
Показаны сообщения с ярлыком GFNI. Показать все сообщения

среда, 21 сентября 2022 г.

GF2p8 Кодирование Рида-Соломона в изоморфном представлении поля

Коды Рида-Соломона считаются через экспоненты и логарифмы по таблице. А у нас с внедренеием инструкций Intel GFNI появилось две другие возможности: использовать инструкцию умножения и аффинные преобразования.

Для примера я взял коды РС из проекта библиотеки генеации QR-кодов. Там используется полином поля G(2^8) P = 0x11D. Кодирование производится в цикле. При переходе к операциям умножения в поле 0x11D, цикл кодрования выглядит так.


	for(i = 0; i < data_length; i++) {
		uint8_t fb =  data[i] ^ Z[0];
		for(j = 0; j < ecc_length-1; j++) {
			Z[j] = Z[j+1] ^ gf2p8_mul(fb, g[j], 0x11D);
		}
		Z[j] = gf2p8_mul(fb, g[j], 0x11D);
	}
	for(i = 0; i < ecc_length; i++) {
		ecc[i] = Z[i];
	}

Преобразование алгоритма для расчета в изоморфном представлении поля с образующим полиномом 0x11B:


	const uint64_t M = 0xFFAACC88F0A0C080;
	const uint64_t Mr= 0xFFAACC88F0A0C080;

	for(i = 0; i < data_length; i++) {
		uint8_t fb =  affine(M,data[i]) ^ Z[0];
		for(j = 0; j < ecc_length-1; j++) {
			Z[j] = Z[j+1] ^ gf2p8_mul(fb, g[j], 0x11B);
		}
		Z[j] = gf2p8_mul(fb, g[j], 0x11B);
	}
	for(i = 0; i < ecc_length; i++) {
		ecc[i] = affine(Mr, Z[i]);
	}

Данное преобразование предполагает так же изоморфные преобразования компонент вектора генератора кодов (g). Векторизация алгоритма может быть выполнена по вложенному циклу.

[1]The comeback of Reed Solomon codes N.Drucker, Sh.Gueron, V.Krasnov

Приложение

Довеском идет код для расчета остатка от деления на 255. Когда в контроллере нет полиномиального умножения, приходится использовать способ расчета умножения и инверсии с использованием таблиц логарифмов и экспонент.
// Multiplication in Galua Fields
int gf_mul(int x, int y) {
    if (x == 0 || y == 0) return 0;
    return exp_[(log_[x] + log_[y])%255];
}
// Inversion in Galua Fields
int gf_inv(int x) {
    return exp_[255 - log_[x]];
}
// Division in Galua Fields
int gf_div(int x, int y) {
    if (x == 0) return 0;
    return exp_[(log_[x] + 255 - log_[y]) % 255];
}
Самая неудобная операция в этом случае - вычисление остатка от деления на 255. Избавился от деления.
// Остаток от деления на 255
uint8_t mod255(unsigned v) {
	uint32_t q = (v*0x1010102UL)>>(24);
	return q;
}
// Остаток от деления на 255
uint8_t mod255(unsigned v) {
	return v + (v>>8);
}

А этот вариант самый простой, применим для сложения логарифмов.

Такая вот оптимизация.

[2] Исходники на Github

вторник, 20 сентября 2022 г.

Крипто-Рубль. Изоморфные преобразования при вычислении Хеш функции Stribog

Для расчета хеш функции ГОСТ Р 34.11-2012 "Stribog" выполняется табличная подстановка в цикле - это общепринятая практика.

Если вернуться к определению алгоритма, его описывают последовательностью матричных операций S, P, L, X:
S - подстановка Sbox для каждого элемента 8x8 байт.
P - транспонирование матрицы 8х8 байт
L - линейное преобразование.
X - исключающее ИЛИ.

Линейное преобразование может быть представлено операцией умножения матрицы коэффициентов 8x8 на матрицу состояния 8x8 в поле GF(2^8) с полиномом 0x171.

LM = {
0x71, 0x05, 0x09, 0xB9, 0x61, 0xA2, 0x27, 0x0E,
0x04, 0x88, 0x5B, 0xB2, 0xE4, 0x36, 0x5F, 0x65,
0x5F, 0xCB, 0xAD, 0x0F, 0xBA, 0x2C, 0x04, 0xA5,
0xE5, 0x01, 0x54, 0xBA, 0x0F, 0x11, 0x2A, 0x76,
0xD4, 0x81, 0x1C, 0xFA, 0x39, 0x5E, 0x15, 0x24,
0x05, 0x71, 0x5E, 0x66, 0x17, 0x1C, 0xD0, 0x02,
0x2D, 0xF1, 0xE7, 0x28, 0x55, 0xA0, 0x4C, 0x9A,
0x0E, 0x02, 0xF6, 0x8A, 0x15, 0x9D, 0x39, 0x71,
Пояснение, как оно соотносится с с матрицей A. Колонка 1 матрицы M - это элемент A[0] = 0x8e20faa72ba0b470,
для которого выполнен разворот порядка следования бит.
8E=>71 20=>04 FA=>5F
Колонка 2 матрицы M - это элемент A[8],
. . .
Колонка 8 матрицы M - это элемент A[56],

Я независимо вывел эту таблицу, из представленнй таблицы подстановки. О возможности подобного представления читал в статье автора [1]. Продемонструю, как выполнить декомпозицию таблиц. Берем два числа: первоя строка, вторая и третья колонка:
A[1] = 0x47107ddd9b505a38
A[2] = 0xad08b0e0c3282d1c
Выделяем по маске 7 бит в каждом байте (A[1] & ~0x0101010101010101) >> 1)^A[2] видим числа 0x8E008E8E8E000000. В каждом байте стоит полином, по которуму производится редуцирование 0x171, с отраженным порядком бит будет 0x8E. Перед нами просто таблица умножения в поле GF(2^8)/[0x171] с отраженным порядком следования бит в байте.

Линейное преобразование может быть представлено таблично.
Cвойство L(a+b) = L(a) + L(b), как и свойство L(a*b) = L(a)*L(b) - действует и для таблицы.

Мы рассматриваем возможность реализации без таблицы, с использованием афинных преобазований и операци умножения в поле GF(2^8). Операция умножения может быть выполнена с использованием векторных инструкций умножения без переносов (Intel), полиномиального умножения (Arm Neon). Афинные преобразования могут выполняться с использованием векторных инструкций без обращения к таблицам.

Оптимизация под векторные инсрукции Arm Neon. 16 байт составляют вектор, используя вектора можно в два действия эмулировать работу операции умножения на константу. Для эмуляции используется свойство линейности преобразования.

Оптимизация под GF-NI:

Парадокс в том, что отечественные алгоритмы Stribog и Kuznyechik можно считать с использованием американской криптографии. Американская криптография использует умножение в поле GF(2^8) по полиному 0x11B. Отечетсвенная криптография использует умножение с редуцированием по полиному 0x171 в алгоритме Stribog, 0x1C3 в алгоритме Kuznyechik. Можно использовать афинные изоморфные преобразования, чтобы отобразить все матрицы и коэффициенты для использования умножения в поле 0x11b.

Матрицы изоморфного преобразования мы научились находить[2]. Для преобразования 0x171=> 0x11B можно использовать матрицы

матрицы изоморфного преобразования 171 => 11B "Stribog"
171=>11B 02 0E M =49A24EE22C984CF0 M-1 =FB44BAF0CC5A0A1C
171=>11B 02 17 M =C30E560C5240D828 M-1 =6D0A141C3A9C2046
171=>11B 02 52 M =B1A27CD0765CD234 M-1 =F5481A5CBE24D86E
171=>11B 02 54 M =29903A0892D47ABC M-1 =E5126608F2EC44F0
171=>11B 02 A0 M =C154448014AEDCC6 M-1 =1B8C164A06F81208
171=>11B 02 B4 M =1590FECC3E8E8CE6 M-1 =9960C6DA5E32485C
171=>11B 02 F6 M =090EFAA096722E1A M-1 =6FD8B46E36428C4A
171=>11B 02 FD M =A7541EDA2602426A M-1 =CB20B646D48660DA
Изоморфное представление матрицы LM в поле 0x11B:
unsigned char LM_11b = {
	0xF6,0x55,0x74,0xE4,0x56,0x3E,0xC1,0x2F,
	0x54,0xDF,0x17,0x9E,0xA9,0x60,0x43,0x02,
	0x43,0x1D,0x10,0x2E,0xEB,0xBB,0x54,0x65,
	0xA8,0x01,0x39,0xEB,0x2E,0xA1,0xE1,0xAD,
	0x93,0xAB,0x81,0x26,0x4E,0x42,0xF5,0xCE,
	0x55,0xF6,0x42,0x0D,0xFB,0x81,0xC7,0x0E,
	0xBA,0x5C,0xA6,0xEF,0x38,0x30,0xEC,0x71,
	0x2F,0x0E,0x07,0xD1,0xF5,0x2A,0x4E,0xF6,

Необходимо таже преобразовать таблицу Sbox в изоморфное представление Sbox_iso = MR·S·RM-1, где MR и RM-1 - матрицы прямого и обратного изоморфного преобразования из поля образованного полиномом 0x171 в поле 0x11b. R - отражение порядка следования бит.


void SPLX_11b(const uint8_t* s, uint8_t* r)
{
	int i,j,k;
	for (i=0; i<8; i++)
	for (j=0; j<8; j++) {
		uint8_t v= 0;
		for (k=0; k<8; k++) 
			v ^= gf2p8_mul(LM_11b[j*8+k], Sbox_iso[s[k*8+i]]);
		r[i*8+j] ^= v;
	}
}

Полученный алгоритм выполняет операции S,P,L,X. В цике выполняется операция умножения в поле GF(2^8)/0x11B. Оптимизация под ARM Neon может использовать полиномиальное умножение с последующим редуцированием.

Я озаглавил сообщение "Крипто-рубль..". Все алгоритмы добычи крипто-валют основаны на вычислении хеш-функций, этот алгоритм удалось выполнить исключительно в 8-битных числах. Это открывает возможность векторизовать алгоритм для параллельного вычисления до 64 хешей. Другая оптимизация - изоморфное преобразование в композитное поле позволяет получить эффективную реализацию алгоритма аппаратно в ASIC или на FPGA.


[1] Algebraic Aspects of the Russian Hash Standard GOST R 34.11-2012, O. Kazymyrov, V. Kazymyrova

среда, 30 июня 2021 г.

GF2p8 Умножение на аффинных преобразованиях

Пока изучал операции умножения в композитных полях, нашел два варианта аппаратной архитектуры умножителя в поле GF(2m).

Большая штука на рисунке -- умножитель в поле 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" признается готовым к использованию.

[0] E. Mastrovito, “VLSI architectures for computation in Galois fields,”Ph.D. dissertation, Linkoping University, Sweden, 1991
-- не читал, но все ссылаются, как на первоисточник идеи.
[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.

воскресенье, 27 июня 2021 г.

GF2p8 Код Грея через аффинные преобразования

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

Прямое преобразование выражается битовой операцией:
Gi = Bi⊕Bi+1.
Обратное преобразование можно выразить точно так же, рекуррентной формулой:
Bi = Gi⊕Bi+1.

Из за рекурсии расчет можно считать неэффективным, требуется цикл.

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

Преобразование в код Грея
M = 03060C183060C08016
10000000 =0x80
11000000 =0xC0
01100000 =0x60
00110000 =0x30
00011000 =0x18
00001100 =0x0C
00000110 =0x06
00000011 =0x03
Обратное преобразование из кода Грея
M-1 = FFFEFCF8F0E0C08016
10000000 =0x80
11000000 =0xC0
11100000 =0xE0
11110000 =0xF0
11111000 =0xF8
11111100 =0xFC
11111110 =0xFE
11111111 =0xFF

Путем перемножения матриц можно убедиться, что это обратные преобразования, произведение дает единичную диагональную матрицу

пятница, 25 июня 2021 г.

GF2p8 Аффинные преобразования в Композитные поля GF((2^4)^2)

Расширение и сужение поля

Идея в том, чтобы вместо операции над 8-и битными значениям использовать операции с 4-х битными значениями. Как осуществить переход?

В теории полей Галуа есть понятие расширение поля. Расширением GF(2^4) образованного полиномом q(x) до GF(2^8) c образующим полиномом p(z) является полином второго порядка GF((2^4)^2) s(y), но действующий с 4х битными коэффициентами. Запишем так:
s(y) = y2+t·y+n. – полином расширения поля, второго порядка,
где t – “trap” n – “norm” коэффициенты 4 бит. Коэффициенты хотелось бы выбрать таким образом, чтобы один из них (t) обратился в единицу.

Необходимо описать, как выглядят операции в поле: умножение и инверсия. Все числа в надстройке поля представляются парами: старшая и младшая часть:
a = a1·y+a0
y - символизирует сдвиг на 4 разряда.

Умножение в композитном поле выглядит следующим образом:
(a1·y+a0)(b1·y+b0) = (a1·b1)·y2 + (a0·b1+b0·a1)·y + a0·b0

Используя редуцирование по полиному y2+t·y+n =0, получаем:
(a1·y+a0)(b1·y+b0) = (a0·b1 + b0·a1 + a1·b1·t)·y + (a0·b0 + a1·b1·n) = r1·y + r0
Результат представлен в виде пары чисел разрядностью 4 бита.

Таким образом, можно выразить умножение, через умножение нижнем поле GF(2^4):
r1 = (a0·b1 + b0·a1 + a1·b1·t) mod q(x)
r0 = (a0·b0 + a1·b1·n) mod q(x)

Также можем выразить операцию инверсии, через инверсию в нижнем поле:
(a1·y + a0)-1 = (d·a1)·y + d·(a1·t + a0),
где d = (a12·n + a1·a0·t + a02)-1 –- инверсия в поле GF(2^4).
r1 = (d·a1) mod q(x)
r0 = d·(a1·t + a0) mod q(x)

Можно проверить справедливость формулы инверсии с использованием операции умножения.

Надо определить коэффициенты {t, n}, выбрать полином в поле GF(2^4).
Всего существует 4 значения n удовлетворяющих требованию t=1 и три нередуцируемх полинома четвертой степени. Для каждого набора параметров можно составить 8 матриц трансформации.

Нередуцируемые полиномы 4-ого порядка:
13: q(x) =x4+x+1
19: q(x) =x4+x3+1
1F: q(x) =x4+x3+x2+x+1

Вариант расчета: t=1, n=1100 – используется в литературе, применительно к полям GF(2^8)/0x11B (AES) и в сочетании с полиномом GF(2^4)/0x13.

Нас интересует составление матриц изоморфного преобразования между полем GF(2^8) и композитным полем GF((2^4)^2)

Для реализации нам понадобятся следующие алгоритмы:
1. умножение и инверсия в поле расширения;
2. умножение и инверсия в поле основания GF(2^4);
3. выбор генератора поля в композитном поле GF((2^4)^2);
4. составление и проверка матрицы трансформации между полем GF(2^8) и композитным полем GF((2^4)^2).

Операции в поле GF(2^4) можно реализовать методом табличной подстановки. Ниже привожу таблицы инверсии и редуцирования. Кроме того, можно реализовать таблицу квадратов и таблицу умножения на константу. Любую функцию одного переменного мы можем записать LUT таблицей (LookUp Table), четыре входа. При описаний логики для FPGA, четырех-битные LUT – оптимальный случай, поскольку, логика на четыре входа и таблица на четыре входа умещается в один логический элемент. При разработке алгоритмов с использованием векторных инструкций выигрыш в том, что вся таблица помещается в вектор 16 байт и для подстановки можно использовать функцию SHUFFLE.

Таблицы инверсии для полиномов 4-ого порядка:
13: 0,1,9,E,D,B,7,6,F,2,C,5,A,4,3,8
19: 0,1,C,8,6,F,4,E,3,D,B,A,2,9,7,5
1F: 0,1,F,A,8,6,5,9,4,7,3,E,D,C,B,2
Таблицы редуцирования старшей части полиномиального умножения:
13: 0,3,6,5,C,F,A,9,B,8,D,E,7,4,1,2
19: 0,9,B,2,F,6,4,D,7,E,C,5,8,1,3,A
1f: 0,F,1,E,2,D,3,C,4,B,5,A,6,9,7,8

Выбор генератора композитного поля можно совместить с генерацией таблицы инверсии, как было показано ранее. В данном случае нам нужно установить соответствие между (a) примитивным элементом поля GF(2^8) и (b) примитивным элементом поля GF((2^4)^2).

Алгоритм подбора генератора поля GF((2^4)^2):
1. Выбираем следующий примитивный элемент (b) поля GF((2^4)^2).
2. Составляем каким-то чудом матрицу прямого и обратного преобразования.
3. Выполняем проверку: произведение прямой и обратной матрицы дает единичную. 
Проверяем действительно ли полученное преобразование – изоморфизм. 
Если условие не выполняется, выбираем следующий примитивный элемент (b),
возвращаемся к шагу 1.
Свойства изоморфизма
    map(a×b) = map(a)×map(b)
    map(a⊕b) = map(a)⊕map(b)

Проверить можно гомоморфизм групповой по отношению к операциям сложения и умножения, для i,j принадлежащих множеству [0..2k-1] должно выполняться:
at = ai+aj => bt = bi+bj
ai+j = ai·aj => bi+j = bi·bj

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

Multiplicative Inverse GF((2^4)^2) 0x13, N=1100
 00 01 09 0E 0D 0B 07 06 0F 02 0C 05 0A 04 03 08
 AA A0 C7 CB 52 57 4F 4B FE F1 3E 3D 82 8A 2D 2F
 55 E9 50 E7 6A 8F 6C 87 B1 45 BA 41 2E 1E 2C 1F
 66 FA F5 60 96 48 4C 9F 5E D8 D5 5B 3F 1B 1A 3C
 BB 2B 7D C9 B0 29 7A C5 35 97 4E 17 36 9E 4A 16
 22 E6 14 56 E8 20 53 15 A9 75 D4 3B 72 A3 38 D9
 33 F4 E5 B7 EB BC 30 FB D3 9C 24 8E 26 86 DE 95
 99 B4 5C A2 A8 59 BF 90 D1 EF 46 C4 C8 42 E1 DC
 CC DA 1C 8B AF 9A 6D 27 C0 D7 1D 83 A5 93 6B 25
 77 BE A4 8D DF 6F 34 49 B5 70 85 AE 69 D2 4D 37
 11 AB 73 5D 92 8C B3 F2 74 58 10 A1 B8 FD 9B 84
 44 28 F3 A6 71 98 E4 63 AC FC 2A 40 65 EA 91 76
 88 D6 E2 F8 7B 47 CA 12 7C 43 C6 13 80 DB EC F7
 EE 78 9D 68 5A 3A C1 89 39 5F 81 CD 7F E0 6E 94
 DD 7E C2 F9 B6 62 51 23 54 21 BD 64 CE F6 D0 79
 FF 19 A7 B2 61 32 ED CF C3 E3 31 67 B9 AD 18 F0

Алгоритм (2) поиска матрицы прямого и обратного преобразования привожу ниже. Идея сводится к тому, чтобы из степеней генератора поля (a) составить диагональную матрицу и применив те же степени к генератору поля (b) сформировать по столбцам матрицу изоморфного преобразования. Для составления обратной матрицы поля меняются местами.

/* Алгоритм составления матрицы изоморфного преобразования
    a - примитивный элемент поля А
    P1 - образующий полином поля A
    b - примитивный элемент поля B 
    P2 - образующий полином поля B
 */
uint64_t map_ab(uint8_t a, uint32_t P1, uint8_t b, uint32_t P2)
{
	uint64_t m=0;
	int i;
	for(i=0;i<8;i++){
		int k=1;
		uint8_t v=1, r=1;
		while (v != (1<<(i))) {
			v = gf2p8_mul(v, a, P1);
			r = gf2p8_mul(r, b, P2);
		}
		m |= (uint64_t)r<<(8*(7-i));
	}
	return m_transpose(m);
}

В представленном алгоритме при выполнении операции умножения для композитного поля применяется арифметика композитного поля GF((2^4)^2), для поля AES арифметика поля GF(2^8), функции gf2p8_mul для композитного поля и для GF(2^8) по сути будут разные.

Матрицы изоморфного преобразования 11B => 13, N=1100
11B=>13 03 20 M =03A85C6870D2ACA0 M-1 =B1B042829AD45E54
11B=>13 03 22 M =737AF0C870D2ACA0 M-1 =A1B062A23A94BE14
11B=>13 03 35 M =91BC4AC4AE720CA0 M-1 =5FD0BAFAC21C2E9C
11B=>13 03 36 M =3FCE4664AE720CA0 M-1 =AFD01A5AE2DCCE5C
11B=>13 03 41 M =5DFC1C18DCAC72A0 M-1 =83700CACA4128692
11B=>13 03 45 M =81506EB8DCAC72A0 M-1 =B370CC6CE432E6B2
11B=>13 03 5B M =05E608CAA20CD2A0 M-1 =259024044C2A36AA
11B=>13 03 5E M =A7EADA6AA20CD2A0 M-1 =759064448C8A560A

Полученный результат позволяет выполнять операции композитного поля на инструкциях GFNI при изоморфном преобразовании из композитного поля в поле AES. И наоборот, позволяет перейти от операций в поле GF(2^8) к операциям меньшей разрядности GF(2^4). Такая техника может давать прирост производительности при реализации алгоритмов шифрования аппаратно и в микроконтроллерах.

В качестве тестового примера для проверки работоспособности методики выполняем генерацию S-Box AES в композитном поле с использованием изоморфного преобразования. При этом операция инверсии выполняется в арфиметике композитного поля.

S(x) = (A·M-1)·(M·x)-1 ⊕ C
A =F1E3C78F1F3E7CF816 - матрица аффинных преобразовний AES. 
C =6316
M =A7EADA6AA20CD2A016, -- матрица изоморфного преобразования
M-1 =759064448C8A560A16, -- обратное изоморфное преобразование
(A·M-1) =2F33DDCF49B6701E16

sbox[x] = affine(m_mul(A, M_), composite_inverse(affine(M, x)))^0x63;

Результат совпадает!

Расчет S-Box CLEFIA в композитных полях

Ранее мы рассматривали пример расчета нелинейного преобразования блочного шифра CLEFIA в изоморфных полях GF(2^8). Хотим закрепить навык выполнения расчетов в композитных полях и использовать технику преобразования в композитное поле, с расчетом инверсии на полиноме GF(2^4)/0x13.

Расчет нелинейного преобразования (биекции) дается следующим выражением:

S(x) = A2·(A1·x ⊕ C1)-1 ⊕ C2
A2,A1 - матрицы аффинных преобразовний.
инверсия считается в поле GF(2^8)/0x11D.
A1 =0x81605C6503015118 C1=0x1E
A2 =0x449002302058410A C2=0x69

Перейдем к расчетам в композитном поле:

S(x) = (A2·M-1)·((M·A1)·x ⊕ (M·C1))-1 ⊕ C2
инверсия считается в композитном поле GF((2^4)^2)/0x13.
M =B9C292428A441A2016
M-1 =9F547C4E5A805C0A16
(M·A1)  =FE297B311D0D060116,  (M·C1) =0x6A
(A2·M-1) =205054DA8048C31A16
// Расчет SBox - фрагмент кода
sbox[x] = affine(m_mul(A2,M_), composite_inverse(affine(m_mul(M,A1), x) 
                                               ^ affine(M,C1))) ^ C2;

Результат совпадает!


[1] D. Canright. A Very Compact S-box for AES
[2] A. Rudra et al. Efficient Rijndael Encryption Implementation with Composite Field Arithmetic

четверг, 24 июня 2021 г.

GF2p8 Аффинные преобразования (ч. 8): CLEFIA

Нашел еще один алгоритм блочного шифрования, на котором можно применить аффинные преобразования GFNI. В алгоритме CLEFIA блок нелинейной подстановки (биекции) Sbox1 формируется методом APA (Affine-Power-Affine). Преобразования (инверсия) выполняются на полиноме P(x)=x8+x4+x3+x2+1 = 0x11D.

S(x) = A2·(A1·x ⊕ C1)-1 ⊕ C2
A2,A1 - матрицы аффинных преобразовний. 
A1 =0x81605C6503015118 C1=0x1E
A2 =0x449002302058410A C2=0x69
CLEFIA S-Box:
 6C E9 C3 DA 0A 9D 4E 3D B8 38 B4 36 0C 34 13 D9
 BF 8F 94 74 E5 9C B7 DC 9E 4F 49 07 B0 2C 98 93
 12 B3 CD EB 41 E7 92 60 E3 3B 27 21 D2 19 E6 0E
 91 3F C7 11 A1 8E 2A BC 2B 0F C5 C8 87 F3 5B 8B
 FB 20 DE F5 84 A7 C6 CE D8 C9 51 65 43 EF A4 53
 25 31 9B 5D 0D 3E E8 D7 80 8A 69 FF 73 0B BA 5C
 6E 62 15 54 30 35 F6 52 A3 28 D3 16 AA FA 32 5E
 CF 78 ED EA 09 58 33 7B 63 46 C1 C0 A9 DF 1E 99
 55 86 C4 04 82 77 39 EC 40 97 90 18 83 DD 59 1F
 9A 24 06 37 A5 7C 64 56 48 D0 85 08 CA 26 61 6F
 7E 71 B6 6A 05 70 A0 D1 45 1C 23 8C 89 EE F0 AD
 7A 2F C2 4B 4D 5A DB 76 67 F4 2D 17 4A B1 CB A8
 B5 3A 47 22 4C 10 D5 72 CC E0 F9 00 FE E2 FD AE
 F8 F1 AB 5F 81 42 1B D6 BE A6 29 44 AF B9 57 F2
 D4 BB 66 75 50 9F 68 02 01 8D 7F 3C BD 88 1A AC
 F7 96 79 E4 6D FC A2 B2 6B 2E E1 03 95 14 7D 1D

Для аппаратного ускорения на иснтрукциях Intel GFNI применил изоморфные преобразования из поля 0x11D в 0x11B и обратно. Инверсия выполняется на полиноме P(x)=x8+x4+x3+x+1 = 0x11B с использованием инструкции gf2p8affineinv

S(x) = (A2·M-1)·((M·A1)·x ⊕ (M·C1))-1 ⊕ C2
A2,A1 - матрицы аффинных преобразовний. 
M, M-1 -- матрицы изоморфного преобразования
(A2·M-1) =0C70AA50A0B83F2216
00100010 =0x22
00111111 =0x3F
10111000 =0xB8
10100000 =0xA0
01010000 =0x50
10101010 =0xAA
01110000 =0x70
00001100 =0x0C
(M·A1) =931C707D4B19491816, (M·C1) = 0x18
00011000 =0x18
01001001 =0x49
00011001 =0x19
01001011 =0x4B
01111101 =0x7D
01110000 =0x70
00011100 =0x1C
10010011 =0x93

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

M = M-1 = FFAACC88F0A0C08016
10000000 =0x80
11000000 =0xC0
10100000 =0xA0
11110000 =0xF0
10001000 =0x88
11001100 =0xCC
10101010 =0xAA
11111111 =0xFF

Есть и другие варианты изоморфных преобразований 0x11D => 0x11B, которые можно применить в данном преобразовании. Результат совпадает!

P1   P2  g1 g2          M                   M-1
11D=>11B 02 03 M =FFAACC88F0A0C080 Mt =FFAACC88F0A0C080
11D=>11B 02 05 M =CFB00A10BC602840 Mt =DDE4F2E008A080AA
11D=>11B 02 11 M =5BD4D0B4F6487068 Mt =D9B2068A4AA0AAE4
11D=>11B 02 1A M =DDEE9CA64E38FC18 Mt =E12C36EE6EA0E4B2
11D=>11B 02 4C M =95D4EA8EEC0C2E2C Mt =7BC0D4F446A078E8
11D=>11B 02 5F M =6FAAD692CAC49EE4 Mt =2378BEF65CA0B22C
11D=>11B 02 E5 M =3BB06E74F85A567A Mt =ADE85E3EDAA02C78
11D=>11B 02 FB M =57EED8E22A228202 Mt =4F803A301CA0E8C0

понедельник, 21 июня 2021 г.

GF2p8 Аффинные преобразования (ч. 7)

Произведение матриц 8x8 бит

Инструкция аффинных преобразований gf2p8affine -- векторная. Аргументом является вектор элементов 64 бит, состоящий из матриц и вектор состоящий из октетов 8x8бит, над которыми выполняется преобразование. Мы хотим использовать эту инструкцию для перемножения матриц 8x8. Произведение матриц позволяет вычислять матрицы аффинного преобразования. В прошлом сообщении было показано, как выполнять транспонирование матриц. Операция аффинного преобразования обладает некоторой симметрией аргументов, которая позволяет менять аргументы местами при условии транспонирования.

gf2p8affine(A,B) = A·BT
gf2p8affine(A,BT) = A·B
gf2p8affine(AT,B) = (B·A)T= AT·BT

Чтобы инструкция gf2p8affine выдавала произведение матриц необходимо транспонировать второй аргумент.

// Произведение матриц A,B 8x8 бит
__m128i m_mul(__m128i A, __m128i B) {
	const __m128i E = _mm_set1_epi64x(0x0102040810204080);
	B = _mm_gf2p8affine_epi64_epi8(E, B, 0);// транспонирование
	A = _mm_gf2p8affine_epi64_epi8(A, B, 0);
	return A;
}

Результатом работы функции - произведение матриц A·B. Каждый аргумент состоит из двух матриц 8x8 по 64 бит, тип uint64x2_t. Параллельно можно выполнить два умножения матриц, независимо. Результат содержит две матрицы 8x8. Аналогично, при использовании инструкции 512 бит (_mm512_gf2p8affine_epi64_epi8), в результате получаем вектор uint64x8_t -- восемь матриц.

пятница, 18 июня 2021 г.

GF2p8 Аффинные преобразования (ч. 6)

Продолжаю тему, понял что некоторые разделы не раскрыты..

Умножение на константу в поле GF(2^8)

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

Матрица умножения состоит из столбцов, в которых располагаются степени числа 2. При умножении констант, должно быть выполнено редуцирование по образующему полиному.
M = [δ·27,δ·26,δ·25,δ·24,δ·23,δ·22,δ·21,δ·20]T

Матрицу нужно транспонировать, если заполнение производится по строкам.

// Алгоритм формирования матрицы умножения на константу в поле GF(2^8)
uint64_t m_mulC(uint8_t c, uint32_t Poly){
	uint64_t m = c;
	int i;
	for (i=1;i<8;i++){
		uint8_t v = gf2p8_mul(c, (1<<i), Poly);
		m = (m<<8) ^ v;
	}
	return m_transpose(m);
}

Редуцирование в поле GF(2^8)

Надо уметь выполнять редуцирование в форме аффинных преобразований - в одно действие. Дана старшая часть умножения без переноса, требуется выполнить редуцирование, т.е. перевести разряды в младщую часть. Редуцирование можно выполнить методом умножения на х8 mod Poly. В данном контексте матрица аффинного преобразования получается из выражения m_mulC(Poly & 0xFF, Poly), редуцирование - это матрица умножения в поле на специальную константу.

Таблица матриц редуцирования
0x11B: M=0xB1D3A6FD4B962C58
0x11D: M=0x71E2B51B478E1C38
0x165: M=0xDDBAA952A495F7EE
0x171: M=0x8D1A34685D37E3C6
0x177: M=0x4DD7E3C6C1CFD3A6
0x1C3: M=0x5BEDDAB468D0FBAD
0x1F5: M=0x2346AF5E9F1D1911

Транспонирование матрицы 8х8 бит в поле GF(2^8)

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

// Транспонирование матриц 8x8 бит
__m128i transpose_8x8 (__m128i m){
    __m128i E = _mm_set1_epi64x(0x0102040810204080ULL);
    return _mm_gf2p8affine_epi64_epi8(E, m, 0);
}
Транспонирование матрицы M = 0x0123456789ABCDEF;
    M             MT
11101111 =EF   11110000 =F0
11001101 =CD   11001100 =CC
10101011 =AB   10101010 =AA
10001001 =89   00000000 =00
01100111 =67   11110000 =F0
01000101 =45   11001100 =CC
00100011 =23   10101010 =AA
00000001 =01   11111111 =FF

Чтобы убедиться, что это действительно транспонирование, применяем повторно. Повторное транспонирование возвращает матрицу к исходному виду. Транспонирование в таком виде нужно, чтобы динамически составлять матрицы умножения на константу. Ниже привожу функцию динамической генератции матрицы аффинного преобразования, которая выполняет умножение на константу в поле GF(2^8)/0x11B

// Составление Матрицы умножения на константу в поле 0x11B
__m128i m_mulC_11B(uint8_t c) {
    __m128i E = _mm_set1_epi64x(0x0102040810204080ULL);
    __m128i m = _mm_gf2p8mul_epi8(E, _mm_set1_epi8(c));
    return _mm_gf2p8affine_epi64_epi8(E, m, 0);// транспонирование 8x8
}
Матрица редуцирования [1B, 0x11B] A= 0xB1D3A6FD4B962C58
Матрица умножения на  [A5, 0x11B] A= 0xABFCF9581A356AD5

GF2p8 Аффинные преобразования (ч. 5): блочный шифр ГОСТ "Кузнечик"

Блочный шифр "Кузнечик". Умножение в изоморфном представлении поля GF(2^8)

Весь цикл сообщений посвященный аффинным преобразованиям мы затеяли ради того чтобы представить этот и подобные ему алгоритмы. Ниже мы рассматриваем преобразование алгоритма шифрования ГОСТ Р 34.12-2015 Кузнечик к изоморфному представлению поля. Целью преобразования является ускорение работы за счет использования инструкций Intel GFNI.

Алгоритм зашифровывания выглядит, как последовательное примение частей преобразования: наложение ключа, биективное нелинейное преобразование S-Box, Линейные преобразования. Линейные преобразования можно свести к умножению в поле GF(2^8)/0x1C3. Можно записать линейные преобразование, как умножение вектора состояния (uint8x16_t, 16 байт) на матрицу коэффициентов (16х16 элементов). Алгоритм линейных преобразований (L) выглядит, как умножение матрицы на вектор:

// Линейное преобразование ГОСТ "Кузнечик"
uint8x16_t LMT[16];
// пересчет матрицы линейного преобразования
for(i=0; i<16; i++)
    LMT_iso[i] = _mm_gf2p8affine_epi64_epi8(LMT[i], M, 0);
// пересчет коэффициентов S-Box
for(i=0; i<256; i++)
    sbox_iso[i]= affine_byte(M, sbox[affine_byte(Mr, i)]);
/*! Линейное преобразование
В цикле выполняется умножение с редуцированием по полиному 0x11B */
uint8x16_t LS(uint8x16_t s) {
    uint8x16_t r = {0};
    int i;
    for(i=0; i<16; i++) {
        uint8_t v = sbox_iso[s[i]];
        r ^= _mm_gf2p8mul_epi8(_mm_set1_epi8(v), LMT_iso[i]);
    }
    return r;
}
/*! Зашифровывание */
uint8x16_t kuzn_encrypt(uint8x16_t* K, const uint8x16_t a)
{
    int i;
    a = _mm_gf2p8affine_epi64_epi8(a, M, 0);
    uint8x16_t s = a ^ K[0];
    for (i=0; i<9; i++)
        s = LS(s) ^ K[i+1];
    return _mm_gf2p8affine_epi64_epi8(s, Mr, 0);
}

Получается простое описание. Преобразование констант из матрицы коэффициентов LMT, как и аффинные преобразования S-Box, можно выполнить статически. Попробую алгоритм выразить формулой. Я убежден, что для программистов проще читать код, чем разбираться в формулах. Я применяю формулы там, где запись дает понимание. В формуле видны все места, куда надо вставить прямое и обратное преобразование. В частности, надо выполнить аффинные преобразования ключа K из представления поля 0x1C3 в 0x11B. Если дана функция двух переменных, то следует выполнить преобразования обеих переменных. Для разработчиков важно уяснить, что прямую и обратную матрицу добавляют в формулу в сочетании M-1·M·А, а потом, используя линейные преобразования, "вдвигают" матрицу внутрь формулы: M-1·А·M. Или же применяют к уравнению (слева) одно и тоже линейное преобразование: к левой и правой частям. В аналитических выкладках, надо четко понимать, какая операция в каком представлении поля и на каком полиноме выполняется.

gfmul[0x1C3](K, x) = M-1gfmul[0x11B](M·K, M·x);
sbox[0x1C3](x) = M-1sbox[0x11B](M·x);
sbox[0x11B](x) = M·sbox[0x1C3](M-1x);

Ключи (K) пересчитываются в представление поля 0x11B на этапе развертывания ключа. После подстановки формулы sbox[0x1C3] в качестве аргумента gfmul[0x1C3], произведение матриц M· M-1 = E дает единицу, сокращается:

gfmul[0x1C3](K, x) = M-1gfmul[0x11B](M·K, (M· M-1)sbox[0x11B](M·a));

Благодаря тому, что gfmul[0x1C3] выполняется в цикле, умножение вектора состояния на обратную и прямую матрицу можно вынести из цикла. В итоге в цикле аффинные преобразования не используются, матрицы изоморфного преобразования применяются только ко входным и выходным данным алгоритма зашифровывания, в то время как умножение производится на полиноме 0x11B.

В таблице представлены матрицы изоморфного аффинного преобразования из представления поля GF(2^8) 0x1C3 в 0x11B и обратно. В расчете можно использовать любую пару.

Таблица преобразований из представления поля 0x1C3 в 0x11B
M =0x5D0CE430CEE6BCD0 M-1 =0xC9248C8EB6BE7C4A
M =0x61D8CC543E9296A4 M-1 =0x359C60B0663A0EDA
M =0x2FA2EA44FA5AD66C M-1 =0x6332D8D6145ED06E
M =0x59A208A62EC29AF4 M-1 =0x61EC0A04B4F2D01C
M =0xED406082D258646E M-1 =0x89F844381A0602F0
M =0x5BD81880DC3CDA0A M-1 =0x494212C2C6360E08
M =0x0340F81A7C8C1EBA M-1 =0x87864834BAD4025C
M =0xC90C4A9E5604C632 M-1 =0x195A202216CC7C46

Идеальный S-Box

Глядя на алгоритм Кузнечик, видно его несовершенство.. Совершенный алгоритм -- тот, который сохраняя высокую степень нелинейности S-Box дает высокую производительность и высокую степень параллельности кода. Для этого нелинейное преобразование S-Box должно выражаться математическими операциями. В статьях предлагают S-Box полученные с использованием операции инверсии (возведения в степень) и аффинных преобразований, (англ. APA -- affine-power-affine).

s(x) = A2·Inv(A1·a(x)⊕C1) ⊕C2

В итоге формула "совершенного алгоритма" для следующего поколения ГОСТ выражается следующим циклом. При этом стандартизировать алгоритм можно с использованием любого из 30 нередуцируемых полиномов поля GF(2^8) и любых немыслимых аффинных преобразований, прежде всего в качестве основы для проектирования S-Box выбираются циркулянты - матрицы построенные на вращении полинома, которые дают обратимые преобразования. Возможно существует декомпозиция, которая представляет S-Box Кузнечик в совершенном виде.

// Предлагаемый алгоритм с использованием APA S-Box и матричного умножения.
uint8x16_t LS(uint8x16_t s) {
    s = _mm_gf2p8affine_epi64_epi8(s, A1, C1);
    s = _mm_gf2p8affineinv_epi64_epi8(s, A2, C2);
    uint8x16_t r = {0};
    for(i=0; i<N; i++) {
    	s = _mm_alignr_epi8(s,s,1);
        r^= _mm_gf2p8mul_epi8(s, LMR[i]);
    }
    return r;
}

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

среда, 16 июня 2021 г.

GF2p8 Аффинные преобразования (ч. 4)

Вычисление CRC на операции умножения в изоморфном представлении поля GF(2^8)

Мы приблизились к тому, ради чего все затевали.. Нет, вычисление CRC-8 не самоцель. Это самый простой пример применения арифметики Галуа, на котором можно проверить методику. Наша цель - изоморфные преобразования поля, чтобы свести все варианты расчета к одному полиному 0x11B.

Хотелось бы донести хотя бы до одного человека, как выводить алгоритмы. На первом шаге берем операцию, на которой работает CRC:
a(x)×x8 mod P(x) -- операция '×'-умножение без переноса или сдвиг.
Это можно представить иначе:
a(x)⊗(x8 mod P(x)) -- операция '⊗'-умножение в поле P(x).
Для полинома 0x11D выражение принимает вид: a(x)⊗0x1D -- операция '⊗'-умножение в поле P(x)=0x11D.

// CRC-8 в поле 0x11D
uint8_t crc = 0xFD;// начальное значение
uint8_t x8  = 0x1D;
for (i=0; i<length; i++) {
    crc = gf2p8_mul(crc ^ data[i], x8, 0x11D);
}
return crc;

M-1·M(a(x)×0x1D mod 0x11D) = M-1((M·a(x)) × (M·0x1D) mod 0x11B) -- операция '×'-сначала выполняется в поле P(x)=0x11D, а после преобразования в поле P(x)=0x11B.

// CRC-8 в изоморфном представлении 0x11D=>0x11B
// матрицы изоморфного преобразования {11D, 02} =>{11B, FB}
uint64_t M =0x57EED8E22A228202;
uint64_t Mr=0x4F803A301CA0E8C0;

uint8_t crc = affine_byte(M, 0xFD);// начальное значение
uint8_t x8  = affine_byte(M, 0x1D);// преобразование 0x11D=>0x11B
for (i=0; i<length; i++) {
    crc^= affine_byte(M, data[i]);
    crc = gf2p8_mul(crc, x8, 0x11B);
}
return affine_byte(Mr, crc);// обратное преобразование 0x11B=>0x11D
FD => 6A        57EED8E22A228202
1D => 11        57EED8E22A228202
DE => 7E        4F803A301CA0E8C0
CRC8I=7E ..ok

Операция вычисляется на двух инструкциях в цикле: gf2p8_mul и affine. Работает!!

GF2p8 Аффинные преобразования (ч. 3)

Эмуляция аффиных преобразований

Нужна эффективная замена для операций gf2p8affine. Каждый алгоpитм, который мы разрабатываем с использованием инструкций Intel GFNI должен исполняться и на платформе, где эти инструкции не представлены.

По определению любое аффинное преобразование обладает возможносью разлжения на табличные подстановки. Мы планируем эмулировать работу инструкции на платформе с поддержкой векторизации. Для работы с векторами мы можем использовать встроенную псевдо-функцию __builtin_shuffle(). Ранее мы использовали фунуцию affine_byte для расчета результата и подстановки аффинных преобразований от константы. Для данного преобразования необходимо составить две таблицы подстановок: для старших и младших разрядов.

// Аффинные преобразования с использованием таблиц 
uint8x16_t gf2p8_affine(uint8x16_t x, const uint64_t M)
{
    static const uint8x16_t lut_hi = {
        affine(M, 0x00),affine(M, 0x10), ...,affine(M, 0xF0)};
    static const uint8x16_t lut_lo = {
        affine(M, 0x00),affine(M, 0x01), ...,affine(M, 0x0F)};
    return shuffle(lut_hi, x>>4) ^ shuffle(lut_lo, x&0x0F)
}

Циркулянт и Анти-циркулянт

Циркулярная матрица - это матрица, полученная циклическим сдвигом некоторого ассоциированного полинома. Такие матрицы используются при составлении S-Box в качестве обратимого линейного преобразования. Мы говорим про матрицу специальной конструкции, которую можно описать, как произведение в поле GF(2^8) с редуцированием по полиному x8+1. Матрицу преобразования (A) можно определить двумя числами (полиномами), C(x) и D(х):

b(x)= C×a(x) ⊕ D mod (x8+1)

Пока что нам попались две такие циркуляторные (мне не нравится слово "циркулянт") матрицы: от AES S-Box и от SM4 S-Box. Мы можем составить все возможные циркуляторные матрицы и найти обратные к ним. Матрицу можно составить по сответствующему полиному степени n-1:
С(x)=с0 + с1·x + с2·x2+...+сn-1·xn-1 и константе D(x).


#define m_rotate(M, n) ((M)<<(8*(8-n)) ^ (M)>>(8*(n)))
uint64_t m_circulant(uint8_t C){
	return 
	  (C & 0x01 ? E: 0) 
	^ (C & 0x02 ? m_rotate(E,1): 0)
	^ (C & 0x04 ? m_rotate(E,2): 0)
	^ (C & 0x08 ? m_rotate(E,3): 0)
	^ (C & 0x10 ? m_rotate(E,4): 0)
	^ (C & 0x20 ? m_rotate(E,5): 0)
	^ (C & 0x40 ? m_rotate(E,6): 0)
	^ (C & 0x80 ? m_rotate(E,7): 0);
}

Операцию m_rotate можно выразить через циклические сдвиги единичной матрицы (E).

Сколько всего циркулянтов? как находить анти-циркулянты? Какой период циркуляции? Зачем нужны циркулянты? -- Всего 128 циркулянтов, включая единичную матрицу E. Все циркулянты имеют порядок вращения Mn=M, как и сам полином Cn mod 0x101 = C, период вращения равнен n=2, 4 или 8. Cn-1=1 -- есть такой признак. Обратной матрицей является полином полученный методом Cr = Cn-2 mod 0x101. Алгоритм нахождения обратного полинома: ищем порядок, находим полином на предыдущем шаге. Из этого описания получаем алгоритм поиска обратного циркулярного преобразования. Вернее, мы находим обратный полином (Cr), такой что Cr×C mod 0x101 =1, из которого составляется обратная циркуляторная матрица методом m_circulant(Cr).

// Алгоритм нахождения обратного циркулянта и порядка повтора
int m_circulant_order(uint8_t C, uint8_t* Cr) {
    uint8_t r= C, r_;
    int i;
    for (i=1; i<8; i++){
        r_= r;
        r = gf2p8_mul(r,C, 0x101);
        if (r==1) {
            if(Cr) *Cr = r_;
            break;
        }
    }
    return i;
}
// Вывод агоритма подбора циркулянтов (фрагмент)
С           M          Cr          M-1 
. . .
EA=> AE5DBA75EAD5AB57, AE=> EAD5AB57AE5DBA75, order=3
EC=> 6EDCB973E6CD9B37, 3B=> B973E6CD9B376EDC, order=7
EF=> EFDFBF7FFEFDFBF7, EF=> EFDFBF7FFEFDFBF7, order=1
F1=> 1F3E7CF8F1E3C78F, A4=> 4A942952A4499225, order=3
F2=> 9E3D7AF4E9D3A74F, 16=> D0A143860D1A3468, order=7
AES:
1F=> F1E3C78F1F3E7CF8, 4A=> A44992254A942952, order=3
SM4:
CB=> A74F9E3D7AF4E9D3, 85=> 43860D1A3468D0A1, order=7

вторник, 15 июня 2021 г.

GF2p8 Аффинные преобразования (продолжение)

В дополнение к предыдущим сообщениям по теме "Аффинные преобразования" и "изоморфные отображения между представлениями полей Галуа", продолжу разбор примеров применения ...

Вычисление CRC-8 с использованием аффинных преобразований

Операция вычисления CRC состоит из редуцирования в поле GF(2^8). В качестве примера используем CRC-8 с полиномом 0x11D. Для получения матрицы преобразования используем опыт полученный в предыдущем сообщении, составляем матрицу из операции сдвиг-редуцирование по полиному 0x11D, повторяем операцию 8 раз.

const uint64_t M = m_exp((E>>(8*1))^m_transpose(0x1D), 8);

M=71E2B51B478E1C3816 -- матрица преобразования. Результат от строки CRC8I("123456789")=0x7E.

// Вычисление CRC-8I с использованием аффинных преобразований
#define CRC8I_INIT 0xFD
#define CRC8I_CHECK 0x7E

for (i=0, crc=CRC8I_INIT; i<9; i++)
	crc = affine_byte(0x71E2B51B478E1C38, crc ^ data[i], 0);

Компилирую с использованием инструкций GFNI..


#include <intrin.h>
__m128i crc8_update(__m128i crc, uint8_t val) {
    __m128i M = _mm_set_epi64x (0x71E2B51B478E1C38, 0x71E2B51B478E1C38);
    __m128i v = _mm_maskz_set1_epi8(1, val);
    return _mm_gf2p8affine_epi64_epi8(v^crc, M, 0);
}

Результат совпадает!

Вычисление S-Box SM4 с использованием изоморфных преобразований

В прошлой части уже писал, как мне удалось преобразовать и сгенерировать S-Box с использовнием аффинных преобразований. Цель состоит в том, чтобы перейти от полинома 0x175 к вычислениям на полиноме 0x11B, что позволит использовать инструкции с аппаратной поддержкой инверсии и умножения в поле.


S[i] = affine_byte(A1, Inv[affine_byte((A1), i, 0xD3)], 0xD3);

// Выбираем изоморфное отображение 0x175=>0x11B
M_sm4  =0x51DE2604863CEA62;// прямая матрица
M_sm4_ =0x71BA0862FCB68CA2;// обратная
A2 = m_mul(A1, M_sm4_);// =0xD72D8E511E6C8B19
A3 = m_mul(M_sm4,  A1);// =0x34AC259E022DBC52
C3 = affine_byte(M_sm4, 0xD3);// =0x65

S[i] = affine_byte(A2, Inv[affine_byte(A3, i, C3)], 0xD3);

Количество преобразований от перехода в поле 11B не изменилось, всего две операции аффинного преобразования, которые могут быть представлены инструкцией gf2p8affineinv и gf2p8affine. В итоге S-Box можно расчитать за два такта для вектора 64 байт.

// Функция подстановки S-Box SM4 Intel GFNI
#include <intrin.h>
__m128i sm4_tau(__m128i state) {
    state = _mm_gf2p8affine_epi64_epi8(state, A3, 0x65);
    state = _mm_gf2p8affineinv_epi64_epi8(state, A2, 0xD3);
    return state;
}

Результат совпадает!

понедельник, 14 июня 2021 г.

GF2p8 Аффинные преобразования

Продолжаю тему затронутую в прошлом сообщении. Задача - научиться использовать на практике аффинные преобразования в поле GF2p8.

Пока что можно сказать, что аффинные преобразования предложенные набором инструкций Intel GF-NI реализуют всевозможные сдвиги и реализуют в себе битовые операции исключающее ИЛИ. Любые функции от одного переменного, которые можно представить, как сдвиги и битовое ИЛИ можно предствить этой операцией.

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

Матрицы - это что-то из области линейной алгебры, над матрицами определены операции: сложения и умножения. Алгеброй называем множество матриц 8x8 и операции определенные на этом множестве. Меня заносит в абстрактную алгебру...
Утверждаем существование нейтрального элемента -- единичной диагональной матрицы (E), такой, что умножение на единичную матрицу дает ту же матрицу E·M = M. Вообще-то не все равно, с какой стороны умножать матрицы... Утверждаем существование нейтрального элемента нулевой матрицы (Z) по отношению к операции сложения матриц, M+Z = M.

Утверждаем существование обратной матрицы для каждого элемента - обратная матрица имеет свойство M·M-1 = E - произведение прямой и обратной матрицы дает нейтральный элемент - единичную диагональную матрицу. Следует также определить обратную операцию и для сложения, M⊕(-M) = Z -- так определим. Обратная матрица по отношению к сложению - та же самая матрица, поскольку базовой операцией для сложения элементов матрицы явдяеься исключающее ИЛИ.

Операции сложения и умножения элементов матрицы определим на множестве состоящем из 0 и 1. Операция сложения - Исключающее ИЛИ, операция умножения - логическое И.

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

Матричное представление:
11111000 =F8
01111100 =7C
00111110 =3E
00011111 =1F
10001111 =8F
11000111 =C7
11100011 =E3
11110001 =F1

0xF1E3C78F1F3E7CF8 - сериализованное представление 64бит.

Нам понадобятся операции выделения строки(8бит) и столбца(8бит) из матрицы 8x8. Кроме того, понадобится функция, которая позволяет установить столбец и строку в матрице. Мы можем простыми сдвиговыми операциями поместить строку (8бит) в матрицу, в сериальном представлении 64бит, а для помещения строки в столбец нужно выполнить транспонирование матрицы после добавления строки.

// Выделение строки матрицы 8 бит
uint8_t m_row(uint64_t M, int n){
    return M>>(8*(7-n));
}
// Выделение колонки матрицы
uint8_t m_column(uint64_t M, int n){
    uint8_t r=0;
    uint8_t mask = 1<<n;
    int i;
    for (i=0;i<8;i++)
        if((M>>(8*(7-i))) & mask) 
            r |= (1<<i);
    return r;
}

Операцию умножения вектора на матрицу назовем аффинными преобразованиями.

// Аффинные преобразования векторов 8бит.
uint8_t affine_byte(uint64_t M, uint8_t byte, uint8_t inv) {
    uint8_t res = 0;
    int i;
    for (i=0; i < 8; i++) {
        // умножить строку на вектор, побитово:
        uint8_t x = m_row(M,i) & byte;
        // выставить бит в векторе результата 
        res |= parity(x)<<i;
    }
    return res ^ inv;
}

Функция parity - исключающее или от байта результата, в GCC я использовал встроенную псевдофункцию __builtin_parity().

При последовательном применении аффинных преобразований матрица воздействует на вектор 8бит. Это записывается, как произведение матриц - последовательное наложение. Матрицы можно перемножать и результат последовательного наложения матриц преобразования - наложение произведения матриц. Определим операцию умножения матриц.

// Умножение матриц
uint64_t m_mul(uint64_t m1, uint64_t m2)
{
    uint64_t r=0;
    int x,y;
    for (y=0; y<8; y++)
      for (x=0; x<8; x++){
        v  = m_row(m1, y) & m_column(m2, x);
        r |= (uint64_t)parity(v) << ((7-y)*8+x);
      }
    return r;
}

Для операции сложения матриц не нужна отдельная функция, потому что сложение матриц - это побитовая операция исключающее ИЛИ. Для порядка определим функцию, но дальше по тексту будем использовать оператор ^(шляпка).

// Сложение матриц
uint64_t m_add(uint64_t m1, uint64_t m2)
{
	return m1 ^ m2;
}

К матрицам нужно относится аккуратно, как к операторам, оператор воздействует на объект справа от оператора. Операторы нельзя менять местами -- последовательность применения операторов важна. M1·M2 не тоже самое, что M2·M1. Однако, сложение матриц обладает свойством M1⊕M2 = M2⊕M1.

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

// Экспонирование матриц
uint64_t m_exp(uint64_t m, int n)
{
    uint64_t r=E;// единичная диагональная матрица
    int i;
    for (i=0;i<n;i++)
        r = m_mul(r,m);
    return r;
}
// транспонирование матрицы 8x8
uint64_t m_transpose(uint64_t m)
{
	uint64_t r = 0;
	int x,y;
	for (y=0; y<8; y++){
		uint8_t row = m>>(8*(7-y));
		for (x=0; x<8; x++){
			if(row & (1<<x))
				r |= 1ULL<<((7-x)*8+y);
		}
	}
	return r;
}

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

Базовые операции

В этом разделе предложим некоторый базис операций и соответствующие им матрицы аффинных преобразований. Представленные ниже матрицы можно использовать для построения более сложных операций. Более сложные операции получаются методом последовательного наложения матриц (перемножение матриц) и сложения матриц. При этом сколько бы мы ни складывали матрицы и сколько бы не накладывали - результирующую операцию можно представить матрицей, одна процессорная инструкция с пропускной способностью две инструкции типа gf2p8affine за такт. Основной способ комбинирования базовых операций - наложть (перемножить) последовательно несколько одинаковых матриц или сложить матрицу одной операции с матрицой другой операции от общего аргумента.

010204081020408016 -- единичная матрица (E), результат без изменения

804020100804020116 -- отраженный порядок бит.

800102040810204016 -- циклический сдвиг влево на разряд

M= E>>(8*n) | E<<(8*(8-n)) -- циклический сдвиг влево на n разрядов.

020408102040800116 -- циклический сдвиг вправо на разряд

M= E<<(8*n) | E>>(8*(8-n)) -- циклический сдвиг вправо на n разрядов.

020408102040808016 -- арифметический сдвиг вправо на разряд с расширением знака

000102040810204016 -- сдвиг влево на разряд

889830E858B162C416 -- сдвиг влево с редуцированием по полиному 0x11B

M=(E>>(8*1)) ^ m_transpose(0x1B)

3849921D03070E1C16 -- сдвиг вправо -уполовинивание по полиному 0x11B

M=(E<<(8*1)) ^ m_transpose((0x11BULL>>1)<<56)

51D022F0946028C016 -- возведение в квадрат с редуцированием по полиному 0x11B.

Сдвиг с редуцированием - это базовая операция, на которой строится умножение в поле. Для каждого образующего полинома будет своя матрица.

Полиномиальное умножение на константу (умножение без переносов)

Можно предложить две операции умножения, первая операция возвращает - старшие 8 бит, вторая - младшие 8 бит.

Младшая часть умножения на константу P(8 бит)
M=(P&0x01? E>> 0: 0)
 ^(P&0x02? E>> 8: 0)
 ^(P&0x04? E>>16: 0)
 ^(P&0x08? E>>24: 0)
 ^(P&0x10? E>>32: 0)
 ^(P&0x20? E>>40: 0)
 ^(P&0x40? E>>48: 0)
 ^(P&0x80? E>>56: 0);
Старшая часть умножения на константу P(9 бит)
M=(P&0x100? E<< 0: 0)
 ^(P&0x080? E<< 8: 0)
 ^(P&0x040? E<<16: 0)
 ^(P&0x020? E<<24: 0)
 ^(P&0x010? E<<32: 0)
 ^(P&0x008? E<<40: 0)
 ^(P&0x004? E<<48: 0)
 ^(P&0x002? E<<56: 0);

Редуцирование по методу Барретта

Редуцирование предполагает использование двух умножений на константу (Ur, Pr), при этом используется только старшая часть умножения без переноса. Возможно в одной матрице совместить умножение на константу и последующее редуцирование. Или можно предложить операцию, которая на вход принимает старшую часть умножения и выдает результат редуцирования, который надо добавить к младшей части.

Умножение на константу с редуцированием по полиному можно выполнить в одну инструкцию.

B = Mpl·Mbh·Mch·A + Mcl·A = (Mpl·Mbh·Mch + Mcl)·A
A - байт (вектор из 8 бит) над которым производится операция.
Mcl - младшая часть умножения на константу
Mch - старшая часть умножения на константу
Mpl - младшая часть результата умножения без переноса на Pr=0x11B 
Mbh - старшая часть умножения на константу Барретта Ur=0x145
Mpl[0x11B]·Mbh[0x145] -- операция редуцирования по полиному 0x11B. 

Это и есть способ получения матрицы:

M = Mpl·Mbh·Mch + Mcl

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

Возведение в квадрат

Почему возможно? Причина прозаичная: (A⊕B)2 = A2 ⊕ B2 -- можно независимо считать квадрат каждого двоичного разряда. Для составления матрицы использую тот же принцип - по разрядом, каждый разряд возводится в квадрат независимо. Для старших разрядов выполняется редуцирование по полиному поля (0x11B).

M2=0x0100000000000000ULL// 1^2 = 0x01
  ^0x0000020000000000ULL// 2^2 = 0x04
  ^0x0000000004000000ULL// 4^2 = 0x10
  ^0x0000000000000800ULL// 8^2 = 0x40
  ^0x1010001010000000ULL//     = 0x1B
  ^0x0000202000202000ULL//     = 0x6C
  ^0x4040004000400040ULL//     = 0xAB
  ^0x0080008080000080ULL;//    = 0x9A
  
M2=0x51D022F0946028C0

Шестнадцатеричные константы справа получены методом сдиг-редуцирование. Сдвиги выполняются независимо по каждому разряду.

Изоморфные отображения между представлениями полей

Представлением поля назовем поля сформированные с использованием двух основных вещей -- это нередуцируемый полином и генератор поля. Геренатором поля может быть любой примитивный элемент поля. Все числа кроме нуля можно получить методом возведения в степень генератора поля. Множество чисел в поле Галуа начинается с нуля b и считаются последовательно только в человеческом представлении. В нечеловеческом представлении можно установить однозначное соответствие между степенями некоторого примитивного элемента одного поля и примитивного элемента другого поля. В частности рассматриваем изоморфизм между полем GF(2^8)/p=0x1C3 -- "ГОСТ Streebog/Кузнечик" и полем GF(2^8)/p=0x11B -- AES. Когда говорим про изоморфное преобразование надо уточнять, какие именно примитивные элементы поля использованы для отображения. Отображений великое множество, примерно половина элементов поля может оказаться примитивными элементами, вариантов отображения много.

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

Кроме того, нам понадобится функция, которая находит решение матричного уравнения вида A*x=C. Поиск матрицы отображения - это решение восьми таких уравнений, для каждого столбца матрицы отображения.

Свойства изоморфизма (отображения представлений поля) -- линейность.
	map(a×b) = map(a)×map(b)
	map(a⊕b) = map(a)⊕map(b)

Это означает, что можно выполнить переход от одного поля в другое, выполнить расчет и вернуться обратно. Это значит, что мы можем перейти из поля Национальной криптографии 1С3 (ГОСТ) в поле иностранной криптографии 11B, использовать аппаратную функицию умножения и поиска обратного числа, затем отобразить результат расчета в поле национальной криптографии. Хочу заметить, что свойство линейности отображения позволяет выполнить эту операцию в два действия с использованием таблиц подстановки и векторной инструкции SHUFFLE. Отображения бывают разные. Представим, что отобразить можно в поле, где операций меньше или сложность операции становится меньше. Например существует отображение в композитные поля типа GF((2^4)^2).

Разные утверждения в конечном поле можно доказывать методом перебора всех элементов поля. Можно решать уравнения методом подбора решения. Вот и матрицы мы будем искать методом подбора. Разыскиваем такое Отображение для двух выбранных полей. Для каждого поля мы выбираем примитивный элемент, проверяем его примитивные свойства, находим матрицу отображения и проверяем, чтобы отображение работало для всех степеней k=0..255. Если отображение не работает, выбираем другой примитивный элемент.

Определение изоморфного отображения представлений поля GF(2^8)
	a - примитивный элемент поля с полиномом A(x)
	b - примитивный элемент поля с полиномом B(x)
	bk = M·ak
	ak = M-1bk
	M·M-1 = E

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

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

/* Алгоритм генерации таблицы инверсии 
 Poly - образующий полином поля GF(2^8), например 0x1C3
 g - генератор поля, например 0x02, 0x03
 Функция возвращает TRUE, если расчет дал 255 элементов поля.
 */
int initialize_inv(uint8_t inv[256], uint32_t g, uint32_t Poly)
{
	int count=0;
	uint8_t p = 1, q = 1;
	uint8_t g_inv = gf2p8_inverse(g, Poly);
	do {
		p = gf2p8_mul(p, g, Poly);
		q = gf2p8_mul(q, g_inv, Poly);
		inv[p] = q;
		if (++count>255) {
			break;
		}
	} while(p != 1);
	inv[0] = 0;
	return (count==255);
}

Среди прочих отображений хочу выделить главное - отображение, которое приводит степени генератора к диагональному виду.
1)Надо найти такие k_i: ak_i = (1<<i), {01 02 04 08 10 20 40 80} в своем поле А.
2)Поместить в i-й столбец матрицы значение bk_i.
Преобразование M переводит число ak_i в число bk_i в поле B. Числа вида (1<<i) испольщзуют только i-й столбец матрицы отображения M. Поэтому Матрица отображения в столбцах содержит сами значения bk_i. Удобнее при этом соствлять матрицу по строкам, затем транспонировать к нормальному виду.

Привожу алгоритм составления матрицы отображения представлений поля.

/* Алгоритм составления матрицы отображения представлений поля
    PolyA - полином поля А,
    a   - генератор поля A,
    PolyB - полином поля B,
    b   - генератор поля B,
    Возвращает матрицу отображения
 */
uint64_t map_ab(uint8_t a, uint32_t PolyA, uint8_t b, uint32_t PplyB)
{
	uint64_t m=0;
	int i;
	for(i=0;i<8;i++){
		int k=1;
		uint8_t v=1, r=1;
		while (v != (1<<(i))) {
			v = gf2p8_mul(v, a, PolyA);
			r = gf2p8_mul(r, b, PplyB);
			k++;
		}
		m |= (uint64_t)r<<(8*(7-i));
	}
	// матрицу надо транспонировать
	return m_transpose(m);
}

Ниже привожу алгоритм подбора подходящего генератора поля для данного поля заданного полиномом и генератором {Poly1, g1} и искомого, заданного полиномом и начальным значением генератора поля {Poly2, g2}. В алгоритме вычисляется прямая (M) и обратная матрица (M-1). Проверка выполняется по критерию, является ли одна матрица обратной для другой, M· M-1 == E. Идея построения такого алгоритма мне досталась с формуа криптографов [JeffaReid], остальные алгоритмы можно считать авторскими. Отладку алгоритма и всей методики выполнял по матрице отображения мажду полиномами 0x175 SM4 и 0x11B AES представленной в патенте[US9800406, S.Guerоn]. В работах S.Guerоn алгоритм поиска отображения не рассматривается, хотя в одной из публикаций дана ссылка на препринт, сам препринт недоступен. Также нигде не говориться о существовании нескольких вариантов отображения между представлениями полей. То что отображений несколько показано ниже, любое из полученных отображений будет работать.

/* Алгоритм подбора подходящего генератора поля */
uint8_t map_ab_test(uint8_t g1, uint32_t Poly1, uint8_t g2, uint32_t Poly2)
{
	int count=0;
	do {
		g2++;
		while (!(initialize_inv(inv, g2, Poly2))) g2++;
		map = map_ab(g1, Poly1, g2, Poly2);
		map1= map_ab(g2, Poly2, g1, Poly1);
		if (m_mul(map1, map)== E) {
        // Отображение найдено!!
        	break;
		}
	} while (1);
	return g2;
}

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

матрица для 1C3 => 11B "Kuznyechik"
 A    B   a  b           M                    M-1
1C3=>11B 02 30 M =5D0CE430CEE6BCD0 M-1 =C9248C8EB6BE7C4A
1C3=>11B 02 70 M =61D8CC543E9296A4 M-1 =359C60B0663A0EDA
1C3=>11B 02 77 M =2FA2EA44FA5AD66C M-1 =6332D8D6145ED06E
1C3=>11B 02 7A M =59A208A62EC29AF4 M-1 =61EC0A04B4F2D01C
1C3=>11B 02 98 M =ED406082D258646E M-1 =89F844381A0602F0
1C3=>11B 02 C1 M =5BD81880DC3CDA0A M-1 =494212C2C6360E08
1C3=>11B 02 C9 M =0340F81A7C8C1EBA M-1 =87864834BAD4025C
1C3=>11B 02 DC M =C90C4A9E5604C632 M-1 =195A202216CC7C46
матрица для 171 => 11B "Stribog"
171=>11B 02 0E M =49A24EE22C984CF0 M-1 =FB44BAF0CC5A0A1C
171=>11B 02 17 M =C30E560C5240D828 M-1 =6D0A141C3A9C2046
171=>11B 02 52 M =B1A27CD0765CD234 M-1 =F5481A5CBE24D86E
171=>11B 02 54 M =29903A0892D47ABC M-1 =E5126608F2EC44F0
171=>11B 02 A0 M =C154448014AEDCC6 M-1 =1B8C164A06F81208
171=>11B 02 B4 M =1590FECC3E8E8CE6 M-1 =9960C6DA5E32485C
171=>11B 02 F6 M =090EFAA096722E1A M-1 =6FD8B46E36428C4A
171=>11B 02 FD M =A7541EDA2602426A M-1 =CB20B646D48660DA
матрица для 1F5 => 11B "SM4"
1F5=>11B 02 23 M =C772D8A4A8EEE4B0 M-1 =ED166E76E6BA48DC
1F5=>11B 02 3E M =3D0EE232C2D63888 M-1 =0B664A2E7A14D8AE
1F5=>11B 02 65 M =0BFCBE30BCDAE684 M-1 =8D14F098CEC61270
1F5=>11B 02 69 M =8BFC9C12C00A4A54 M-1 =51B45C94BC666070
1F5=>11B 02 86 M =01DEF6D40840181E M-1 =011ADA10501620A2
1F5=>11B 02 8E M =110E9E4EE0589406 M-1 =A9C64682A8B40AAE
1F5=>11B 02 CE M =6572562A78CC1692 M-1 =43B61CA4EA1A44DC
1F5=>11B 02 D6 M =51DE2604863CEA62 M-1 =71BA0862FCB68CA2
матрица для 165 => 11B "BelT"
165=>11B 02 18 M =DD04D04E96781490 M-1 =1B92021C42FA84C2
165=>11B 02 48 M =CD6480F638C48E2C M-1 =FBE27CDAEC265804
165=>11B 02 57 M =979A3EC0E25C52B4 M-1 =99FAD04A86E23C34
165=>11B 02 5B M =5F127C6AD66C3684 M-1 =6F260EF0244C928E
165=>11B 02 A5 M =2704C658680E9CE6 M-1 =CB4C026E9C84FAB0
165=>11B 02 E7 M =D3120A1C4892567A M-1 =6D580E5C5A924C22
165=>11B 02 EB M =7B9AA05E944A32A2 M-1 =E584D046F83CE238
165=>11B 02 F5 M =BF647E082E5A06B2 M-1 =F53C7C08325826D6
матрица для 11D => 11B "Kalyna"
11D=>11B 02 03 M =FFAACC88F0A0C080 M-1 =FFAACC88F0A0C080
11D=>11B 02 05 M =CFB00A10BC602840 M-1 =DDE4F2E008A080AA
11D=>11B 02 11 M =5BD4D0B4F6487068 M-1 =D9B2068A4AA0AAE4
11D=>11B 02 1A M =DDEE9CA64E38FC18 M-1 =E12C36EE6EA0E4B2
11D=>11B 02 4C M =95D4EA8EEC0C2E2C M-1 =7BC0D4F446A078E8
11D=>11B 02 5F M =6FAAD692CAC49EE4 M-1 =2378BEF65CA0B22C
11D=>11B 02 E5 M =3BB06E74F85A567A M-1 =ADE85E3EDAA02C78
11D=>11B 02 FB M =57EED8E22A228202 M-1 =4F803A301CA0E8C0

Приведенные выше преобразования относятся к отображениям между представлениями, изоморфные. А как насчет возможности отобразить представление поля на себя?

Отображение представления поля на себя
1C3=>1C3 02 02 M =0102040810204080 M-1 =0102040810204080
1C3=>1C3 02 04 M =D1B0C26084C0D830 M-1 =0B24924A2AAAA282
1C3=>1C3 02 10 M =BD745818F2E80C44 M-1 =65388CCCC4460CA6
1C3=>1C3 02 1F M =A34EE6EA26B84054 M-1 =A34EE6EA26B84054
1C3=>1C3 02 96 M =91CA5AF8B214D89E M-1 =075E3CE41C48A21A
1C3=>1C3 02 A0 M =65388CCCC4460CA6 M-1 =BD745818F2E80C44
1C3=>1C3 02 C3 M =075E3CE41C48A21A M-1 =91CA5AF8B214D89E
1C3=>1C3 02 FD M =0B24924A2AAAA282 M-1 =D1B0C26084C0D830

Первое отображение, понятно - на себя (Диагональная матрица E), второе - дает квадраты чисел.. Это все квадраты и корни квадратные в поле 0x1С3: 022=04, 042=10, 102=C3, C32=1F, 1F2=96, 962=A0, A02=FD, FD2=02. Всего существует 8 таких отображений, по этой же причине существует по 8шт отображений в каждое представление поля.

Чего не удается получить аффинными преобразованиями:
Инверсия числа - найти обратное не получается с использованием только аффинных преобразований.

Пример 1. Линейные преобразования AES S-Box

Матрица аффинных преобразорваний AES S-Box состоит из матриц циклических сдвигов

A = E ^ (E<<56 ^ E>>8) ^ (E<<48 ^ E>>16) ^ (E<<40 ^ E>>24) ^ (E<<32 ^ E>>32);
где E - единичная диагональная матрица.

A =F1E3C78F1F3E7CF816 -- результат сложения матриц.

Обратная матрица аффинных преобразорваний AES InvS-Box
A-1= (E<<56 ^ E>>8) ^ (E<<40 ^ E>>24) ^ (E<<16 ^ E>>48);

A-1=A44992254A94295216 -- обратная матрица

Прямая и обратная Матрица AES [NIST FIPS 197]
   |1 1 1 1 1 0 0 0| =F8            |0 1 0 1 0 0 1 0| =52
   |0 1 1 1 1 1 0 0| =7C            |0 0 1 0 1 0 0 1| =29
   |0 0 1 1 1 1 1 0| =3E            |1 0 0 1 0 1 0 0| =94
A =|0 0 0 1 1 1 1 1| =1F    A^{-1} =|0 1 0 0 1 0 1 0| =4A
   |1 0 0 0 1 1 1 1| =8F            |0 0 1 0 0 1 0 1| =25
   |1 1 0 0 0 1 1 1| =C7            |1 0 0 1 0 0 1 0| =92
   |1 1 1 0 0 0 1 1| =E3            |0 1 0 0 1 0 0 1| =49
   |1 1 1 1 0 0 0 1| =F1            |1 0 1 0 0 1 0 0| =A4
   
C =(1 1 0 0 0 1 1 0) =63    A^{-1}С=(0 0 0 0 0 1 0 1) =05

Проверим, что синтезированная матрица дает ожидаемый результат...


b = affine_byte(0xF1E3C78F1F3E7CF8, gf2p8_inverse(a), 0x63);

a = gf2p8_inverse(affine_byte(0xA44992254A942952, b, 0x05));

Пример 2. Линейные преобразования SM4

A1 =A74F9E3D7AF4E9D316 -- матрица аффинных преобразований для S-Box SM4.

   |1 1 0 1 0 0 1 1| =D3
   |1 1 1 0 1 0 0 1| =E9
   |1 1 1 1 0 1 0 0| =F4
A1=|0 1 1 1 1 0 1 0| =7A
   |0 0 1 1 1 1 0 1| =3D
   |1 0 0 1 1 1 1 0| =9E
   |0 1 0 0 1 1 1 1| =4F
   |1 0 1 0 0 1 1 1| =A7

C1=(1 1 0 1 0 0 1 1) =D3
// функция генерации таблицы S-Box SM4
C1 = 0xD3;
SBox_SM4[i] = affine_byte(A1, Inv_SM4[affine_byte(A1, i, C1)], C1);

В этом примере Inv_SM4[256] - таблица инверсии по полиному 0х175, получена функцией initialize_inv().

Выражение, которое должно было получиться
s(x) = A2·Inv(A1·x ⊕ C1) ⊕ C2.

В моем случае матрица A2 получилась равной A1 и С2 равно С1. Вероятнее всего где-то не совпал порядок бит, у меня или в работе S.Gueron, но мое выражение для SBox оказалось проще. В примере 3. рассмотрим вопрос оторажения порядка следования бит.

Пример 3. Отражение и вращение матрицы

Есть множество алгоритмов, где используется отраженный порядок бит. В матричном представлении нам достаточно поменять строки матрицы местами, чтобы подстроиться под заданный порядок бит на выходе операции. Или надо менять столбцы местами, чтобы подстроится под заданный порядок бит на входе. Иногда это сложно представить сразу. Мне кажется проще это делать методом последовательного наложения матриц. Представим себе матрицу, которая меняет положение бит в векторе на обратное, назовем матрицу R=8040201008040201. В общем случае замес может быть более сложный, поэтому представляем себе матрицу замеса R и обратную матрицу которая возвращает биты на место R-1. Есть несколько вариантов применения:

y = R·x -- меняет порядок следования бит.
y = R·A·x -- меняет порядок бит на выходе операции А.
y = A·R·x -- меняет порядок бит на входе операции А.
y = R-1·A·R·x -- меняет порядок на входе, на выходе возвращает к исходному.
   T=0x1122334455667788
    T             T*R            R*T           R*T*R
10001000 =88   00010001 =11   00010001 =11   10001000 =88
01110111 =77   11101110 =EE   00100010 =22   01000100 =44
01100110 =66   01100110 =66   00110011 =33   11001100 =CC
01010101 =55   10101010 =AA   01000100 =44   00100010 =22
01000100 =44   00100010 =22   01010101 =55   10101010 =AA
00110011 =33   11001100 =CC   01100110 =66   01100110 =66
00100010 =22   01000100 =44   01110111 =77   11101110 =EE
00010001 =11   10001000 =88   10001000 =88   00010001 =11

Для нашей матрицы R обратной операцией является тот же самый разворот R-1=R. В примере видно, что матрица отражается по вертикали или горизонтали, последнее преобразование разворачивает матрицу на 180градусов.

Пример 4. Выполнение операций в изоморфном представлении поля

y = Inv(x, Poly1);
y =M-1·Inv(M·x , Poly2);

Предлагается проделать операцию инверсии и убедиться, что результат идентичен для нескольких представлений поля. Матрица M выполняет изоморфноего преобразование из представления с полиномом 0x1C3 в представление с полиномом 0x11B, операция инверсии выполняется в поле 0x11B. Также можно проверить работу с отображением в то же представление поля 0x1C3.

Для x=0xA5: Inv(0xA5, 0x1C3)==0x5B
Poly1 Poly2 g1  g2 
 1C3 =>1C3, 02, 1F, M =A34EE6EA26B84054 M-1 =A34EE6EA26B84054
 1C3 =>1C3, 02, C3, M =075E3CE41C48A21A M-1 =91CA5AF8B214D89E
 1C3 =>11B, 02, 30, M =5D0CE430CEE6BCD0 M-1 =C9248C8EB6BE7C4A
 1C3 =>11B, 02, DC, M =C90C4A9E5604C632 M-1 =195A202216CC7C46

У меня работает! Результат идентичный: 0x5B.

пятница, 11 июня 2021 г.

GF2p8 Представления конечных полей

Речь пойдет про математику, и приложение новых инструкций Intel GFNI для криптографии и оптимизации разных алгоритмов.

Во-первых, никто не понимает, как пользоваться инструкциями GF-NI - нечеловеческая математика. Расширение состоит из трех инструкций, которые расчитывают матричные операции над полиномами: M*x+C, M*x-1+C и полиномиальное умножение с редуцированием A⊙B mod 0x11B. В систему команд введены инструкции для работы всего с одним полиномом (0x11B), в то время как этот полином используется мало в каких алгоритмах.

Первая операция M*x+C (аффинные преобразования gf2p8affine) не привязана никак к полиному(0x11B), можно сказать про полином вида x8+1, по которому выполняется редуцирование.

Операцию циклического сдвига A<<<n можно представить, как умножение и редуцирование по полиному x8+1. Запишется так: A×2n mod 0x101. Боже, откуда такое чувство, что меня никто не понимает?! Полиномиальное представление двоичных чисел - это просто термин такой, буквой x обозначается сдвиг на разряд. Вот что, например, написано про полиномиальное умножение в стандарте AES.. Rijndael S-box можно представить, как аффинное преобразование от обратного числа и записать в виде набора цикличекских сдвигов.
s = b ⊕ (b<<<1) ⊕ (b<<<2) ⊕ (b<<<3) ⊕ (b<<<4) ⊕ 6316

Это же утверждение записать через операцию полиномиального умножения (умнжения без переносов) и операцию редуцирования по полиному.

S = A×3110 mod 24710.

Это чтобы всех запутать сделано: числа даны в десятичном виде, операции в поле - умножение полиномиальное и редцирование по полиному x8+1.

Первое, что нашлось по теме использования инструкции - много циклических сдвигов можно представить, как аффинное преобразование [0x80.pl].

Из всего это следует, что подстановку AES S-box можно выполнить единственной инструкией gf2p8affineinv: S = A*inv(x)+C. Надо понять, как записать и протестировать эту возможность. Матрица преобразования записана в стандарте AES.

// Матрица аффинного преобразования, сдвиги: 0,1,2,3,4
|11111000|
|01111100|
|00111110|
|00011111|
|10001111|
|11000111|
|11100011|
|11110001|
// Матрица Обратного аффинного преобразования, сдвиги: 1,3,6
|01010010| 
|00101001|
|10010100|
|01001010|
|00100101|
|10010010|
|01001001|
|10100100|

Сдвиги в прямой матрице можно представить полиномом 7й степени (x7+x3+x2+x+1)=0xF1 сдвиги в обратной (x7+x5+x2)=0xA4. Для расчета обратного преобразования предлагается использовать утверждение (x7+x3+x2+x+1)^3 mod (x8+1)==(x7+x5+x2). Методом уножения на себя два раза из полинома 7й степени можно получить обратный ему полином.

Таким образом, утверждается, что имея полином A 7й степени, можно получить обратный ему полином B 7й степени, для которых выполняется условие A×B mod (x8+1) = 1. Не уверен, что это правило работает, но оставил.

// Нахождение обратного полинома 7й степени.
void inverse_af2(uint32_t poly)
{
	uint32_t mask = 0xFF;
	uint32_t p = poly & mask;
	uint32_t r;
	r = pmul(p, p);// умножение полиномиальное
	r = (r ^ (r>>8))& mask;// редуцирование
	r = pmul(r, p);// умножение полиномиальное
	r = (r ^ (r>>8))& mask;// редуцирование
	return r;
}

Навык поиска обратного значения нам понадобился, чтобы составлять прямые и обратные матрицы аффинных преобразований.

Обратный полином можно найти и методом подбора, перебрать все элементы поля, пока не найдется удовлетворяющий условию -- произведение на обратный элемент дает единицу.

// Операция умножения в поле GF(2^8) с полиномом Poly
uint32_t gf2p8_mul(uint32_t p, uint32_t i, uint32_t Poly){
    uint32_t r=0;
    while (i) {
        if (i&1)r^=p;
        if (p&0x80){
            p=(p<<1) ^ Poly;
        } else
            p=(p<<1);
        i>>=1;
    }
    return r;
}

Следующий инструмент разработки - составление таблиц преобразования, прямых и обратных. Тут надо вспомнить некоторое свойство. После g^(p-1) = g. Число g - генератор поля выбирается таким образом чтобы его степени проходили по всем элеементам - это, например, число 3. Методом уножения числа на 3 можно перебрать все элементы поля.

// Инициализация таблицы инверсии inv(x) == x-1
static void initialize_inv(uint8_t inv[256]) 
{
    uint8_t p = 1, q = 1;
    do {
        p = gf2p8_mul(p, 0x03, Poly);// умножить p на 3
        q = gf2p8_mul(q, 0xF6, Poly);// делить   q на 3
        inv[p] = q;
    } while p != 1);
    inv[0] = 0;
}

В цикле используется умножение на обратное число (0xF6), вместо деления на 3. Произведение чисел p и q единица, одно для другого на каждом цикле остается обратным числом. За счет выбора генератора, перебираются все числа, кроме нуля. В результате получаем таблицу обратных чисел.

Таблица обратных чисел для Poly=11B
 00 01 8D F6 CB 52 7B D1 E8 4F 29 C0 B0 E1 E5 C7
 74 B4 AA 4B 99 2B 60 5F 58 3F FD CC FF 40 EE B2
 3A 6E 5A F1 55 4D A8 C9 C1 0A 98 15 30 44 A2 C2
 2C 45 92 6C F3 39 66 42 F2 35 20 6F 77 BB 59 19
 1D FE 37 67 2D 31 F5 69 A7 64 AB 13 54 25 E9 09
 ED 5C 05 CA 4C 24 87 BF 18 3E 22 F0 51 EC 61 17
 16 5E AF D3 49 A6 36 43 F4 47 91 DF 33 93 21 3B
 79 B7 97 85 10 B5 BA 3C B6 70 D0 06 A1 FA 81 82
 83 7E 7F 80 96 73 BE 56 9B 9E 95 D9 F7 02 B9 A4
 DE 6A 32 6D D8 8A 84 72 2A 14 9F 88 F9 DC 89 9A
 FB 7C 2E C3 8F B8 65 48 26 C8 12 4A CE E7 D2 62
 0C E0 1F EF 11 75 78 71 A5 8E 76 3D BD BC 86 57
 0B 28 2F A3 DA D4 E4 0F A9 27 53 04 1B FC AC E6
 7A 07 AE 63 C5 DB E2 EA 94 8B C4 D5 9D F8 90 6B
 B1 0D D6 EB C6 0E CF AD 08 4E D7 E3 5D 50 1E B3
 5B 23 38 34 68 46 03 8C DD 9C 7D A0 CD 1A 41 1C

Заметим, что аналогичным образом можно составить таблицу экспоненты и логарифма числа, в качестве генератора используется число 2. Т.е. табличное представление умножения и инверсии через экспонирование и логарифирование чисел, с отображением на конечное поле 255 с обычной модульной арифметикой -- это нас не пугает, потому что традиционно используется в расчетах кодов Рида-Соломона. Интересно, можно ли аффинными преобразованиями получить экспоненту или логарифм числа, чтобы заменить таблицы на аффинные преобразования?

Вернемся к генерации таблиц.

S = A×inv(x)⊕C = AF(inv[x])⊕ 6316

Всего существует 30шт нередуцируемых полиномов восьмого порядка, которые можно использовать для построения полей GF(2^8). Утверждается, что все поля образованные нередуцируемыми полиномами 8-ого порядка изоморфны, т.е. имеют однозначное отображение одного множества чисел на другое множество. Это можно не доказывать, и так очевидное утверждение.

Другое утверждение -- представления полей GF(2^8) с разными полиномами можно отображать используя аффинные преобразования с матрицей 8x8 бит. Вот только как эти матрицы получать? Погуглить в интернете?! Нашел работы 2019-2021 года по разработке новых алгоритмов для пост-квантовой криптографии, нашел оптимизации таблиц подстановок типа S-box на основе таких преобразований. От Intel традиционно, внедрение новых инструкций Intel для криптографии продвигает Shay Gueron и группа товарищей, по его работам, патентам с его участием и статьям можно проследить развитие направления.

[US9800406] US Патент описывет преобразование представления GF(2^8) для вычисления S-Box китайского национального алгоритма блочного шифра SM4. Образующий поле AES полином x8 + x4 + x3 + x + 1=0x11B, в SM4 используется образующий полином x8 + x7 + x6 + x5 + x4 + x2 + 1=0x1F5. Преобразования Sbox записывается следующим образом.


SBoxSM4(x) = A2*inv(A1*x⊕C1)⊕C2
SBoxAES(x) = A*inv(x) + 6316
Матрицы SM4
   |1 0 1 0 0 1 1 1|    |1 1 0 0 1 0 1 1|
   |0 1 0 0 1 1 1 1|    |1 0 0 1 0 1 1 1|
A1=|1 0 0 1 1 1 1 0| A2=|0 0 1 0 1 1 1 1|
   |0 0 1 1 1 1 0 1|    |0 1 0 1 1 1 1 0|
   |0 1 1 1 1 0 1 0|    |1 0 1 1 1 1 0 0|
   |1 1 1 1 0 1 0 0|    |0 1 1 1 1 0 0 1|
   |1 1 1 0 1 0 0 1|    |1 1 1 1 0 0 1 0|
   |1 1 0 1 0 0 1 1|    |1 1 1 0 0 1 0 1|
   
C1=(1 1 0 0 1 0 1 1)TC2=(1 1 0 1 0 0 1 1)T
Матрица AES [NIST FIPS 197]
   |1 0 0 0 1 1 1 1|
   |1 1 0 0 0 1 1 1|
   |1 1 1 0 0 0 1 1|
A =|1 1 1 1 0 0 0 1|
   |1 1 1 1 1 0 0 0|
   |0 1 1 1 1 1 0 0|
   |0 0 1 1 1 1 1 0|
   |0 0 0 1 1 1 1 1|

C =(1 1 0 0 0 1 1 0)T

Патент представил матрицу преобразования элементов из поля SM4(0x1F5) в поле AES(0x11B)

Матрица отображения
  |0 1 1 0 0 0 1 0|
  |1 1 1 0 1 0 1 0|
  |0 0 1 1 1 1 0 0|
M=|1 0 0 0 0 1 1 0|
  |0 0 0 0 0 1 0 0|
  |0 0 1 0 0 1 1 0|
  |1 1 0 1 1 1 1 0|
  |0 1 0 1 0 0 0 1|

Отображение из SM4 в AES можно описать выражением

y = M*(A1*x⊕C1)

Обратное отображение из AES в SM4

s = A2*(M-1*A-1*u)⊕C2

Матрицы 8x8 можно перемножать между собой и сокращать выражения.

Матрицы отображения
     |01010010|             |11001011|
     |10111100|             |10011010|
     |00101101|             |00001010|
M*A1=|00000010|, A2*M-1*A-1=|10110100|
     |10011110|             |11000111|
     |00100101|             |10101100|
     |10101100|             |10000111|
     |00110100|             |01001110|

M*C1=0x65               C2 =0xD3

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

Для разработки необходимы методы: 1) Поиск матриц M для перехода между представлениями, 2) Декомпозиция S-Box 3) перемножение матриц, поиск обратных матриц, произведение матрицы на вектор - линейная алгебра. Из своего образования вспоминаю LU-разложение матриц, и свойства линейности матриц, которое позволяет представлять аффинные преобразования в виде перестановок, т.е снизить разрядность операции. Надо вспоминать линейную алгебру! Довольно много разных алгоритмов используют различные полиномы. Эти инструкции позволяют существенно ускуорить их работу, но пока нет опыта использования и надо затратить усилия на освоение математики и разработку алгоритмов синтеза таблиц.

Ниже привожу таблицу нередуцируемых полиномов..

Таблица нередуцируемых полиномов (30 шт):
x^8 + x^7 + x^5 + x^4 + 1
x^8 + x^6 + x^5 + x^4 + 1
x^8 + x^7 + x^5 + x^3 + 1
x^8 + x^6 + x^5 + x^3 + 1
x^8 + x^5 + x^4 + x^3 + 1
x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + 1
x^8 + x^6 + x^5 + x^2 + 1
x^8 + x^7 + x^6 + x^5 + x^4 + x^2 + 1
x^8 + x^7 + x^3 + x^2 + 1
x^8 + x^6 + x^3 + x^2 + 1
x^8 + x^5 + x^3 + x^2 + 1
x^8 + x^4 + x^3 + x^2 + 1
x^8 + x^7 + x^6 + x^4 + x^3 + x^2 + 1
x^8 + x^7 + x^5 + x^4 + x^3 + x^2 + 1
x^8 + x^7 + x^6 + x + 1
x^8 + x^7 + x^5 + x + 1
x^8 + x^6 + x^5 + x + 1
x^8 + x^7 + x^6 + x^5 + x^4 + x + 1
x^8 + x^7 + x^3 + x + 1
x^8 + x^5 + x^3 + x + 1
x^8 + x^4 + x^3 + x + 1 - AES
x^8 + x^6 + x^5 + x^4 + x^3 + x + 1
x^8 + x^7 + x^2 + x + 1
x^8 + x^7 + x^6 + x^5 + x^2 + x + 1
x^8 + x^7 + x^6 + x^4 + x^2 + x + 1
x^8 + x^6 + x^5 + x^4 + x^2 + x + 1
x^8 + x^7 + x^6 + x^3 + x^2 + x + 1
x^8 + x^7 + x^4 + x^3 + x^2 + x + 1
x^8 + x^6 + x^4 + x^3 + x^2 + x + 1
x^8 + x^5 + x^4 + x^3 + x^2 + x + 1

Следующая статья по теме - закрепление. Материала много для одного сообщения, но надо дополнить.

[S.Gueron et al. Speed up over the Rainbow] "All the representations of GF(2^8) are isomorphic. Therefore, it is possible to pass from one representation to another by means of multiplying by an 8 ×8-bit matrix"