[Из песочницы] Конечное поле GF(256) и немного магии14.03.2016 13:18

Введение


Будучи студентом я посещаю занятия по криптографии. И разумеется этот курс не мог обойти вниманием стандарт AES.

При реализации данного алгоритма встает вопрос о реализации полей GF (2^8), что будет освещено в данной статье. Будут рассмотрены: битовая магия для умножения элементов поля, шаблоны для генерации таблиц замен на этапе компиляции.

Вторая часть предполагает, что читатель имеет доступ к компилятору с поддержкой C++14. Первая часть будет написана в стиле Си.

Что собой представляют конечные поля
https://ru.wikipedia.org/wiki/Конечное поле

Поле GF (p)

Сперва рассмотрим как образуется поле с простым числом элементов GF (p).

Его элементы — числа \{0, 1, ..., p-1\}. Операции сложения и умножения — сложение и умножение по модулю p.

Например при p=7:

2 + 6 = 8%7 = 1
4×3 = 12%7 = 5

Поле GF (p^q)

На основе полей GF (p) строятся более общие поля GF (p^q), p — простое, q — натуральное.

Элементы таких полей — многочлены над полем GF (p):

\{0, 1, 2, \ldots , p-1, x, \ldots (p-1)x, ... , x^{q-1}, \ldots , (p-1)x^{q-1}\}

Сложение в данном поле будет непосредственно сложением данных многочленов.

Например при p=2, q=3:

(x + 1) + (x^2 + 1) = (x^2 + x)

С умножением чуть сложнее. Для того что бы определить его требуется многочлен Q (x) степени q, неприводимый над GF (p) (не существует двух многочленов меньшей степени, в произведении дающих данный). К счастью для любых p и q такой многочлен существует.
Когда многочлен Q (x) выбран, для того что бы найти произведение двух элементов a*b поля нужно:

  1. Найди произведение их многочленов a (x)*b (x)
  2. Найти остаток от деления этого произведения на многочлен Q (x). Это и есть a*b

Например при p=2, q=3:

Q (x) = x^3 + x + 1 (Неприводимый многочлен)

a = x^2 + 1, b = x^2
a (x) * b (x) = x^4 + x^2
(x^4 + x^2) % Q (x) = (x^4 + x^2 — x*Q (x)) % Q (x) = (x^3 + x) % Q (x) = (x^3 + x — Q (x)) % Q (x) = 1 = a * b

При p=3, q=2:

Q (x) = x^2 + x + 2 (Неприводимый многочлен)

a = x+2
b = 2x + 2
a (x) * b (x) = 2x^2 + (6%3) * x + (4%3) = 2x^2 + 1
(2x^2 + 1) % Q (x) = (2x^2 + 1 — 2Q (x)) % Q (x) = ((-2%3)x + (-3%3)) % Q (x) = x % Q (x) = x


Реализация


Итак, требуется реализовать следующие операции в поле GF (256) над многочленом x^8 + x^4 + x^3 + x + 1:

  1. Умножение
  2. Нахождение обратного


Начнем с произведение элементов поля:

Первое что приходит на ум — реализовать произведение многочленов наивным алгоритмом: 

uint16_t polynomeMul(uint8_t a, uint8_t b){  
    unsigned rez = 0;
    for(int i=0; i<8; ++i){ 
        rez ^= a * b&(1<


После чего написать функцию нахождения остатка от деления на многочлен.

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

А точно ли нельзя? Посмотрим что нам мешает реализовать произведение многочленов через простое произведение.

По формуле для произведения многочленов имеем:

(\sum_i a_i*x^i) * (\sum_i b_i*x^i) = \sum_i \sum_{n=0}^{i} a_n b_{i-n} x^i

Для произведения двух чисел в двоичной системе счисление получаем почти аналогиное выражение:

(\sum_i a_i*2^i) * (\sum_i b_i*2^i) = \sum_i \sum_{n=0}^{i} a_n b_{i-n} 2^i

Разница состоит в том, что для многочленов выражение \sum_{n=0}^{i} a_n b_{i-n} определяет коэффициент при x^i.

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

Как избавиться от переполнения? Очень просто. Рассмотрим следующую запись:

(\sum_i a_i*2^{4i}) * (\sum_i b_i*2^{4i}) = \sum_i \sum_{n=0}^{i} a_n b_{i-n} 2^{4i}

Поскольку при перемножении многочленов степени не выше 7 нельзя получить многочлен степени выше 14, разряду будет соответствовать сумма не более чем из 15 нулей и единиц (на самом деле легко убедиться что не более 8), а значит переполнение невозможно. Остается только преобразовать непосредственную сумму в сумму по модулю 2, выделив младший бит.

И так, если представлять полином как число, в котором каждому блоку соответствует блок из 4 бит, то произведение можно написать следующим образом: 

uint64_t polynomeMul(uint64_t a, uint64_t b){ 
    return (a*b) & 0x1111111111111111; // & 0b0001000100010001...
}


Теперь разберем произведение как элементов поля Галуа.

Посмотрим на полином x^8 + x^4 + x^3 + x + 1. В выбранном представлении он выглядит как q = 0×100011011. Невооруженным глазом видно большое количество нулевых блоков сразу после старшего блока. При умножении Q (x) на полином вида x^n(a_0 + a_1*x + a_2*x^2 + a_3*x^4) мы получим полином:

x^{n+8}(a_0 + a_1*x + a_2*x^2 + a_3*x^4) + \sum_i^{n+7}b_i x^i

или полином, старшие блоки которого есть a_0, a_1, a_2, a_3. Этим мы и воспользуемся для написания функции умножения: 

uint64_t galuaMul(uint64_t a, uint64_t b){ 
    uint64_t mul = polynomeMul(a, b);
    const uint64_t basePolynome = 0x100011011;
    mul ^= polynomeMul(basePolynome, mul>>48<<16); // Вычитаем произведение Q(x) и старших блоков имеющегося произведения (сдвинутых на 4 блока), оставляя ненулевыми не более 12 младших блоков
    mul ^= polynomeMul(basePolynome, mul>>32); // Таким же образом получаем многочлен не более 8 степени, сохраняя остаток от деления. После этой операции именно остаток будет записан в mul.
    return mul;
}


С умножением разобрались. Теперь нужно научиться находить обратный элемент по нему.

Вспоминаем, что поле без сложения и 0 образует группу из 255 элементов. Отсюда получаем, что для любого элемента x найдется число r, равное размеру подгруппы, образованной этим элементом, такое что x^r = 1. Так как порядок подгруппы — делитель порядка группы, x^{255} = x^{kr} = (x^r)^k = 1, что с свою очередь дает нам что x * x^{254} = 1. Тогда согласно определению обратного элемента имеем x^{-1} = x^{254}: 

uint64_t galuaPow(uint64_t a, unsigned n){
    // Быстрое возведение в степень.
    if(n==0){
        return 1;
    }else if(n%2 == 0){
        return galuaPow(galuaMul(a, a), n/2); // (a*a)^(n/2)
    }else{
        uint64_t square = galuaMul(a, a);
        return galuaMul(galuaPow(square, n/2), a); // a * (a*a)^[n/2] 
    }
}
uint64_t galuaInverse(uint64_t a){
    return galuaPow(a, 254);
}


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

Обозначим байт как 0bABCDEFGH. Первое что приходит в голову это умножение на 0b1001001 трех младших бит:

0bFGH * 0b1001001 = 0bFGHFGHFGH

0bFGHFGHFGH | 0b100010001 = 0bF000G000H, или три младших бита встали на свои места.

Аналогично проделывается над средней тройкой бит и старшей парой. Трюк был придуман. Но три умножения это как то слишком много. Нельзя ли делать тоже самое в раз по 4 бита? Рассмотрев многочисленные выборки бит я сумел найти лишь одну четверку с которой это работает:

image

Обратите внимание на младшие биты 7, 6, 1, 0 блоков. Для них характерны наличие нужного бита на своем месте и что не менее важно, невозможность переполнения за счет младших (относительно данных) бит.

Как было сказано, двух парных четверок бит я не нашел. Неудача? Не совсем. Если мы способны поставить семь из восьми бит на свои места использовав 2 умножения, мы можем поставить все 8, водрузив последний на свое место простым сдвигом.

uint64_t extendToGalua(uint8_t a){ 
    return (a & 0xC3) * 0x249249 & 0x11000011 |
           (a & 0x1C) * 0x1240 & 0x00011100|
           (a & 0x20) << 15;
}


С сжатием обратно дела обстоят проще. Следующее умножение показывает как сжать 4 бита:

image

С учетом этого код приобретает следующий вид: 

uint8_t shrinkFromGalua(uint64_t a){ 
    return (a & 0x11110000) * 0x249 >> 21 & 0xF0 |
           (a & 0x00001111) * 0x249 >> 9 & 0x0F;
}


Пусть считает компилятор


Зачем использовать довольно дорогие преобразования байт когда можно воспользоваться просчитанной заранее таблицей? В данной секции я поясню как просчитать её на этапе компиляции при помощи магии шаблонов на примере таблицы обратных элементов по сложению.

Первым делом всем ранее написанным функциям добавим спецификатор constexpr (с этого момента компиляция потребует поддержки с++14). Это позволит использовать данные функции в качестве аргументов шаблонов.

// Вспомогательная функция для укорачивания шаблонов
contexpr static uint8_t inverse(uint8_t x){
    return shrinkFromGalua(galuaInverse(extendToGalua(x)));
}

template
class GaluaTable{
public:
static constexpr auto& data = GaluaTable::data;
}



template
class GaluaTable<0, Data…>
public:
static constexpr uint8_t data[] = {Data…};
}



template
constexpr uint8_t GaluaTable<0, Data…>::data[];


Рассмотрим что происходит при попытке использовать GaluaTable: data.

Компилятор находит соответствующую специализацию шаблона, в которой data определена как GaluaTable<255, inverse(255)>:: data. Она в свою очередь определена через GaluaTable<254, inverse(254), inverse(255)>:: data и так далее.

На каждом шагу при этом мы имеем шаблон вида: GaluaTable. И так пока m не достигнет 0.

Когда m достигает 0 компилятору удается найти более конкретную специализацию шаблона (А компилятор всегда предпочитает более конкретную). Тут то и завершается рекурсивное задание классов и из последовательности в Data… создается сам массив, заимствованный всеми предыдущими классами.

Data… на это шаге будет ничем иным как inverse (0), inverse (1), …, inverse (255), что нам и было нужно.

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

Надеюсь, статья была чем-либо полезной.

© Habrahabr.ru