Программная реализация умножения в полях Гаула

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

И так, я отвлёкся. В детстве-юношестве для помехоустойчивого кодирования можно было бы применить контроль чётности по матричному методу, но сейчас это не серьёзно. Полистав интернет я решил остановиться на кодировании по методу Рида-Соломона. Алгоритм не то, чтобы совсем новый, его ещё в первых CD применяли, но при этом, насколько мне известно, не потерявший своей актуальности и на данный момент. В этой статье о самих кодах Рида-Соломона я не буду сильно распространяться – это за меня на Хабре сделали много раз и много кто. Здесь я хочу описать реализацию алгоритма умножения в GF[256]. Тем не менее, чтобы не заставлять читателя прыгать по ссылкам, кратенько опишу зачем вообще приходится иметь дело
с полями Гаула.

Избыточное кодирование Рида-Соломона это о матрицах. Мы используем матрицы при кодировании и при декодировании. В этих процессах присутствуют как операции умножения матриц, так и операции нахождения обратных матриц — деления. Деление здесь требуется не приблизительное-целочисленное, а самое настоящее, точное. А точное деление для компьютера нерешаемая задача: один разделить на три это ноль целых и бесконечное количество троек после запятой, но память при этом 640 килобайт хватит каждому.

Гаул жил в начале 19 века, жил очень мало (21 год, вообще про его личность история интересная, но сейчас не об этом). Его работы очень сильно пригодились при цифровом кодировании. Конечное поле Гаула — это набор чисел от нуля до n. Суть и интересность этих полей в том, что для элементов этого набора можно определить операции сложения-умножения-вычитания-деления такие, что результат операции будет находиться в самом этом поле. Например, возьмём набор (поле):

$\$\$display\$\$[0,1,2,3,4,5,6,7]\$\$display\$\$$

2+2=4 в этом поле также, как и в поле привычных нам, действительных чисел. А вот 5+6 равно не 11, как мы привыкли, а 5+6=3, и это замечательно! 3 входит в это поле (вообще сложение и вычитание для конкретно этого поля это просто побитовый XOR). Для умножения и деления тоже выполняется правило: результат от умножения или деления любых двух чисел из поля (набора… Далее буду говорить только «поле») будет находиться в этом поле.

Вообще, количество элементов поля это число не произвольное. Оно называется «характеристикой поля» и может быть либо простым числом, либо степенью простого числа (её называют «порядком поля», основание этой степени называют «основанием поля»). К счастью, количество чисел, которое можно зашифровать с помощью одного байта равно 256, а это два (простое число) в степени 8, и мы можем принять множество всех возможных значений байта за конечное поле Гаула. По умному оно называется GF[256], GF значит Galois Field.

Так вот. Если использовать арифметику в поле Гаула над байтами, то при делении матриц, состоящих из байтов никогда не получится чисел, которые нельзя записать в один байт, а значит то, как реализовывать «в железе» (я использую stm32 в основном) алгоритмы кодирования-декодирования становится немного понятнее.

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

$\$\$display\$\$A+B\; =\; A\; xor\; B\$\$display\$\$$

$\$\$display\$\$A-B\; =\; A\; xor\; B\$\$display\$\$$

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

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

Пример:

$\$\$display\$\$5=101\_2=1*x^2+0*x^1+1*x^0=x^2+1\$\$display\$\$$

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

$\$\$display\$\$x^2+x^2=0\$\$display\$\$$

$\$\$display\$\$x^2+x^2+x^2=x^2\$\$display\$\$$

$\$\$display\$\$x^2+x^2+x^2+x^2=0\$\$display\$\$$

То есть коэффициенты складываются по модулю 2. Пример умножения в GF[8]:

$\$\$display\$\$6*3\; =\; 110\_2*11\_2=\; (1*x^2+1*x^1+0*x^0)*(1*x^1+1*x^0)=\$\$display\$\$$

$\$\$display\$\$\; =(x^2+x)(x+1)=x^2*x+x^2*1+x*x+x*1=\$\$display\$\$$

$\$\$display\$\$=x^3+x^2+x^2+x\$\$display\$\$$

Далее «складываем» (в кавычках — потому, что по модулю 2) подобные члены, и получаем $\$inline\$x^3+x\$inline\$$, что мы переводим в обычное число $\$inline\$1010\_2=10\$inline\$$.

Так как это умножение мы подразумевали в GF[8], то число 10 в результате должно вызвать недоумение: «Как же так? 10 нету в наборе [0,1,2,3,4,5,6,7]». И да. Тут кроется ещё ложка дёгтя. Если получается такая ситуация (результат перемножения не в поле), то далее мы проводим операцию деления по правилам полиномов. Причём делим мы не на что попало, а на определённый, так называемый «порождающий» полином. Откуда он берётся? Почему делится только сам на себя? И прочие вопросы я оставлю без ответа (не потому, что жадный, а потому что не хочу вводить в заблуждение, сам плохо понимаю). Скажу лишь, что для конкретного поля Гаула есть один или несколько таких подходящих порождающих полиномов. Для GF[8] их всего два: 11 и 13, то есть 1011 и 1101 или $\$inline\$x^3+x+1\$inline\$$ и $\$inline\$x^3+x^2+1\$inline\$$

И так, чтобы умножить 6 на 3 в поле GF[8] недостаточно лишь перемножить два многочлена. После перемножения надо результат $\$inline\$x^3+x\$inline\$$ поделить на один из возможных порождающих многочленов. Здесь мы выберем $\$inline\$x^3+x+1\$inline\$$. Правила деления многочленов столбиком здесь похожи на те, которые мы должны помнить из школьной программы, но во-первых «должны», а во вторых они всё же отличаются. Поясню на нашем примере. Берём старшую степень делимого, в нашем случае это $\$inline\$x^3\$inline\$$ и делим на старшую степень делителя (здесь это также $\$inline\$x^3\$inline\$$). Записываем результат 1. Далее получившийся результат умножаем на весь делитель (получается $\$inline\$x^3+x+1\$inline\$$). Это мы записываем под делимым, и проводим вычитание (а в поле GF[8] вычитание и сложение есть одно и то же) то есть $\$inline\$(x^3+x)\; +(x^3+x+1)\$inline\$$. Сложение проводим по модулю 2, этому нас обязывает работа в поле, и получаем 1. Далее этот результат мы опять бы делили на порождающий многочлен таким же образом, если бы степень получившегося результата (после вычитания) была бы больше или равной степени делителя (порождающего многочлена). А так как степень меньше, то мы за результат деления принимаем остаток от деления, то есть результат вычитания(сложения – здесь это одно и тоже), то есть 1.

Всё. 6 * 3 = 1 для GF[8] c порождающим полиномом 11.

Казалось бы всё непросто. С такими сложностями нужно много-много кода и суперкомпьютеры и вообще проще хранить таблицу умножения как массив 256х256 (ну или 8х8, смотря какое поле мы используем), Но всё не так плохо. У полей Гаула есть ещё одно замечательное свойство, и связано оно с операцией возведения в степень. Так как возведение в степень — это просто последовательное умножение, то результат возведения в степень так же является элементом поля Гаула. Также для любого поля верно утверждение, что в нём есть минимум один элемент, степени которого содержат в себе всё поле (кроме 0). То есть, если взять для GF[8] число 2, то
$\$inline\$2^0=1\$inline\$$      $\$inline\$2^7=1\$inline\$$
$\$inline\$2^1=2\$inline\$$      $\$inline\$2^8=2\$inline\$$
$\$inline\$2^2=4\$inline\$$      $\$inline\$2^9=4\$inline\$$
$\$inline\$2^3=3\$inline\$$      $\$inline\$2^\{10\}=3\$inline\$$
$\$inline\$2^4=6\$inline\$$      $\$inline\$2^\{11\}=6\$inline\$$
$\$inline\$2^5=7\$inline\$$      $\$inline\$2^\{12\}=7\$inline\$$
$\$inline\$2^6=5\$inline\$$      $\$inline\$2^\{13\}=5\$inline\$$

Если присмотреться, то можно обратить внимание, что любое число от 1 до 7 можно однозначно представить как 2 в какой либо степени. Так же свойство операции возведения в степень в этом поле таково, что $\$inline\$2^7=2^0\$inline\$$ ещё $\$inline\$2^8=2^1\$inline\$$, $\$inline\$2^9=2^2\$inline\$$ и так далее. Это значит, что 2 в любой степени может быть представлено в виде двойки в степени от нуля до 6 с помощью операции получения остатка от деления на 7 в случае положительной степени и простой формулы $\$inline\$2^\{-x\}=2^\{(7-(x\backslash :\; mod\backslash :7))\}\$inline\$$, если показатель — отрицательное число

Собственно, если вспомнить свойство умножения степеней, то умножение чисел многократно упрощается. И чтобы умножить 6 на 3 мы теперь смотрим, в какой степени двойка равна 6 и в какой степени двойка равна 3, складываем степени и смотрим результат — два в какой-то степени, которую можно представть как 2 в степени от 0 до 6 пример:

$\$\$display\$\$6*3=2^4*2^3=2^\{(4+3)\}\; =\; 2^7\; =\; 2^0\; =\; 1\; \$\$display\$\$$

С делением тоже всё становится предельно понятно:

$\$\$display\$\$3/6=2^3/2^4=2^\{(3-4)\}\; =\; 2^\{-1\}\; =\; 2^\{(7-(1\backslash :\; mod\backslash :7))\}\; =\; 2^6\; =\; 5\; \$\$display\$\$$

И вроде бы вот оно! счастье программиста — реализация арифметики в поле Гаула — пара строчек кода, не нужно заморачиваться с полиномами… Да и быстродействие такого кода будет высоким, но тут я столкнулся с ещё одной проблемой: Таблицу степеней двойки в поле GF[8] с порождающим полиномом 11 найти легко. Даже на этом ресурсе есть. А вот таблицу степеней для GF[256] я с нахрапа не нагуглил, так что решил её составить самостоятельно, C# мне в помощь. Пришлось реализовать алгоритм умножения по правилам полиномов.

Привожу рабочую функцию перемножения для GF[2^n] c произвольным полиномом. Ограничение — я испольую 32-битную арифметику, так что n должно быть меньше 16. Также здесь добавлена функция поиска номера старшего бита числа

private uint GetLeadBitNum(UInt32 Val) {
    int BitNum = 31;
    uint CmpVal = 1u << BitNum;
    while (Val < CmpVal) {
        CmpVal >>= 1;
        BitNum--;
    }
    return (uint)BitNum;
}
private uint Galois_b2_ext_mult(uint m1, uint m2, uint Poly) {
    if (0 == m1 || 0 == m2) { return 0; }
    uint m1_tmp = m1;
    uint m2_tmp;
    uint m1_bit_num = 0;
    
    //Перемножение двух полиномов, при использовании арифметики по модулю 2 достаточно простое занятие.
    //перебираем единички и нолики (для каждого бита первого числа перебираем все биты второго (или наоборот)), складываем номера позиций битов,
    //но не всегда, а только когда оба перебираемых бита - единицы, и инвертируем бит результата под номером, равном сумме позиций для данного шага перебора
    //(инверсия - это прибавление единицы по модулю 2)
    uint PolyMultRez = 0;

    while (m1_tmp != 0) {
        uint bit_m1 = (m1_tmp & 1u) == 0u ? 0u : 1u;
        m1_tmp = m1_tmp >> 1;
        m2_tmp = m2;
        uint m2_bit_num;
        m2_bit_num = 0;
        while (m2_tmp != 0) {
            uint bit_m2 = (m2_tmp & 1u) == 0u ? 0u : 1u;
            m2_tmp = m2_tmp >> 1;
            if ((bit_m1 != 0) && (bit_m2 != 0)) {
                int BitNum = (int)(m2_bit_num + m1_bit_num);
                PolyMultRez ^= 1u << BitNum;
            }
            m2_bit_num = m2_bit_num + 1;
        }
        m1_bit_num = m1_bit_num + 1;
    }

    //Тут есть результат умножения полиномов PolyMultRez. Осталось найти остаток от деления на выбранный порождающий полином.
    //Деление полиномов происходит так: Берём старшую степень делимого, и вычитаем старшую степень делителя. 
    //Получаем число - степень частного
    //Теперь перемножаем, а по сути, просто прибавляем к каждой степени делителя степень получившегося частного
    //и повторяем всё по кругу, пока степень делимого не окажется меньше степени делителя
    uint TmpDivisor_lead_bit_n;
    uint TmpQuotient;
    uint TmpDivisor = Poly;
    uint TmpDividend = PolyMultRez;
    uint TmpDividend_LeadBitNum;
    uint TmpMult_bitNum;
    uint TmpMult_rez;

    TmpDividend_LeadBitNum = GetLeadBitNum(TmpDividend);
    TmpDivisor_lead_bit_n = GetLeadBitNum(TmpDivisor);

    while (TmpDividend_LeadBitNum >= TmpDivisor_lead_bit_n) {

        TmpQuotient = (TmpDividend_LeadBitNum - TmpDivisor_lead_bit_n);

        TmpMult_bitNum = 0;
        TmpMult_rez = 0;
        while (TmpDivisor != 0) {
            uint bit_TmpMult = (TmpDivisor & 1u) == 0u ? 0u : 1u;
            TmpDivisor >>= 1;
            TmpMult_rez ^= bit_TmpMult << (int)(TmpQuotient + TmpMult_bitNum);
            TmpMult_bitNum = TmpMult_bitNum + 1;
        }
        TmpDividend = TmpDividend ^ TmpMult_rez;
        TmpDivisor = Poly;
        TmpDividend_LeadBitNum = GetLeadBitNum(TmpDividend);
    }            
    //Результат умножения числел есть остаток от деления произведения многочленов на порождающий полином.
    return TmpDividend;
}

Теперь, с помощью приведённой выше функции можно составить таблицу степеней двойки для нужного мне GF[256] и написать модуль арифметики Гаула для stm32 с использованием двух таблиц по 256 — одна для прямого, вторая для обратного преобразования числа в его степень. К самой реализации кодирования Рида-Соломона я ещё даже не приступал, но, когда будет готово, думаю поделюсь здесь же. Надеюсь, получится короче.