Алгоритм блочного шифра "Кузнечик"

Занимался оптимизацией алгоритма шифрования Кузнечик [ГОСТ Р 34.12-2015]. Последние действия могу объяснить, а откуда получены промежуточные, забываю. Надо зафиксировать знание. Чтобы восстановить логику, пишу статью. Цель статьи -- показать, как можно оптимизировать алгоритмы в конечном поле.

Цикл алгоритма состоит из трех преобразований (LSX): наложение колюча(X), табличная подстановка(S) и некоторые линейные действия в конечном поле GF(2^8) с полиномом 0x1С3.

Преобразование L состоит из множества умножений и сложений в конечном поле. Операцию L можно исключить и сделать табличной. А наоборот, нельзя ли вычислить L честно и в тоже время быстро?

L состоит из 16 циклов R замешивания и умножения. В результате можно написать произведение матрицы коэффициентов (LM) 16х16 элементов на вектор из 16 элементов: LM*(a) -- линейная алгебра в поле GF. Если в нашем распоряжении есть векторная инструкция умножения в конечном поле, то задача может быть решена.

В архитектуре ARM Neon есть инструкция полиномиального векторного умножения (vpmull_p8), которая позволяет выполнить 8 умножений параллельно, результат - вектор 128 бит состоящий из 8 произведений выполненных без переносов(полиномиальных). Чтобы вернуть результат умножения в поле, надо выполнить операцию редуцирования. Операция матричного умножения в поле GF(2^8)

// Операция  матричного умножения в поле GF(2^8)
// редуцирование вынесено из цикла
#include <arm_neon.h>
poly16x8_t v0 = {0};
poly16x8_t v1 = {0};
int i;
for(i=0; i<16; i++) {
    // размножаем значение на все элементы
    poly8x8_t a8 = (poly8x8_t)vdup_n_u8(sbox[a[i]]);
    // загружаем столбцы матрицы коэффициентов
    poly8x16_t p8 = LMT[i];
    // умножаем младшую часть вектора poly8x8
    v0 ^= vmull_p8(a8, vget_low_p8 (p8));
    // умножаем старшую часть вектора poly8x8
    v1 ^= vmull_p8(a8, vget_high_p8(p8));
}

В результате вектор {v0,v1} содержит результат матричного умножения, 16 элементов по 16 бит в двух векторных регистрах, который далее следует редуцировать до 8 бит, чтобы вернуть в поле GF(2^8).

Редуцирование выполняется по методу Баретта. В алгоритме редуцирования используется обратное число - константа баретта (0x1B5) и сам полином (0x1С3).

Редуцирование по Баретту напоминает остаток от деления.
q = a/P = a*B
r = a - q*P

В нашем случае нужно заменить операции на соответствующие операции в поле GF

Алгоритм №__. Редуцирование по Баретту B=0x1B5, P=0x1C3.
На Входе: [a1:a0]
    [q1:q0] := a1 ⊙ B -- операция умножения
    [r1:r0] := [a1:a0] ⊕ q1 ⊙ P
На выходе: r0

// редуцирование вынесли из цикла
    poly8x8_t t0;
	poly16x8_t t;
	const poly8x8_t B5 = vdup_n_p8(0xB5);
	const poly8x8_t C3 = vdup_n_p8(0xC3);
    // сдвиг вправо элементов вектора 16 бит на 8 бит с заужением
    t0 = (poly8x8_t) vshrn_n_u16((uint16x8_t)v0, 8);
    t  = vmull_p8(t0, B5);// 
    t0^= (poly8x8_t) vshrn_n_u16((uint16x8_t)t, 8);
    v0^= vmull_p8(t0, C3);// v0 := v0 - t0*P

Операция редуцирования выполняется над каждым элементом вектора v0 и v1. Результат - вектор poly16x8_t. Операция выполняется над вектором v1 - старшая часть вектора 16х16. Финально необходимо преобразовать вектор обратно к разрядности poly8x16_t

Аналогичный алгоритм может быть построен на операции умножения без переносов на платформе Intel, с использованием инструкции PCLMUL

//Матричное умножение в поле GF
    for(i=0;i<16;i+=1){
        a16[0] = sbox[a[i]];
        poly64x2_t L = (poly64x2_t)UNPACKLBW128((int8x16_t)LMT[i], Z);
        poly64x2_t H = (poly64x2_t)UNPACKHBW128((int8x16_t)LMT[i], Z);
        v0 ^= CL_MUL128(a16, L, 0x00);
        v1 ^= CL_MUL128(a16, L, 0x10);
        v2 ^= CL_MUL128(a16, H, 0x00);
        v3 ^= CL_MUL128(a16, H, 0x10);
    }

Привожу табличный алгоритм, чтобы было с чем сравнивать

//Векторное воплощение операции LSX
uint8x16_t LSX(uint8x16_t a, uint8x16_t K)
{
    poly64x2_t r0={0};
    poly64x2_t r1={0};
    int i;
    poly64x2_t bb= (poly64x2_t)(a ^ K);
    uint64_t b0 = bb[0];
    uint64_t b1 = bb[1];
#pragma GCC unroll 8
    for(i=0;i<8;i++) {
        uint64_t aa = a[i];
        r0 ^= LS[i  ][(b0>>(i*8))&0xFF];
        r1 ^= LS[i+8][(b1>>(i*8))&0xFF];
    }
    return (uint8x16_t)(r0 ^ r1);
}

Результат практический: удалось получить производительность на Intel Core (10th Gen) до 20 тактов на байт при шифровании MGM. В то время, как производительность табличного метода - в три раза выше. Теоретически можно пытаться увеличить производительность алгоритма за счет использования векторных инструкций VPCLMUL с разрядностью 256 и 512 бит. На платформе Arm Neon ожидаю сравнимую производительность с табличным методом, поскольку в цикле используется меньше инструкций. Данных по производительности на ARM пока нет.

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