Конечное поле GF(256) и немного магии

Введение


Будучи студентом я посещаю занятия по криптографии. И разумеется этот курс не мог обойти вниманием стандарт 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<<i) // Прибавляем по модулю 2 многочлен b_i * (x^i) * a(x)
    }
    return rez;
}

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

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

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

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

(\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 = 0x100011011. Невооруженным глазом видно большое количество нулевых блоков сразу после старшего блока. При умножении 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 galoisMul(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 galoisPow(uint64_t a, unsigned n){
    // Быстрое возведение в степень.
    if(n==0){
        return 1;
    }else if(n%2 == 0){
        return galoisPow(galoisMul(a, a), n/2); // (a*a)^(n/2)
    }else{
        uint64_t square = galoisMul(a, a);
        return galoisMul(galoisPow(square, n/2), a); // a * (a*a)^[n/2] 
    }
}
uint64_t galoisInverse(uint64_t a){
    return galoisPow(a, 254);
}

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

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

0bFGH * 0b1001001 = 0bFGHFGHFGH

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

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

image

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

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

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

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

image

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

uint8_t shrinkFromGalois(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 shrinkFromGalois(GaloisInverse(extendToGalois(x)));
}
<p>template<int N, int… Data>
class GaloisTable{
public:
static constexpr auto& data = GaloisTable<N-1, inverse(N-1), Data…>::data;
}</p>

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

<p>template<int… Data>
constexpr uint8_t GaloisTable<0, Data…>::data[];</p>

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

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

На каждом шагу при этом мы имеем шаблон вида: GaloisTable<m, inverse(m), inverse(m+1), …, inverse(255)>. И так пока m не достигнет 0.

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

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

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

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

Update: Префикс Galua заменен на Galois.
AdBlock has stolen the banner, but banners are not teeth — they will be back

More
Ads

Comments 10

    0
    Крутой вариант шаблонной магии С++! :)
      0
      Функцию galuaMul можно упростить:

      uint64_t galuaMul(uint64_t a, uint64_t b){ 
          uint64_t mul = polynomeMul(a, b);
          const uint64_t basePolynome = 0x100011011;
          return mul ^ polynomeMul(basePolynome, mul>>32);
      }

      Не важно сколько нулей и единиц в basePolynome, при делении многочленов неполное частное зависит только от старшего одночлена. В формате программы, неполное частное от деления mul на basePolynome — это mul>>32.
        +1
        Оптимизировать — так до конца. Число операций galuaMul при нахождении обратного элемента можно уменьшить с 13 до 11 если вместо общего алгоритма воспользоваться разложением 254 = 2 * (1 + 2 * 3 * 3 * (1 + 2 * 3)):

        uint64_t galuaInverse(uint64_t a){
            uint64_t a2 = galuaMul(a, a);
            uint64_t a6 = galuaMul(a2, galuaMul(a2, a2));
            uint64_t a7 = galuaMul(a6, a);
        
            uint64_t a21 = galuaMul(a7, galuaMul(a7, a7));
            uint64_t a63 = galuaMul(a21, galuaMul(a21, a21));
            uint64_t a126 = galuaMul(a63, a63);
            uint64_t a127 = galuaMul(a126, a);
        
            return galuaMul(a126, a126);
        }
          +2
          Скажите, чем вас не устроила таблица логарифмов?
            0
            Присоединяюсь к вопросу. Обычно вычисляют таблицу логарифмов. Либо в compile-time, либо можно написать генератор.
            Умножение — это потенцированная сумма логарифмов, по модулю N-1. Недостаток — само умножение в compile-time загнать сложнее (но в принципе возможно).
            0
            "Galua" — это от души. Бедный Galois.
              0
              Не все знакомы с французским
              0
              А еще, чтобы не писать сишные GaloisMul, можно слепить класс с переопределенными арифметическими операциями. Хороший компилятор умеет экземпляр такого класса размещать в регистрах, если нужно. Делал так в рабочем проекте. Всё это крутилось на платформе CUDA с умопомрачительной скоростью :)
                –1
                >Итак, требуется реализовать следующие операции в поле GF(256) над многочленом x^8 + x^4 + x^3 + x + 1:
                Обозначение поля странное. Если числовое поле вычетов, то GF(256), то в скобках должно быть простое число (не 256), если поле расширения, то в скобках степень простого GF(2^8). Не все вещи можно игнорировать, если не хотите демонстрировать свою неграмотность
                Поле задается не только неприводимым многочленом, но и примитивным элементом.
                У Вас он какой? Не встретил о нем в тексте ни одного слова.
                  0

                  GF(n) — это обозначение для поля Галуа из n элементов. Поскольку для каждого n оно такое единственное (с точностью до изоморфизма и если вообще существует), нет никакой разницы в каком виде записано число n.


                  Можно даже GF(255+1) написать, это будет так же корректно.

                Only users with full accounts can post comments. Log in, please.