GF2p8 Кодирование Рида-Соломона в изоморфном представлении поля

Коды Рида-Соломона считаются через экспоненты и логарифмы по таблице. А у нас с внедренеием инструкций Intel GFNI появилось две другие возможности: использовать инструкцию умножения и аффинные преобразования.

Для примера я взял коды РС из проекта библиотеки генеации QR-кодов. Там используется полином поля G(2^8) P = 0x11D. Кодирование производится в цикле. При переходе к операциям умножения в поле 0x11D, цикл кодрования выглядит так.

	for(i = 0; i < data_length; i++) {
		uint8_t fb =  data[i] ^ Z[0];
		for(j = 0; j < ecc_length-1; j++) {
			Z[j] = Z[j+1] ^ gf2p8_mul(fb, g[j], 0x11D);
		}
		Z[j] = gf2p8_mul(fb, g[j], 0x11D);
	}
	for(i = 0; i < ecc_length; i++) {
		ecc[i] = Z[i];
	}

Преобразование алгоритма для расчета в изоморфном представлении поля с образующим полиномом 0x11B:

	const uint64_t M = 0xFFAACC88F0A0C080;
	const uint64_t Mr= 0xFFAACC88F0A0C080;

	for(i = 0; i < data_length; i++) {
		uint8_t fb =  affine(M,data[i]) ^ Z[0];
		for(j = 0; j < ecc_length-1; j++) {
			Z[j] = Z[j+1] ^ gf2p8_mul(fb, g[j], 0x11B);
		}
		Z[j] = gf2p8_mul(fb, g[j], 0x11B);
	}
	for(i = 0; i < ecc_length; i++) {
		ecc[i] = affine(Mr, Z[i]);
	}

Данное преобразование предполагает так же изоморфные преобразования компонент вектора генератора кодов (g). Векторизация алгоритма может быть выполнена по вложенному циклу.

[1]The comeback of Reed Solomon codes N.Drucker, Sh.Gueron, V.Krasnov

Приложение

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

// Multiplication in Galua Fields
int gf_mul(int x, int y) {
    if (x == 0 || y == 0) return 0;
    return exp_[(log_[x] + log_[y])%255];
}
// Inversion in Galua Fields
int gf_inv(int x) {
    return exp_[255 - log_[x]];
}
// Division in Galua Fields
int gf_div(int x, int y) {
    if (x == 0) return 0;
    return exp_[(log_[x] + 255 - log_[y]) % 255];
}

Самая неудобная операция в этом случае - вычисление остатка от деления на 255. Избавился от деления.

// Остаток от деления на 255
uint8_t mod255(unsigned v) {
	uint32_t q = (v*0x1010102UL)>>(24);
	return q;
}

// Остаток от деления на 255
uint8_t mod255(unsigned v) {
	return v + (v>>8);
}

А этот вариант самый простой, применим для сложения логарифмов.

Такая вот оптимизация.

[2] Исходники на Github

Нахождение обратного числа

Одна из способностей, которая проявилась у меня с детства...

Пишу программы под всякую странную аппаратуру. В странной аппаратуре плохо считается операция деления. Например, берем самую совершенную, казалось бы, архитектуру Intel, в ней операция деления выполняется [AF1] очень медленно. Все арифметические инструкции, включая умножение, запускаются по 4 шт за такт, а вот деление только одно и продолжается целых 6 тактов.
Core i7 Coffee Lake: Умеет паралелить множество инструкций, а операция умножения делается только на одном из 8 конвейеров, имеет задержку 4 такта, на прохождение (Latency) и может подряд выполняться с производительностью 1 такт на умножение (Throughput). Я не беру векторные инструкции, беру простую задачу поделить одно число на другое. Деление вызывает настоящий ступор в системе, 26 тактов на загрузку инструкции (Latency) и 6 тактов зависание на конвейере.
Берем менее совершенную архитектуру, ARM Cortex, в документации указано число циклов на операцию умножения и деления, 2 и 12 соответственно. В 6 раз операция деления медленнее операции умножения. А еще есть контроллеры, которые вообще не имеют аппаратного деления. Эмуляция деления это сотни тактов процессора.

Вот, к примеру, мне понадобилось штамп времени разложить на число, месяц и год. Потом, часы, минуты и секунды. Начинаем считать: остаток от деления штампа времени на (24*60*60) - это секунды от начала дня. Потом делим на (60*60), потом на 60 ... сплошные операции деления на константу. Или берем задачу форматирования числа, надо на каждом десятичном разряде выполнить деление на 10 с остатком. Вот тут и приходит здравая мысль - это оптимизация, которая может выполняться на этапе компиляции. Целочисленное деление на константу можно заменить на умножение со сдвигом. Фактически мы заменяем деление на умножение на дробь. Не все числа "рациональные", не все обратные числа можно представить в виде дроби с достаточной точностью. Приближение первое: мы хотим представить константу в виде дроби 1/B ≈ C/2^n. Первое что пришло в голову n=32, брать только старшие 32 бита от целочисленного умножения с разрядность 32х32бит, но это так не работает. Надо умножать на самое большое число, для которого не возникает переполнение, использовать максимально возможное n, чтобы константа С годилась для всех целых чисел. Число С должно содержать 1(единицу) в старшем разряде, не меньше 0x80000000UL.
Оптимизация деления на константу на языке Си запишется так:
uint32_t Q = (A*C)>>n,
где С и n -- константы.
Для нахождения констант предлагаю алгоритм:

uint32_t div_c(uint32_t b) {// находим число С
    uint32_t C;
    int k = __builtin_ctz(b);// число ноликов сзади
    b>>=k;
    n = 63 - __builtin_clz(b);
    C = ((1ULL<<(n))/b)+1;// сдвигаем на число ноликов спереди
    n+=k;
    return C;
}

Проверить справедливость замены можно методом перебора чисел, это занимает не много времени.

uint32_t verify(uint32_t B, uint32_t C, int n) {
    uint32_t A;
    for (A=(~0UL); A!=0; --A){// все числа 32 бит
        uint32_t Q = ((uint64_t)A*C)>>32>>(n-32);
        if (Q != (A/B)) {
            break;
        }
    }
    return A;
}

Примеры таких констант Q=A/B = (A*C)>>n:
 

Q = (A* 0xCCCСCСCDULL)>>35; // Q=A/10 кратно /5 работает для всех A
Q = (A* 0xCCCСCСCDULL)>>34; // Q=A/5 работает для всех A
Q = (A* 0xAAAAAAABULL)>>33; // Q=A/3 работает для всех A

B=  3: 0xAAAAAAAB n=33
B=  5: 0xCCCСCСCD n=34
B=  7: 0x92492493 n=34 Q=A/7 для 31 бит.
B= 60: 0x88888889 n=37
B= 10: 0xCCCСCСCD n=35
B=100: 0xA3D70A3E n=38
B=1000: 0x83126E98 n=41
B=3600: 0x91A2B3C5 n=43
B=86400: 0xC22E4507 n=48
B=1000000: 0x8637BD06 n=51
B=365: 0xB38CF9B1 n=40 Q=A/365 для 31 бит.
B=153: 0xD62B80D7 n=39
B=255: 0x80808081 n=39

Некоторые константы нельзя представить в виде 32 битного умножения, точности не хватает. Для 31 битных целых чисел всегда хватает точности.

-- Магия чисел!

[AF1] https://agner.org/optimize/instruction_tables.pdf