Group 1 (2).png
Эварист Галуа, умер на дуэли на 21 году жизни, его работы являются фундаментом современной алгебры, а Intel в 2021 году выпустил набор процессорных расширений для работы с полями Галуа.

"Попросите Якоби или Гаусса публично высказать своё мнение — не о истинности, а о важности этих теорем. Позже, я надеюсь, найдутся люди, которым будет выгодно разобраться во всём этом хаосе."

Этими словами заканчивалось письмо Эвариста Галуа, написанное для своего друга Огюста Шевалье за два дня до его смерти от полученных на дуэли ран на 21 году жизни. Ни Якоби, ни Гаусс в его теоремах не разобрались, зато спустя 15 лет разобрался Жозеф Лиувилль и опубликовал работы Галуа, ставшие впоследствии фундаментом современной алгебры, известные сейчас как теория Галуа. В статье расскажу про одну из частей этой теории - поля Галуа, получившая настолько повсеместное применение в криптографии и избыточном кодировании, что Intel и AMD выпустили набор процессорных расширений для эффективной реализации операций над этими полями.

Заметка! Если вам довелось использовать/реализовывать поля Галуа, то большая часть статьи для вас скорее всего будет не интересна, но возможно в последних разделах будет что-то для вас новое.

Что в математике называют полем?

Простыми словами поле - это множество чисел, которые можно складывать, вычитать, умножать и делить, все операции при этом обратимые. Например полями являются рациональные, вещественные и комплексные числа, но например целые и натуральные не являются из-за невозможности деления в них (обычное деление целых чисел приводит к рациональным числам, а деление с остатком необратимо). Более формально поле - это множество \mathbb{F} и две бинарные операции +,\times:\mathbb{F}\times\mathbb{F}\rightarrow\mathbb{F}(каждая операции сопоставляет паре элементов множества другой элемент), удовлетворяющее следующим свойствам

  • Ассоциативность

a+(b+c)=(a+b)+c \\ a \times (b\times c)=(a\times b)\times c
  • Коммутативность

a+b=b+a \\ a\times b=b\times a
  • Дистрибутивность

a \times (b+c)=a\times b + a\times c
  • Существование нуля (нейтральный элемент относительно сложения)

a+0=a, ~~ a\times 0=0 \ \forall a\in\mathbb{F}
  • Существование обратного по сложению элемента

\forall a\in\mathbb{F}~\exists b\in\mathbb{F}: a+b=0
  • Существование единицы (нейтральной элемент относительно умножения)

a\times 1=a~~\forall a\in\mathbb{F}
  • Существование обратного по умножению элемента

\forall a\neq 0\in\mathbb{F}~\exists b\in\mathbb{F}: a\times b=1

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

Модульная арифметика

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

Лемма 1 (о делении с остатком). Для любых целых чисел a, b существуют единственная пара целых чисел q, r таких, что 0\leq r<bи

a=qb+r

Лемма стандартная, оставлю без доказательства, с её помощью однозначно задаётся операция деления с остатком: q- частное, r - остаток. Стоит отметить, что такое определение принято в математике, но в программирование обычно деление немного отличается, остаток может получиться отрицательным например при делении отрицательного числа на положительное, в этом случае можно получить положительный остаток если прибавить b.

Если обозначить за a \bmod bостаток от деления a на b, то оказывается, что

ac\bmod b=(a\bmod b)(c\bmod b)\bmod b

ну или если проговорить это словами, то "остаток от деления произведения чисел равен остатку от деления произведения остатков этих чисел". Если всё еще не совсем понятно, то вот частный пример: последняя цифра произведения чисел такая же как у произведения последних цифр этих чисел. Доказывается довольно просто, по лемме о делении с остатком

a=qb+r \\ c = eb+f

и соответственно

ac=(qeb+qf+er)b+rf

Левое слагаемое правой части делится на b, поэтому не влияет на остаток. Аналогичное свойство работает и с суммой. В итоге получается, что остатки можно складывать и умножать. А можно ли делить? Оказывается что можно, но не всегда.

Лемма 2 (тождество Безу). если a, b взаимно простые, то существуют целые x, y такие, что

ax+by=1

Доказывается конструктивно через расширенный алгоритм Евклида. Что это даёт? А вот что

Следствие. Если p - простое число, 0\leq a < p, то существует 0\leq x<p такое, что

ax\bmod p=1

Для доказательства нужно просто применить тождество Безу к a, p. Возвращаясь к аксиомам поля это даёт нам последнее недостающее звено: обратный элемент по умножению.

Пример. Поле остатков по модулю 7 состоит из элементов \{0, 1, 2, 3, 4, 5, 6\} c таблицей умножения

\begin{array}{ccccccc} \times & 0 & 1 & 2 & 3 & 4 & 5 & 6 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 0 & 1 & 2 & 3 & 4 & 5 & 6 \\ 2 & 0 & 2 & 4 & 6 & 1 & 3 & 5 \\ 3 & 0 & 3 & 6 & 2 & 5 & 1 & 4 \\ 4 & 0 & 4  & 1 & 5 & 2 & 6 & 3 \\ 5 & 0 & 5 & 3 & 1 & 6 & 4 & 2 \\ 6 & 0 & 6 & 5  & 4 & 3 & 2 & 1 \end{array}

А обратные по умножению элементы разбиваются на 1=1\times 1=2\times 4=3\times 5=6\times 6. Таким образом это полноценное поле.

От древней Греции до наших дней

Считается, что описание модульной арифметики в современном виде сделал Гаусс, однако очень многие сопутствующие факты были известны задолго до него:

  • Расширенный алгоритм Евклида для решения диофантовых уравнений. Авторство нахождение наибольшего общего делителя однозначно отдают Евклиду, но вот с его расширенной версией не всё так понятно, но в целом справедливо будет назвать его авторами Евклида и Диофанта. Известно ровно два эффективных метода для нахождения обратного элемента в поле остатков по модулю -- расширенный алгоритм Евклида и возведение в степени на основании малой теоремы ферма и теоремы Эйлера, Евклид побеждает по производительности (пишите в комментариях это не так или если на самом деле есть алгоритмы получше).

  • Теоремы Сунь Цзы (китайская теорема об остатках) о том, что на интервале [0, p_1p_2\ldots p_n)существует единственное решение линейной системы сравнений

a_1\bmod p_1=b_1 \\ a_2\bmod p_2=b_2 \\ \vdots \\ a_n\bmod  p_n=b_n
  • Теорема датируется 13-м веком, на просторах интернета можно найти следующее применение, о достоверности которого судить не готов: у генерала батальон примерно на 1000 человек, после очередного боя есть потери, нужно их быстро оценить, пересчитывать каждого по очереди долго, вместо этого можно быстро построить солдат сначала в 3 шеренги, потом в 5, и наконец в 7, попутно запоминая сколько оставалось в последнем ряду. Если потерь было не больше 3\cdot 5\cdot 7=105, то этой информации достаточно для того, чтобы эти потери точно посчитать. В современности же на этой теореме основаны например системы остаточных классов: если нужно провести какие-то вычисления на больших числах, то можно сделать провести эти вычисления по модулю различных простых чисел, а потом объединить результат -- в некоторых ситуациях таким образом можно ускорить вычисления, особенно если исходный алгоритм плохо распараллелен.

  • Малая теорема Ферма о том, что для простого числа p и не делящегося на p числа a выполняется равенство

a^{p-1}\bmod p=1
  • Теорема Эйлера о том, что если обозначить \varphi(b)- количество натуральных чисел меньше b, являющихся с b взаимно простыми, то

a^{\varphi(b)}\bmod b=1

Все эти факты лежат в основе современной криптографии, в частности алгоритм Диффи-Хеллмана использует возведение в степени по модулю и опирается на сложность решения задачи дискретного логарифмирования в качестве гарантии надёжности шифрования, а алгоритм RSA опирается на теорему Эйлера и факт о том, что для вычисления \varphi(b) для составного числа b нужно уметь раскладывать его на множители -- гораздо более сложная задача, чем возведение в степень. Алгоритм Шёнхаге-Штрассена для умножения чисел основан на адаптации алгоритма Кули-Тьюки быстрого преобразования Фурье для конечных полей по модулю простых Ферма - это простые числа вида 2^m+1, удобство которых в существовании первообразного корня из единицы \omega степени 2^m.

А что собственно придумал Галуа?

Один из вопросов, которым задался Галуа - какие еще бывают поля с конечным числом элементов кроме полей остатков по простому модулю? Во-первых, стоит отметить, что остатки по составному модулю не образуют поля, например у 2 нет обратного по модулю 4, в общем случае делители модуля (кроме единицы) не будут иметь обратного элемента. Галуа доказал, что существуют поля, состоящие из p^n элементов, где p простое, а n натуральное; доказал, что других нет; и наконец доказал, что два поля, состоящие из одного числа элементов структурно идентичны. Более формально

Теорема Галуа о конечных полях. Для любого простого числа p и натурального числа n

  • Существует поле \mathbb{F}, состоящее ровно из p^n элементов

  • Если натуральное числа q не является степенью простого числа, то не существует поля, состоящего из q элементов

  • Для любого другого поля \overline{\mathbb{F}} также состоящего из p^n элементов существует изоморфизм \varphi: \overline{\mathbb{F}}\rightarrow \mathbb{F} такой, что

\varphi(a+b)=\varphi(a)+\varphi(b) \\ \varphi(a\times b)=\varphi(a)\times\varphi(b)

Здесь я ограничусь только частичным доказательством первого пункта, так как он нам позже понадобится. Итак, для n=1мы знаем, что поле существует - поле остатков по модулю, обозначим его F. Теперь допустим n>1, давайте посмотрим на множество многочленов степени меньше n

F_n[x]=\{a_0+a_1x+a_2x^2+\ldots +a_{n-1}x^{n-1}~|~a_i\in F\}

Определим сумму как обычно, т.е. отдельно по коэффициентам

\sum_{i=0}^{n-1}a_ix^i+\sum_{i=0}^{n-1}b_ix^i=\sum_{i=0}^{n-1}(a_i+b_i)x^i

А вот для умножения в F_n[x] сделаем хитрость: возьмем некоторый многочлен g(x) степени ровно n и при умножении в F_n[x]будем делать обычное умножение многочленов, после чего брать остаток при делении на g(x)

\sum_{i=0}^{n-1}a_ix^i\times \sum_{i=0}^{n-1}b_ix^i=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_ib_jx^{i+j}\bmod g(x).

По аналогии с числами, остаток при делении a(x)на b(x) - это единственный многочлен r(x) степени меньше b(x)такой, что для некоторого q(x)

a(x)=q(x)b(x)+r(x)

Так как остаток при делении на многочлен g(x) степени n будет иметь степени меньше n, то результат такого умножения остается в множестве F_n[x]. Корректность всех свойств поля кроме существования обратного элемента легко проверяется, давайте разбираться с ним.

Возьмем произвольный элемент a\in F_n[x]и посмотрим на последовательность степеней a: 1, a, a^2, a^3, \ldots В этой последовательности рано или поздно должно произойти одно из трёх:

  • Очередная степень станет равной нулю и все последующие тоже обнуляться

  • Очередная степень станет равной единице и дальше последовательность зациклится

  • Очередной степень будет равна другой уже присутствующей ранее, но не единице, последовательность зациклица, единица будет только на первом месте

Сначала посмотрим на последнюю ситуацию. Допустим повторилась степень a^k=a^l, k<l, раз это первое повторение, то значит, что a^{l-k}\neq 1, и при этом a^ka^{l-k}=a^kили a^k(a^{l-k}-1)=0. Оставим это пока и посмотрим на первую ситуацию, допустим оказалось, что a^k=0, т.е. в частности получается, что a\times a^{k-1}=0. В обоих случаях мы получаем обычное произведение двух многочленов степени меньше n делится на g(x). А теперь давайте представим, что изначально мы выбрали g(x) таким образом, что у него нет делителей в F_n[x], тогда вышеописанные ситуации становятся невозможными. Остаётся второй случай когда очередная степени равна единице, т.е. a^k=1или a\times a^{k-1}=1, т.е. a^{k-1}- обратный элемент к a. Таким образом, чтобы конструкция образовывала поле нужен многочлен g(x), не имеющий делителей в F_n[x]- оказывается, что такой многочлен всегда существует, доказательство этого факта приводить не буду.

Конечные поля из q элементов часто обозначают как GF(q) от Galois field.

Поля Галуа в современности: Rijndael и Reed-Solomon

Есть два крупных распространённых применения полей Галуа, которые сводятся к тому, что естественные для компьютера единицы информации - последовательности из k байт - можно интерпретировать в виде элемента GF(2^{8k}) и проводить над ним какие-то преобразования, которые можно делать только в поле.

Advanced Encryption Standard (AES). Принятый NIST стандарт шифрования, основанный на методе Rijndale, особо ничего про него сказать не могу кроме того, что одна из его составных частей - это несколько побайтовых умножений в GF(256). Проконсультировался с коллегами алгебраистами, так если честно и не понял почему и зачем именно такое преобразование используется, напишите в комментариях если знаете зачем. И тем не менее факт в том, что часть AVX-512 сделана специально под этот стандарт для работы с GF(256).

Reed-Solomon error correction. Общая идея довольно простая. Представьте, что ваши данные можер представить в виде прямой на плоскости, и вы хотите эти данные передать по ненадёжной сети, где время данные могут потеряться. Чтобы передать вашу прямую достаточно двух точек на этой прямой - из них саму прямую можно однозначно восстановить. Но так как у нас могут потеряться данные, то можно вместо двух точек взять скажем 5 (разных разумеется), тогда чтобы восстановить прямую достаточно, чтобы не потерялись любые две из них. Аналогично например с кругом, который однозначно задаётся тремя точками.

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

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

Оказывается, что есть обобщение такого подхода: кривая задаётся уравнением y=a_0+a_1x+a_2x^2+\ldots a_{n-1}x^{n-1}, то из интерполяционной теоремы мы знаем что любые n точек на этой кривой однозначно выделят её среди других кривых в этом классе. Соответственно алгоритм передачи данных будет выглядеть следующим образом

  • Хотим передать последовательность a_0, a_1, \ldots, a_{n-1}

  • Строим многочлен P(x)=a_0+a_1x+\ldots+a_{n-1}x^{n-1}

  • Вычисляем и отправляем P(1), P(2), \ldots, P(n+k)

  • На стороне получателя получили величины с предыдущего шага, что-то потерялось

  • Если получили хотя бы n из значений, то проводим интерполяцию и получаем P(x), его коэффициенты - это передаваемое сообщение

  • Замечание! Так как передаём мы только P(i), но не сами i, то чтобы провести интерполяцию на последнем шаге необходимо знать соответствие i\rightarrow P(i), которое обычно известно косвенно из протокола.

Вышеописанная схема почти соответствует простейшей схеме алгоритма Рида-Соломона. Проблема заключается в том, что если всё вышеперечисленное мы делаем над обычными вещественными числами, то из-за ошибок округления эту схему практически нигде не применить. Последний кусочек пазла заключается в том, что всё это можно проделывать в конечном поле, где проблемы ошибок округления нет. В частности последовательность из n байт может быть однозначно интерпретирована как коэффициенты многочлена степени меньше n, P(i) является элементом поля и соответственно также легко представляется в виде байта. Итоге мы получаем протокол, который берёт последовательность из n байт, превращает её в последовательность из n+k байт, передаёт её по сети и восстанавливает исходную последовательность если дошли любые из n переданных байт.

Особенности реализации полей размера 8k

В основном я буду писать про наиболее распространённое GF(2^8), но обычно верно и для других полей, какие-то уточняющие моменты буду подсвечивать отдельно.

Итак, для построения поля GF(2^8) будем использовать представление через многочлены, описанное два раздела назад. Нам понадобится

1. Выбрать способ представления байта в виде элемента GF(256)

Напомню, что элементы GF(2^8) - это многочлены степени меньше 8 с коэффициентами 0 или 1, так это же просто 8 битов, давайте хранить коэффициент перед степенью i в i-ом бите байта, например

uint8_t element = 0b00110101; // x^5+x^4+x^2+1

2. Научиться складывать и вычитать

uint8_t Add(uint8_t a, uint8_t b) {
  return a ^ b;
}

uint8_t Subtract(uint8_t a, uint8_t b) {
  return a ^ b;
}

Так, сложение и вычитание у нас определялось просто как обычное сложение и вычитание для многочленов, т.е. по каждой степени по отдельности, а в GF(2) это просто XOR, соответственно реализация сложения/вычитания очень проста

Ну и сразу отметим, что в 0 является нейтральным по отношению с сложению, а любой элемент является обратным к самому себе по сложению так как x ^ x = 0.

3. Научиться умножать

Так, вот здесь немного хитрее, для начала нам нужен неприводимый в GF(2) многочлен степени 8. За нас их уже посчитали добрые алгебраисты, их всего 30 выбирайте любой, я возьму вот этот x^8+x^4+x^3+x+1, обозначим его младшие степени в виде отдельного элемента.

const uint8_t irreducible_poly = 0x1b; // x^4+x^3+x+1

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

uint8_t Multiply(uint8_t a, uint8_t b) {
  uint8_t result = 0;
  while (a) {
    result ^= b * (a & 1);
    a >>= 1;
    // b << 1 соответствует домножению на многочлена на x
    // (b >> 7) соответветствует проверки, что b
    // многочлен степени 7.
    // ^ (irreducible_poly * (b >> 7)) соответствует
    // взятию по модулю в этом случае. В противном случаем
    // для взятия по модулю ничего не нужно делать
    b = (b << 1) ^ (irreducible_poly * (b >> 7));
  }
  return result;
}

Здесь хочется отметить, что все вычисления производятся над 8 битами, стоит отметить, что b<<1 в реальности вылезет на один бит, но он обрубится при присвоении. Если не совсем понятно почему этот код делает ровно то, что заявлено, добавлю еще аналогичный код обычного умножения в столбик двух чисел по модулю

int MultiplyMod(int a, int b, int mod) {
  int result = 0;
  while (a) {
    result += b * (a & 1);
    a >>= 1;
    b = b << 1; // b <-- 2b
    // До предыдущей операции выполнялось 
    // b < mod, а значит 2b < 2mod-1
    // и для mod достаточно одного вычитания
    if (b > mod)
      b -= mod;
  }
  return result;
}

4. Научиться находить обратный элемент и делить

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

Лемма 3 (следствия из теоремы Лагранжа для групп). Для любого ненулевого элемента поля a\in F_{p^n} выполняется

a^{p^n-1}=1

Более того, существует элемент \alphaтакой, что степени 1, \alpha, \alpha^2, \ldots, \alpha^{p^n-2}различны и соответствуют всем элементам поля F_{p^n} кроме нуля (такой элемент принято называть примитивным элементом поля).

Таким образом, найти a^{-1} в поле GF(2^8) можно например как a^{254}.

uint8_t Inverse(uint8_t a) {
  uint8_t result = 1;
  uint8_t pow = 254;
  while (pow) {
    if (pow & 1)
      result = Multiply(result, a);
    a = Multiply(a, a);
    pow >>= 1;
  }
  return result;
}

На этом базовая реализация заканчивается …

LUT для умножения и обращения

Для начала стоит отметить, что таблицу умножения для GF(2^{8k}) можно просто посчитать в виде двумерного массива размера 2^{16k}, что в целом годится для k=1, но для k=2 уже неприменимо. Но есть более элегантный способ, у нас есть примитивный элемент \alpha, для GF(2^8) с порождающим многочленом x^8+x^4+x^3+x+1 это например x+1 (3 в битовой записи), любой элемент представим в виде степени примитивного, соответственно если a=\alpha^n, b=\alpha^m, то a\times b=\alpha^{n+m}. И при этом мы знаем, что \alpha^{255}=1. Из этих свойств возникает следующая идея

  • Подсчитаем таблицу степеней \alpha, т.е. \exp[i]=\alpha^i

  • Подсчитаем таблицу логарифмов \alpha, т.е. \alpha^{\log[x]}=x.

  • Получаем a\times b=\exp[\log[a]+\log[b]\bmod255]

uint8_t exp[256];
uint8_t log[255];
const uint8_t primitive_element = 3;

uint8_t InitTables() {
  uint8_t x = 1;
  for (size_t i = 0; i < 255; i++) {
    exp[i] = x;
    log[x] = i;
    x = Multiply(x, primitive_element);
  }
  exp[255] = 1;
}

uint8_t MultiplyLUT(uint8_t a, uint8_t b) {
  if (a == 0 || b == 0)
    return 0;
  auto p = log[a] + log[b];
  // Во всех случаях кроме log[a]+log[b]=255
  // получаем эквивалент p % 255, если же
  // log[a]+log[b]=255, то пользуемся тем, что
  // exp[255] = 1
  return exp[(p & 255) + (p >> 8)];
}

Такие таблицы занимают 2\times 255байт, аналогичные для GF(2^{16})занимала бы 2\times 2\times (2^{16}-1)байт.

Наконец стоит отметить, что для обращения эти таблица тоже можно использовать

uint8_t InverseLUT(uint8_t a) {
  if (a == 0)
    return 0;
  return exp[255 - log[a]];
}

Построение через башни

Указанная выше конструкция с построением через многочлены работает для любого базового поля. Мы делали расширение GF(2) до GF(2^8), но так же например можно делать расширение GF(2^8) до GF(2^{8k}), для этого нужно проделать аналогичную конструкцию с многочленами степени меньше k с коэффициентами из GF(2^8). Технически удобно строить поля последовательностью расширений, например каждое из полей GF(2), GF(2^2), GF(2^4), GF(2^8), GF(2^{16})можно получить из предыдущего расширением многочленами первой степени. Вот пример как можно получить последнее расширение. Нам понадобится неприводимый многочлен степени 2 над GF(2^8), например многочлен вида

x^2+x+\delta

в половине случае является неприводимым (верно для GF(2^k)), возьмём \delta=32. Элементы нашего расширения имеют вид

a_0+a_1x\in GF(2^{16}), ~a_0,a_1\in GF(2^8)

Со сложением как и раньше все просто - это побитовый XOR. А что с умножением? Допустим пытаемся умножить a на b, для этого как и раньше умножить многочлены как обычно, а потом поделить с остатком на x^2+x+\delta, здесь удобно использовать трюк, что это деление эквивалентно замене x^2=x+\delta.

\begin{array}{rl} (a_0+a_1x)(b_0+b_1x)&=a_0b_0+(a_0b_1+a_1b_0)x+a_1b_1x^2\\&=a_0b_0+(a_0b_1+a_1b_0)x+a_1x_1(x+\delta) \\&=a_0b_0+a_1b_1\delta+(a_0b_1+a_1b_0+a_1b_1)x \end{array}

А вот как это будет выглядеть в коде

uint16_t Multiply(uint16_t a, uint16_t b) {
  // a = a_0 + a_1x, b = b_0 + b_1x
  // все отдельные элементы из GF(256)
  uint8_t a_0 = a & 255;
  uint8_t a_1 = a >> 8;
  uint8_t b_0 = b & 255;
  uint8_t b_1 = b >> 8;

  auto t = gf_2_8::MultiplyLUT(a_1, b_1);
  auto low_bits =
      gf_2_8::Add(gf_2_8::MultiplyLUT(a_0, b_0), gf_2_8::MultiplyLUT(t, delta));
  auto high_bits = gf_2_8::Add(
      gf_2_8::Add(gf_2_8::MultiplyLUT(a_0, b_1), gf_2_8::MultiplyLUT(a_1, b_0)),
      t);
  return low_bits + (high_bits << 8);
}

Алгоритм нахождения обратного элемента Itoh–Tsujii

Можно нахождение обратного элемента можно сделать через возведение в степень, но есть алгоритм, позволяющий сократить количество умножений. В общем виде идея выглядит так: в поле GF(q^n). Если обозначить r=(q^n-1)/(q-1), то

1=a^{q^n-1}=a\cdot \underbrace{a^{r-1}a^{r(q-2)}}_{a^{-1}}

Вот здесь есть два интересных момента:

  • a^r \in GF(q). Без доказательства

  • Так как a^r \in GF(q), то a^{r(q-2)}=(a^r)^{-1}, инверсия происходит в поле GF(q)

Применительно к полям GF(2^8) и GF(2^{16}) получаем следующий алгоритм

  • 2^{16}-1=255\cdot 257, r=257

  • Вычисляем a^{256}и a^{257}в GF(2^{16})

  • Вычисляем (a^{257})^{-1}в GF(2^8)

  • a^{-1}=a^{256}\cdot (a^{257})^{-1}

Конкретно в данном случае мы уменьшили количество умножений в GF(2^{16})почти в 4 раза: вычисление влоб

a^{2^{16}-2}=\underbrace{a^{2^1}\cdot a^{2^2}\cdot \ldots\cdot a^{2^{15}}}_{14+15~умножений}

А алгоритму Itoh-Tsujii нужно 8 умножений для подсчета a^{2^8}и еще одно для a^{257}.

uint16_t InvIT(uint16_t a) {
  uint16_t a_r = a;
  for (size_t i = 0; i < 8; ++i) {
    a_r = Multiply(a_r, a_r);
  }
  uint8_t a_r1 = Multiply(a_r, a);
  // Подсчет обратного элемента из GF(256)
  return Multiply(a_r, gf_2_8::Inv(a_r1));
}

Intel Galois Field New Instruction (GFNI)

Из-за популярности и важности GF(2^8) интел выпустили набор расширений для вычислений в этом поле, произошло это кстати не так давно, в 2021 году. Общий принцип этого расширения заключается в следующем. Поле GF(2^8)в том числе является линейных пространством размера 8 над полем GF(2), функция умножения на элемент a

f_a(x)=a\times x

линейна по свойству поля. Из этого следует, что эта функция может быть представлена в виде домножения на матрицу размера 8\times 8. Столбцы этой матрицы можно получить например домножением a на базис 1, x, x^2, x^3, x^4, x^5, x^6, x^7. Следующий код демонстрирует сам концепт, но подобная реализация средствами С++ не будет также быстра как GFNI

void InitGFNI(void) {
  for (int16_t y = 0; y < 256; ++y) {
    gfni_matrix[y] = 0;
    element_t row = y;
    for (size_t i = 0, shift = 0; i < 8; ++i, shift += 8) {
      gfni_matrix[y] |= ((uint64_t)row << shift);
      row = (row << 1) ^ ((row >> 7) * irreducible_poly);
    }
  }
}

element_t MultiplyGFNI(element_t a, element_t b) {
  return ((a & 1) * (gfni_matrix[b] & 255)) ^
         (((a >> 1) & 1) * ((gfni_matrix[b] >> 8) & 255)) ^
         (((a >> 2) & 1) * ((gfni_matrix[b] >> 16) & 255)) ^
         (((a >> 3) & 1) * ((gfni_matrix[b] >> 24) & 255)) ^
         (((a >> 4) & 1) * ((gfni_matrix[b] >> 32) & 255)) ^
         (((a >> 5) & 1) * ((gfni_matrix[b] >> 40) & 255)) ^
         (((a >> 6) & 1) * ((gfni_matrix[b] >> 48) & 255)) ^
         (((a >> 7) & 1) * ((gfni_matrix[b] >> 56) & 255));
}

Интересно, что сами разработчики Intel в технической документации предлагают в качестве применения для такого подхода например перетасовывание битов в числе. Забавно, что эту идею реализовали в Clang 19 где __builtin_bitreverse64 компилируется в инструкцию vgf2p8affineqb в случае её доступности.

Ссылки


Друзья и коллеги! С удовольствием хотел бы прорекламировать CS Space — открытый научный клуб по CS-related темам; идейных последователей питерского Computer Science ClubCS Center), расформировавшегося после известных событий. Ребята организуют крутые лекции и курсы по CS от профессионалов своего дела, да еще и помогают мне с написанием научно-популярных статей!

Сайт сообщества: csspace.io
Telegram-канал: t.me/csspace

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

Комментарии (0)