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

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

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.