Вариация на тему Рида-Соломона

Введение

В одном из проектов столкнулся с задачей кодирования данных с целью восстановления потерянных пакетов. Поскольку обработка пакетов осуществлялась полностью на цифровом уровне без доступа к информации от аналогового приемника (hard-decision), то я решил использовать код Рида-Соломона (РС). Обработка пакетов осуществлялась на контроллере esp32-s3, который среди прочего имеет возможность работы с векторами. И необходимо иметь большую силу воли, чтобы не воспользоваться этой интересной возможностью для ускорения вычисления. Собственно эта краткая статья посвящена адаптации и модификации кода РС для возможности использования векторных операций на этом контроллере.

Описание проблемы и путь решения

Основная вычислительная сложность применения РС обусловлена сложностью умножения в поле Галуа. Насколько я понимаю, основной подход к реализации этого умножения оснван на использовании табличного метода — либо простая таблица всех попарных умножений, либо совместное использование таблицы логарифма/степени базового элемента. К сожалению, табличные методы не векторизуются на esp32-s3, поскольку обращение к памяти не относится к списку доступных векторных операций. А выполнение умножения через побитный цикл сведет к нулю выйгрыш от векторных операций.

Я не буду углубляться в построение РС. Общая задача систематического кодера — используя подходящую матрицу коэффициентов, умножить эту матрицу на вектор входных слов. В результате получается вектор проверочных символов, который при передаче добавляется к исходным словам. Задача декодера, используя все полученные слова, исправить/восстановить исходный вектор входных слов. В моем случае я использовал все 8 битные слова.

Для адаптации РС требуется следует опираться на конкретные требования и особенности моей задачи. Во-первых, мне не требуется искать испорченные данные, я знаю номера потерянных пакетов данных, и требуется лишь их восстановить. Это значительно упрощает декодер и тем самым расширяет возможности вариации кодировщика. Во-вторых, мне не требуется кодировать большие размеры последовательностей. Приемлемые значения — 16 информационных слов и 2 проверочных. Исходя из этих упрощений, можно самому подобрать коэффициенты умножения кодировщика. Использование полинома Лагранжа или деления определенных многочленов в классическом РС продиктовано упрощением алгоритма поиска ошибки, для моей задачи этого не требуется. Поэтому единственным требованием остается — ранг матрицы коэффициентов с добавлением единичных и нулевых коэффициентов проверочных слов должен равняться 2. Тогда любые 2 потерянных слова можно восстановить, так любая подсистема размера 2 будет иметь однозначное решение. Интересно отметить, что это же требование автоматически приводит и к возможности найти и исправить одну ошибку в последовательности. Если попытаться исправить правильно полученное слово, то не получится одновременно выполнить 2 проверочных условия. Но в случае РС перебор линейно зависит от количества проверочных слов, а для общего случая произвольных коэффициентов — полный перебор всевозможных вариантов ошибки.

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

$\begin{pmatrix} 2 & 3 & 1 & 1 & 2 & 1 & 3 & 1 & 5 & 1 & 8 & 1 & 4 & 1 & 9 & 2 \\ 5 & 2 & 1 & 2 & 1 & 3 & 1 & 5 & 1 & 8 & 1 & 4 & 1 & 9 & 1 & 3 \end{pmatrix}$

Не сложно проверить, что ранг соответствующей расширенной матрицы

$\begin{pmatrix} 2 & 3 & 1 & 1 & 2 & 1 & 3 & 1 & 5 & 1 & 8 & 1 & 4 & 1 & 9 & 2 & 1 & 0\\ 5 & 2 & 1 & 2 & 1 & 3 & 1 & 5 & 1 & 8 & 1 & 4 & 1 & 9 & 1 & 3 & 0 & 1 \end{pmatrix}$

равен 2. В отличии от РС эти коэффициенты я подбирал руками, что для 16 слов данных и 2 проверочных сделать не сложно. Используя идею алгоритма Горнера, можно обойтись лишь умножением на 2 и сложением. Первое проверочное слово вычисляется как:

$\begin{aligned} d_{1}+d_{2}+d_{3}+d_{5}+d_{6}+d_{7}+d_{8}+d_{9}+d_{11}+d_{13}+d_{14} \\ + 2\cdot\bigl(d_{0}+d_{1}+d_{4}+d_{6}+d_{15}+2\cdot(d_{8}+d_{12}+2\cdot(d_{10}+d_{14}))\bigr) \end{aligned}$

а второе

$\begin{aligned} d_{0}+d_{2}+d_{4}+d_{5}+d_{6}+d_{7}+d_{8}+d_{10}+d_{12}+d_{13}+d_{14}+d_{15} \\ + 2\cdot\bigl(d_{1}+d_{3}+d_{5}+d_{15}+2\cdot(d_{0}+d_{7}+d_{11}+2\cdot(d_{9}+d_{13}))\bigr) \end{aligned}$

Сложение здесь выполняется как операция XOR для обычных чисел, а умножение на 2 — отдельная функция, которая должна быть векторизована.

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

if(loss_num==2 || loss_num==3 || loss_num==5 || loss_num==7 || loss_num==9 || loss_num==11 || loss_num==13){    word[loss_num]=subtracted_code_word1;}else if(loss_num==15 || loss_num==1){    word[loss_num]=subtracted_code_word1^subtracted_code_word2;}else if(loss_num==0){    word[loss_num]=m2(subtracted_code_word1)^subtracted_code_word2;}else{    word[loss_num]=subtracted_code_word2;}

где, m2() — функция умножения на 2.

Некоторые детали реализации и сравнение скорости выполнения

Для работы с векторами на esp32-s3 я использовал simd библиотеку. Некоторые особенности этой библиотеки на данный момент (либо я что-то не так понял) — поддержка лишь знаковых типов, отсутствие поддержки операций сдвигов, попытке компиляции без assert и проверок приводит к ошибкам компиляции. Но тем не менее ее достаточно для реализации умножения на 2. Базовая скалярная процедура умножения на 2, на которой я основывался, выглядит так:

uint8_t m2(uint8_t x){    return (x << 1) ^ ((x >> 7) * 0x1D);}

Сдвиги можно получить соответствующим применением функции vec_mul_scalar() из библиотеки «vector_basic_functions.h». Итоговый векторный вариант функции умножения на 2 выглядит следующим образом:

void m2(vector_t *v){    vec_mul_scalar(v, 1, &tmpv, 7);    vec_mul_scalar(&tmpv, -29, &tmpv, 0);    vec_mul_scalar(v, 2, v, 0);    vec_xor(v, &tmpv, v);}

tmpv — внешний буферный вектор подходящего размера для промежуточных вычислений. Хотя он делает эту функцию нереинтерабельной, но думаю это не критично. Ресурс векторных регистров весьма ограничен, и активное параллельное использование и контекстное сохранение может снизить скорость вычисления,за которое я боролся. Кодирование выполнялось в callback и о распараллеливанием я не собирался заниматься.

Остальные процедуры тривиальны. Для умножения двух чисел в поле Галуа я использовал таблицы логарифма и степеней 2. Любопытно сравнить скорости кодировщика при векторизации и без. Сравнение будем проводить для кодирования блока данных в 20 кбайт. Это размер блоков, которые мне надо было передавать.

С векторизованными операциями на кодирование тратится 176 мкс. Если использовать функцию с побайтным умножением на 2 , то время составляет 2867 мкс. Компилятор был настроен в режиме optimize for perfomance (-O2). Интересно отметить, что время работы выполнения побайтового кодирования зависит от последовательности выполнения операций. Если работать с массивами данных как с векторами, то время кодирования сокращается до 1881 мкс. Возможно компилятор частично сам векторизует операции. Ниже представлены обе побайтовые функции:

void CalcCode_2867us(uint8_t *TBUF, int size){    for(int i=0; i<size; i++){        TBUF[16*TSIZE+i]=TBUF[17*TSIZE+i]=0;        for(int j=0; j<sizeof(coeff11)/sizeof(coeff11[0]); j++) TBUF[16*TSIZE+i]^=TBUF[coeff11[j]*TSIZE+i];         TBUF[16*TSIZE+i]=m2s(TBUF[16*TSIZE+i]);        for(int j=0; j<sizeof(coeff12)/sizeof(coeff12[0]); j++) TBUF[16*TSIZE+i]^=TBUF[coeff12[j]*TSIZE+i];         TBUF[16*TSIZE+i]=m2s(TBUF[16*TSIZE+i]);        for(int j=0; j<sizeof(coeff13)/sizeof(coeff13[0]); j++) TBUF[16*TSIZE+i]^=TBUF[coeff13[j]*TSIZE+i];         TBUF[16*TSIZE+i]=m2s(TBUF[16*TSIZE+i]);        for(int j=0; j<sizeof(coeff14)/sizeof(coeff14[0]); j++)  TBUF[16*TSIZE+i]^=TBUF[coeff14[j]*TSIZE+i];                 for(int j=0; j<sizeof(coeff21)/sizeof(coeff21[0]); j++)  TBUF[17*TSIZE+i]^=TBUF[coeff21[j]*TSIZE+i];         TBUF[17*TSIZE+i]=m2s(TBUF[17*TSIZE+i]);        for(int j=0; j<sizeof(coeff22)/sizeof(coeff22[0]); j++)  TBUF[17*TSIZE+i]^=TBUF[coeff22[j]*TSIZE+i];         TBUF[17*TSIZE+i]=m2s(TBUF[17*TSIZE+i]);        for(int j=0; j<sizeof(coeff23)/sizeof(coeff23[0]); j++)  TBUF[17*TSIZE+i]^=TBUF[coeff23[j]*TSIZE+i];         TBUF[17*TSIZE+i]=m2s(TBUF[17*TSIZE+i]);        for(int j=0; j<sizeof(coeff24)/sizeof(coeff24[0]); j++)  TBUF[17*TSIZE+i]^=TBUF[coeff24[j]*TSIZE+i];     }}

void CalcCode_1881us(uint8_t *TBUF, int size){    for(int i=0; i<size; i++) TBUF[16*TSIZE+i]=TBUF[17*TSIZE+i]=0;            for(int j=0; j<sizeof(coeff11)/sizeof(coeff11[0]); j++) for(int i=0; i<size; i++) TBUF[16*TSIZE+i]^=TBUF[coeff11[j]*TSIZE+i];     for(int i=0; i<size; i++) TBUF[16*TSIZE+i]=m2s(TBUF[16*TSIZE+i]);    for(int j=0; j<sizeof(coeff12)/sizeof(coeff12[0]); j++) for(int i=0; i<size; i++) TBUF[16*TSIZE+i]^=TBUF[coeff12[j]*TSIZE+i];     for(int i=0; i<size; i++) TBUF[16*TSIZE+i]=m2s(TBUF[16*TSIZE+i]);    for(int j=0; j<sizeof(coeff13)/sizeof(coeff13[0]); j++) for(int i=0; i<size; i++) TBUF[16*TSIZE+i]^=TBUF[coeff13[j]*TSIZE+i];     for(int i=0; i<size; i++) TBUF[16*TSIZE+i]=m2s(TBUF[16*TSIZE+i]);    for(int j=0; j<sizeof(coeff14)/sizeof(coeff14[0]); j++) for(int i=0; i<size; i++) TBUF[16*TSIZE+i]^=TBUF[coeff14[j]*TSIZE+i];         for(int j=0; j<sizeof(coeff21)/sizeof(coeff21[0]); j++)  for(int i=0; i<size; i++) TBUF[17*TSIZE+i]^=TBUF[coeff21[j]*TSIZE+i];     for(int i=0; i<size; i++)TBUF[17*TSIZE+i]=m2s(TBUF[17*TSIZE+i]);    for(int j=0; j<sizeof(coeff22)/sizeof(coeff22[0]); j++)  for(int i=0; i<size; i++)TBUF[17*TSIZE+i]^=TBUF[coeff22[j]*TSIZE+i];     for(int i=0; i<size; i++)TBUF[17*TSIZE+i]=m2s(TBUF[17*TSIZE+i]);    for(int j=0; j<sizeof(coeff23)/sizeof(coeff23[0]); j++) for(int i=0; i<size; i++) TBUF[17*TSIZE+i]^=TBUF[coeff23[j]*TSIZE+i];     for(int i=0; i<size; i++)TBUF[17*TSIZE+i]=m2s(TBUF[17*TSIZE+i]);    for(int j=0; j<sizeof(coeff24)/sizeof(coeff24[0]); j++) for(int i=0; i<size; i++) TBUF[17*TSIZE+i]^=TBUF[coeff24[j]*TSIZE+i];     }

coeff.. — массивы значений, какие данные используются в схеме Горнера на соответствующем уровне умножения, m2s – байтная функция умножения на 2. Возможно, переписав побайтную процедуру в более грамотном виде, компилятор сможет еще лучше оптимизировать ее выполнение, но векторный подход понятнее и прямолинейнее.

Заключение

В этой короткой заметке была рассмотрена адаптация мощного универсального алгоритма коррекции ошибки РС под конкретную задачу и под конкретные аппаратные возможности контроллера. В результате был получен вполне ощутимый выигрыш во времени вычисления и достойная порция дофамина. Надеюсь, представленные здесь некоторые аспекты работы с esp32-s3 будут полезны программистам, стремящимся по максимум использовать аппаратные возможности железа.

ссылка на оригинал статьи https://habr.com/ru/articles/1033246/