MGM -- (Multilinear Galois Mode) векторная реализация

Вышли новые рекомендации к Российской национальной криптографии, режим шифрования MGM, который обещают использовать как основной и пожалуй единственный, если не считать рекомендации к протоколу CRISP. Режим MGM использует полиномиальное умножение в поле Галуа, по полиному x¹²⁸ + x⁷ + x² + x¹ + 1

Опишу порядок создания реализации алгоритма MGM. До начала работы надо иметь работающую реализацию алгоритма зашифровывания "Кузнечик" и тестовые вектора от MGM. Тестовые вектора нашел в draft стандарте IETF [Multilinear Galois Mode (MGM)]

Нужны глубокие познания в области умножения полиномов. Мы хотим использовать векторные инструкции полиномиального умножения, умножения без переносов.
На архитектуре Intel обещают операцию pclmul -- умножение без переноса двух чисел 64 бит, результат. Эта инструкция присутствует практически во всех новых процессорах Intel. И ожидаем появления векторной инструкции pclmul, которая позволяет производить одновременно 4 умножения 64бит на векторе 512 бит (AVX512+VPCLMULQDQ). Эта инструкция будет доступна на серверах Ice Lake в следующем году (2020). Но, алгоритмы на ее основе можно отладить уже сегодня с использованием эмулятора платформы Intel SDE. Я скачал SDE и убедился, что на нем можно эмулировать конкретные семейства процессоров Intel. Я использую компилятор GCC 9.2 под MSYS64, который хорошо осведомлен о новых процессорных инструкциях Intel и поддерживает целевую платформу Ice Lake.

Я не думаю, что векторные инструкции сильно ускорят процесс, потому что Intel обещает только один или два блока AVX512, в то время как каждое ядро способно выполнять четыре инструкции одновременно. Есть подозрение, что четыре инструкции 128 бит могут выполняться с такой же скоростью, что и одна инструкция AVX512 -- это пессимистический взгляд на вещи. Я не знаю, как устроен процессор на самом деле, но есть подозрение, что блоки AVX умеют объединяться в два блока AVX2 или в один блок AVX512. Возможно что в новых процессорах AVX2 -- неделимый блок, два блока AVX2 умеют объединяться в один блок AVX512. Если бы я разрабатывал процессоры Intel, я бы так и сделал, для экономии ресурсов. Почитал сообщения на форуме Intel, подтверждают, что быстродействие векторных инструкций AVX512 в два раза ниже чем быстродействие AVX2 на процессорах Xeon Scalable. В спецификации на процессор указывается число блоков AVX512. Подтверждают, что один блок AVX512 имеет такую же производительность на такт, как две инструкции AVX2. Но тактовая частота блоков AVX512 снижается от загрузки, она чуть ниже чем у ядра. Поэтому производительность получается даже ниже, чем на  инструкциях AVX2. Исходя из этих утверждений я хочу разработать алгоритм под векторные инструкции AVX2+VPCLMULQDQ, два умножения на инструкцию.


Начнем разработку с умножения без переносов.

Для начала решаем еще более простую задачу.
Пишем алгоритм умножения без переносов 64 бит на 64бит. с последующим редуцированием 128 бит в 64 бит. Редуцирование выполняется по полиному x⁶⁴ + x⁴ + x³ + x¹ + 1

uint64_t gfmul64(uint64_t a, uint64_t b)
{
 poly64x2_t r = CL_MUL128((poly64x2_t){a}, (poly64x2_t){b}, 0x00);
 uint64_t x1 = r[1];
 uint64_t x0 = r[0];
 x1 = x1 ^ x1>>63 ^ x1>>61 ^ x1>>60;
 return x0 ^ x1 ^ x1<<1 ^ x1<<3 ^ x1<<4;
}

Второй вариант алгоритма пишем с редуцированием на операциях умножения без переносов.

uint64_t gfmul64(uint64_t a, uint64_t b)
{
 poly64x2_t r = CL_MUL128((poly64x2_t){a}, (poly64x2_t){b}, 0x00);
// Редукция Барретт'а:
 poly64x2_t t = CL_MUL128(r, (poly64x2_t){0x1BULL}, 0x01) ^ r;
 r ^= CL_MUL128(t, (poly64x2_t){0x1BULL}, 0x01);
 return r[0];
}

В этом варианте я использую две константы 0x1B для редуцирования по методу Барретта. Константы высчитывал честно по алгоритму, который представлен в предыдущем сообщении. На самом деле первый вариант умножения gfmul64 тоже основан на редуцировании по методу Барретта только умножение разложено на операции сдвига и сложения в поле.
Проверяю лаконичность высказываний компилятора:

        vpclmulqdq  $0,  %xmm1, %xmm0, %xmm1
        vpclmulqdq  $17, %xmm2, %xmm1, %xmm0
        vpxor       %xmm1, %xmm0, %xmm0
        vpclmulqdq  $17, %xmm2, %xmm0, %xmm0
        vpxor       %xmm1, %xmm0, %xmm0

Я боюсь ошибиться. Поэтому сделал еще и третий вариант побитовое умножение.

uint64_t gfmul64_2(uint64_t a, uint64_t b)
{
    uint64_t r = 0;
    int i;
    for (i=0; i< 64; i++){
        if (b & (1ULL<<i)){
            r ^= a;
        }
        a = (a<<1) ^ (((int64_t)a>>63) & 0x1BULL);
    }
    return r;
}

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

Из популярной литературы я знаю алгоритм редуцирования на основе сдвигов [Carry-Less Multiplication and Its Usage for Computing the GCM Mode, алгоритм 4].

Теперь опишем умножение без переносов 128 бит на 128 бит. с последующим редуцированием 256 бит в 128 бит. Редуцирование выполняется по полиному x¹²⁸ + x⁷ + x² + x¹ + 1

Начнем с референсной реализации, образцовой, на которую ссылается Intel, Sun OpenSolaris, кто с кого списывал:

/// (Copyright 2009 Sun Microsystems, Inc.)
void gfmul(__m128i x_in_m, __m128i y_m, uint64_t *res)
{
   uint64_t R = { 0xe100000000000000ULL };
   struct aes_block z = { 0, 0 };
   struct aes_block v;
   uint64_t x;
   int i, j;
   uint64_t *x_in=(uint64_t*)&x_in_m;
   uint64_t *y=(uint64_t*)&y_m;
   v.a=y[1];
   v.b=y[0];
   for (j = 1; j>=0; j--) {
      x = x_in[j];
      for (i = 0; i < 64; i++, x <<= 1) {
         if (x & 0x8000000000000000ULL) {
            z.a ^= v.a;
            z.b ^= v.b;
         }
         if (v.b & 1ULL) {
             v.b = (v.a << 63)|(v.b >> 1);
             v.a = (v.a >> 1) ^ R;
         } else {
             v.b = (v.a << 63)|(v.b >> 1);
             v.a =  v.a >> 1;
         }
      }
   }
   res[0] = z.b;
   res[1] = z.a;
}

Этот алгоритм так себе описан. Он медленный, 128 циклов. Вроде люди пытались описать с использованием векторных типов и псевдо-функций векторных инструкций (wmmintrin), но что-то пошло не так. Я даже знаю, что пошло не так, компилятор Sun не имел никаких векторных расширений и это был черти какой год, GCC тогда тоже не умел их использовать, а Intel умел только через явное использование этих самых псевдо-функций. Пробую описать тоже самое, теми средствами, которые понятны современному компилятору GCC 9.+.
Надо заметить, что в GCM используется обратный порядок следования бит. А нам нужен прямой порядок.

Алгоритм 1 полиномиальное умножение 128x128 бит
Шаг 1: рассчитываем умножение без переноса для двух чисел: 
  A = [A1:A0]; B= [B1:B0];
Покомпонентный результат умножения будет:
  A0•B0 = [C1:C0], A1•B1 = [D1:D0], A0•B1 = [E1:E0], A1•B0 = [F1:F0]
Шаг 2: составляем результат 256-бит:
 [A1:A0]•[B1:B0] = [D1:F1⊕E1⊕D0:F0⊕E0⊕C1:C0]

На шаге 2 можно применить способ умножения Карацубы

Алгоритм 2 полиномиальное умножение 128x128 бит методом Карацубы
Шаг 1: Покомпонентный результат умножения будет:
  A0•B0 = [C1:C0], A1•B1 = [D1:D0], (A0⊕A1)•(B0⊕B1) = [E1:E0]
Шаг 2: составляем результат 256-бит:
 [A1:A0]•[B1:B0] = [D1:D0⊕D1⊕E1⊕C1:D0⊕E0⊕C0⊕C1:C0]

Если считать, что инструкция умножения без переносов выполняется за такт, за то же время, что и операция Исключающее ИЛИ, то экономии не будет. Это можно проверить.
Описываем умножение по Алгоритму 1:

/// Алгоритм 1 полиномиального умножения 128x128
poly64x4_t CL_MUL256(poly64x2_t p, poly64x2_t v)
{
    poly64x2_t q  = CL_MUL128(p, v, 0x10) 
                  ^ CL_MUL128(p, v, 0x01);// Карацуба отдыхает
    poly64x2_t r0 = CL_MUL128(p, v, 0x00) ^ SLL128(q, 64);
    poly64x2_t r1 = CL_MUL128(p, v, 0x11) ^ SRL128(q, 64);
    return (poly64x4_t){r0[0], r0[1], r1[0], r1[1]};
}

Тут я использую встроенные псевдо-функции, которые призваны генерировать векторные инструкции.

/// SRL128 - сдвиг числа 128бит вправо на указанное число бит
static inline poly64x2_t SRL128(v2du a, const int bits) {
    return (poly64x2_t)__builtin_ia32_psrldqi128((v2di)a, bits);
}

/// SRL128 - сдвиг числа 128бит влево на указанное число бит
static inline poly64x2_t SLL128(poly64x2_t a, const int bits) {
    return (poly64x2_t)__builtin_ia32_pslldqi128((v2di)a, bits);
}

/// CL_MUL128 - полиномиальное умножение 64х64 бит
static inline poly64x2_t CL_MUL128(poly64x2_t a, poly64x2_t b, const int c) {
    return (poly64x2_t)__builtin_ia32_pclmulqdq128 ((v2di)a,(v2di)b,c);
}

Из сатьи взял алгоритм со сдвигами для редуцирования по полиному x¹²⁸ + x⁷ + x² + x + 1

Алгоритм 4 редуцирование по полиному x128
Обозначим входные данные [X3:X2:X1:X0] с разбивкой числа 256 бит по 64 бит.
Шаг 1: Сдвигаем X3 на 63, 62 и 57-бит вправо.
  A = X3>>63, B = X3>>62, С = X3>>57
Шаг 2: D = X2 ⊕ A ⊕ B ⊕ C
Шаг 3: Выполняем свдиги над числам 128 бит
 [E1:E0] = [X3:D]<<7
 [F1:F0] = [X3:D]<<2
 [G1:G0] = [X3:D]<<1
Шаг 4: Сложение (XOR) над числми [E1:E0], [F1:F0], [G1:G0] и [X3:D]. 
 [H1:H0] = [X3 ⊕ E1 ⊕ F1 ⊕ G1 : D ⊕ E0 ⊕ F0 ⊕ G0]
Шаг 5: Возвращаем результат [X1 ⊕ H1: X0 ⊕ H0].

Пояснения к алгоритму 4. Мы можем сделать тоже самое с использованием умножения. Константа Барретта при использовании редуцирования будет равна самому полиному, сдвиги {0,1,2,7} соответствуют умножению на число 0x87, хвост от полинома (x⁷ + x² + x + 1).

Алгоритм 5 редуцирование по полиному x128
Обозначим входные данные [X3:X2:X1:X0] с разбивкой числа 256 бит на 64 бит.
Шаг 1: Умножаем старшую часть числа X на константу Барретта (0x87).
  [A3:A2:A1:A0] = [X3:X2]•0x87
Шаг 2: B = A ⊕ X
Шаг 3: Умножаем старшую часть числа B на полином (0x87)
  [C3:C2:C1:C0] = [B3:B2]•0x87
Шаг 4: D = C ⊕ X
Шаг 5: Возвращаем результат [D1:D0].

Благодаря равенству констант алгоритм можно упростить, сократив число умножений с четырех до трех.

Алгоритм 6 редуцирование по полиному x128
Обозначим входные данные [X3:X2:X1:X0] с разбивкой числа 256 бит на 64 бит.
Шаг 1: Умножаем старшую часть числа X на константу Барретта (0x87).
  [0:A2:A1:A0] = [X3:X2]•0x87
Шаг 2: Умножаем старшую часть числа B на полином (0x87)
  [0: 0: 0:C0] = [0 :A2]•0x87
Шаг 3: D = C ⊕ A ⊕ X
Шаг 4: Возвращаем результат [D1:D0].

Совместим алгоритм №6 с алгоритмом №1 умножения полиномов.

Алгоритм 7 редуцирование по полиному x128 + 0x87
Шаг 1:
  X2•0x87 = [A1:A0], X3•0x87 = [B1:B0]
Шаг 2:
  B1•0x87 = [C1:C0]
Шаг 3: Возвращаем [B0⊕A1⊕С1⊕X1:A0⊕C0⊕X0]

Можно заметить, что заменой умножения на сдвиги можно получить алгоритм №4.

// Алгоритм №7: Редуцирование по полиному x128 + 0x87
poly64x2_t gf128_reduction(poly64x2_t r0, poly64x2_t r1)
{
 const poly64x2_t Px ={0x87};
 poly64x2_t b  = CL_MUL128(r1,Px, 0x01);
 poly64x2_t a  = CL_MUL128(r1,Px, 0x00);
 poly64x2_t c  = CL_MUL128( b,Px, 0x01);
 return r0 ^ a ^ c ^ SLL128(b,64);
}

Проверяем результат компиляции:

// $ gcc -O3 -march=native -o - -S test.c
        vpclmulqdq      $1, %xmm2, %xmm0, %xmm3
        vpclmulqdq      $0, %xmm2, %xmm0, %xmm0
        vpclmulqdq      $1, %xmm2, %xmm3, %xmm2
        vpxor   %xmm0, %xmm2, %xmm2
        vpslldq $8, %xmm3, %xmm3
        vpxor   %xmm3, %xmm1, %xmm3
        vpxor   %xmm3, %xmm2, %xmm0

-- ничего лишнего.
Проверяем тестовые вектора, взял из статьи.

Умножение в конечном поле GF(2128), P(x) = x128 + x7 + x2 + x + 1
a = 0x7b5b54657374566563746f725d53475d
b = 0x48692853686179295b477565726f6e5d
gfmul128 (a, b) 
  = 0x040229a09a5ed12e7e4e10da323506d2

Результат. Для алгоритма MGM разработал эффективную реализацию умножения в конечном поле GF(2⁶⁴) и GF(2¹²⁸). Умножение состоит из двух частей - умножение без переноса и редуцирование. В цикле достаточно выполнять только умножение без переноса, а этап редуцирования можно вынести из цикла.

Несколько комментариев к статье

Операция PCLMUL -- тормозная. Она долго грузится, но исполняется за такт. На разных платформах она выполняется по-разному. Хотелось бы сюда вставить табличку, из которой видно насколько эффективно использовать эту инструкцию. Быстродействие складывается из двух параметров: Задержка и Исполнение. Инструкции могут долго загружаться, но выполняться одна за другой. Например, для новых серверов я ожидаю быстродействие 6 тактов задержка на исполнение и 1 такт на исполнение. Это значит, что группа инструкций PCLMUL не связанных между собой могут грузиться 6 тактов, и исполняться по 1-ому такту на команду. В некоторых случаях процессор выполняет инструкции в несколько потоков, тога пиковая производительность получается выше. Следить за тем чтобы команды эффективно грузились и переставлять порядок исполнения - это дело компилятора и процессора. НО зачастую уменьшив взаимосвязь между параметрами алгоритма можно добиться большей производительности.

Для понимания, когда стиот применять инструкции, а когда это совсем не эффективно, привожу таблицу по производительности операции PCLMUL на разных платформах. Таблицу нашел в итнернетах, источник - [Agner Fog. Lists of instruction latencies, throughputs..] Таблицы таймингов (задержки и производительность инструкций) для некторых моделей процессоров Intel представлены в
[#]Intel(c) Intel® 64 and IA-32 Architectures Optimization Reference Manual, приложение С.3 [Appendix C.3, “Latency and Throughput”].

Процессор      Задержка+обработка 
-- AMD -- 
Bulldozer     12 7
Piledriver    12 7
Steamroller   11 7
Excavator      5 5
Ryzen          4 2
Jaguar         3 1
-- Intel --
Silvermont    10 10 -- ужасно медленно!!
Nehalem       12 8
Ivy Bridge    14 8
Goldmont       6 4
Haswell        7 2
KnightsLanding 6 2
Skylake        7 1
SkylakeX       7 1
Broadwell      5 1

Можно утверждать, что в современных процессорах инструкция выполняется достаточно быстро, чтобы вытеснить другие решения. Все познается в сравнении, другие решения могут сроится на таблицах, операциях PXOR и PSHUFB. Например, на платформе Skylake операция PSHUFB выполняется за 1+1 такт, операция PXOR за 1+0.5 такта.

[22.12.2019]

Нашел возможность дальнейшей оптимизации Алгоритма редуцирования (№7), сократил число умножений до двух.

Алгоритм №8 редуцирование по полиному x128 + 0x87
Шаг 1:
  X3•0x87 = [B1:B0]
Шаг 2:
  (B1⊕X2)•0x87 = [D1:D0]
Шаг 3: Возвращаем [B0⊕D1⊕X1:D0⊕X0]

// Алгоритм №8: Редуцирование по полиному x128 + 0x87
poly64x2_t gf128_reduction(poly64x2_t r0, poly64x2_t r1)
{
 const poly64x2_t Px ={0x87};
 poly64x2_t b = CL_MUL128(r1,Px, 0x01) ^ SLL128(r1,64);
 poly64x2_t d = CL_MUL128(b, Px, 0x01) ^ SLL128(b, 64);
 return r0 ^ d;
}

Возвращаясь к режиму GCM, приведу мой вариант редуцирования для отраженного порядка следования бит.

// Алгоритм №9: Редуцирование по полиному x128 + 0x87
poly64x2_t gf128r_reduction(poly64x2_t r0, poly64x2_t r1)
{
 const poly64x2_t Px ={0x87};
 poly64x2_t b = CL_MUL128(r1,Px, 0x01) ^ SLL128(r1,64);
 poly64x2_t d = CL_MUL128(b, Px, 0x01) ^ SLL128(b, 64);
 return r0 ^ d;
}

Алгоритм №9 я сам вывел. Однако, до меня в примерно в такой форме его использовали в OpenSSL () -- этот факт обнаружился уже позже. Алгоритм в OpenSSL сложно распознать, поскольку он дан уже в ассемблерном коде. Я его распознал там по сигнатуре, схожая последовательность инструкций PCLMUL+PALIGNER, как в моем случае.

[23.12.2019]

Сравнил быстродействие варианта умножения методом Карацубы по Алгоритму №2 (три инструкции PCLMUL) и Алгоритма №1. Существенного выиграша не обнаружил. Думаю следует использовать вариант без Карацубы, по Алгоритму №1, на новых процессорах этот вариант будет быстрее.

Блочный шифр ГОСТ "Кузнечик" -- оптимизация

Речь идет об оптимизации алгоритма блочного шифра ГОСТ "Кузнечик". Я давно отложил задачу на сладкое, а тут время появилось ее доделать.
В основе шифра Кузнечик лежат две операции. 1) нелинейное биективное преобразование - табличная подстановка. 2) Что-то вполне линейное, в смысле алгебра линейная в конечном поле, на основе арифметики Галуа, с редуцированием по полиному 0x1С3.

Я увидел возможность преобразовать алгоритм, сделать его быстрее.
Итак структура алгоритма такова [ГОСТ 34.12-2015]:
E[K] = X[K₁₀]LSX[K₉]... LSX[K₂]LSX[K₁](a)
Я пропущу описание того, что сделано в первой серии. Я умею выражаться длинными векторами, для этого использую векторные расширения языка Си. Блочный шифр -- это преобразование вектора состоящего из 16 элементов разрядностью 8 бит. Для этого использую описание типа:

typedef uint8_t uint8x16_t __attribute__((__vector_size__(16)));

Алгоритм защифровывания выглядит так

uint8x16_t kuzn_encrypt(KuznCtx* ctx, const uint8x16_t a)
{
    uint8x16_t S = a ^ ctx->K[0];
    int i;
    for (i=0; i<9; i++){
        S = LS(S) ^ ctx->K[i+1];
    }
    return S;
}

Шифрование -- это вот такой алгоритм, внутри которого есть функция LS, всего девять циклов по развернутому ключу {K₁,K₂, ...K₁₀}.
Зашифровывание -- это последовательное применение преобразований:
X[K₁], S, L, X[K₂], S, L, ... X[K₉], S, L, X[K₁₀]
Преобразование X[K](a) = K ⊕ a
Преобразование S -- это таблица подстановок -- "нелинейное биективное преобразование".

uint8x16_t LS(uint8x16_t a)
{
    int i;
    for(i=0; i< 16;i++) a[i] = sbox[a[i]]; // -- подстановка, 
    uint8x16_t v = R16(a); // -- то самое преобразование L
    return v;
}

Вот добрались до сути. КАК оптимизировать преобразование R¹⁶?
Преобразование R выполняется над векторами 16 элементов по 8 бит и содержит 16 циклов.
Кроме всяких пермутаций это преобразование содержит много операций умножения в конечном поле GF(2) с полиномом p(x) = x⁸ + x⁷ + x⁶ + x + 1, или Px = 0x1C3.
Путем долгих мучительных преобразований можно записать, что каждый элемент входного вектора зависит от каждого элемента выходного вектора. Всего 16 наборов констант для каждого элемента. v_i = SUM(a_j * LMT_ij) . LMT - это матрица коэффициентов, которую по ходу пришлось пересчитывать и транспонировать, чтобы преобразование стало выглядеть линейно.

uint8x16_t GMULv16(uint8x16_t L, uint8x16_t a);
uint8x16_t R16(uint8x16_t a)
{
    int i;
    uint8x16_t v = {0};
    for(i=0; i< 16;i++) {
        v ^= GMULv16(LMT[i], a[i]); // -- векторное умножение в поле Px = 0x1C3
    }
    return v;
}

Операйця GMULv16 выполняет умножение каждого элемента вектора с редуцированием по полиному 0x1C3.
На этом этапе я остановился в прошлый раз. Большим достижением казалость то, что алгоритм удалось сделать без ветвления и с ипользованием векторных инструкций. GMULv16 можно заменить на 16 таблиц по 256 значений по 16 байт каждое, 64кБайт. Но это не наш путь.

В цикле можно использовать умножение без переноса, а финальное редуцирование выполнять за пределами цикла.
v_i = SUM((a_j * LMT_ij) mod P(x)) =  SUM(a_j * LMT_ij) mod P(x)

В результате получаем такой алгоритм:

uint8x16_t R16(uint8x16_t a)
{
    int i;
    uint16x16_t vh = {0};// -- 16 элементов разрядностью 16 бит. 256 бит.
    for(i=0; i< 16;i++) {
        vh ^= CL_MULv16(LMT[i], a[i]); // -- векторное умножение без переносов
    }
    uint8x16_t v = RR (vh, Px);// -- финальное редуцирование по полиному P(x),
    return v;
}

Размер таблицы - 16 элементов по 16 байт, 256 байт. В цикле получаем умножение 128 бит на 8 бит. Это умножение можно представить 4я инструкциями умножения полиномов в поле разрядностью 64 бит.

Может быть есть компромис между размером таблицы и безумным умножением?
Я так понимаю, что большие таблицы медленно обрабатываются, потому что в кеш память не влезают. Если мы делаем алгоритм для контроллера, то таблицы 65кБайт не лезут во флеш память, не лезут в память оперативную, сильно увеличивают объем кода. Обработка таблицы будет идти со скоростью работы памяти, а не со скоростью ядра процессора. Для встроенных приложений нужен компактный код. На нашей целевой архитектуре есть инструкция полиномиального умножения (ARM Cortex-A8+ NEON, ARMv8, Intel Core+AVX), но в контроллерах ARMv7e-m ничего такого нет, вынужденно используем таблицы. Нам понадобятся оба варианта алгоритма. Таблица не лезит в память, поэтому я сделал две таблицы, одна - это таблица умножения на числа 0,1,2,3...15, младшие биты числа. Вторая тоже таблица умножения на числа 0,16,32,64...256, - умножение на старшие биты числа.

uint8x16_t R16(uint8x16_t a)
{
    int i; 
    uint8x16_t v = {0}; 
    uint8x16_t a0 = a & 0xF,  a1 = (a >>4) & 0xF; 
    for(i=0; i< 16;i++) { 
        v ^= GMULv16_L[i][a0[i]] ^ GMULv16_H[i][a1[i]]; 
    }
    return v; 
}

В этом варианте умножения используются 32 таблицы по 16 значений по 16 байт, 8кБайт. Надо проверить, какой из двух вариантов быстрее: с подстановкой 8кБайт или с безумным умножением...
Проверил. Сильно зависит от процессора. На 3-м поколении Intel Core-i7 умножение без переносов работает крайне медленно, в пять раз медленнее чем табличный вариант. На 7-м поколении Core-i7 умножение работает в два раза медленнее чем табличный вариант. Разница есть. Табличный вариант позволяет получить скорость зашифровывания около 8 Мбайт/сек. Это все равно медленно.

[11.12.2019] Дальше я исследую возможность переноса алгоритма с умножением полиномов на векторные инструкции ARM NEON. У меня под рукой есть кросс-компилятор GNU ARM Embedded Toolchain, GCC. И есть под рукой плата BeagleBone, на которой могу запустить полученный код. Платформа содержит процессор ARM Cortex-A8 с архитектурой ARMv7-A и поддерживает инструкции NEON.
В системе команд NEON имеется инструкция полиномиального умножения VMULL.P8, которая оперирует с векторами poly8x8_t -- 8 элементов по 8бит. Инструкция производит результат - вектор 128 бит, 8 элементов по 16 бит. Подробнее см. [NEON™ Version: 1.0 Programmer’s Guide] на сайте ARM.

Мой код выглядит так:

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

В коде использую псевдофункции, которые описаны в заголовке и в результате генерируют инструкции NEON. Исследую возможные комбинации ключей компилятора, чтобы задействовать инструкции NEON:

$ arm-eabi-gcc -print-multi-lib

Нахожу подходящую комбинацию ключей:

$ arm-eabi-gcc -march=armv7-a+fp -mcpu=cortex-a8 -mfpu=neon-vfpv4 \
  -mfloat-abi=hard -O3 -S -o - kuzn.c

Содержимое цикла преобразцется в набор инструкций

        vdup.8   d18, r3        -- размножить значение 8x8
        vld1.64 {d24-d25}, [r1:64]!  -- загрузить LMT 8x8x2
        vmull.p8 q13, d18, d24  -- умножить вектора 8x8 
        vmull.p8 q9,  d18, d25  -- умножить вектора 8x8
        veor     q8,  q8,  q13  -- исключающее ИЛИ 16x8
        veor     q11, q11, q9   -- исключающее ИЛИ 16x8

Ничего лишнего, одно умножение - одна инструкция процессора. Все инструкции используют вектора 8x16, умножение проивзодится над векторами 8x8, результат умножения имеет размерность 16x8. Не путать - умножение проиводится без переносов.
Можно ли считать этот результат хорошим. Не знаю. Объективно надо замерить скорость и сравнить результат. Код получился компактный, но не факт что инструкция умножения будет выполняться быстрее, чем обращение к памяти. Делаю 1Миллион циклов зашифровывания и замеряю время -- получаю результат = 2.5 секунды, на платформе beaglebone Texas Instrumens Sitara AM335x. Это получается скорость шифрования 6.4Мбайта/сек. Замеряю результат табличного варианта алгоритма зашифровывания = 4.2 секунды. Мы выиграли время и место. Дальше за счет разворачивания цила удалось повысить скорость зашифровывания при использовании умножения до 7.4Мбайт/сек. Это ХОРОШИЙ результат!

Резюме 2019/2020

Ищу работу!

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

Люблю программировать. Умею писать операционные системы, умею писать криптографию, умею на время создавать безумное количество электроники для автоматизации, чтобы сразу работало.
Умею вести производство военной техники, умею разрабатывать военную технику с гарантией 7 лет и сроком службы до 40 лет.
Умею заниматься снабжением военных и гражданских проектов практически любых. Умею работать со спец. счетами и гос контрактами. Умею отбиваться от налоговых и прокурорских проверок, хотя это просто везение, научился выживать в условиях травли. Умею сохранять данные в облаках и на выделенных серваках, шифровать и нарезать мелко-мелко большие данные. Умею жонглировать серверами, когда их выносит ФСБ или МВД. Умею писать обработку протоколов ФТС. Умею писать блок чейны и архивы. Умею разбирать архивы документов XML со скоростью ... с любой скоростью. Имею опыт разработки корпоративных информационных систем на этих технологиях.

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

Почему-то не умею извлечь прибыль из этого.
Скучаю, ищу работу...Безумно устал от всей суеты. В отпуск хочу, в увольнительную. Хочу программы писать со скуки и воспитывать программистов. Достаточно одного.

Резюме 2011

Ищу работу!
Чего умею?

• Умею быть преподавателем. Лет 5 отработал преподавателем электроники. Создал свой собственный курс по микросхемотехнике в Политехническом университете.
• Могу управлять небольшим коллективом разработчиков.
• Умею проектировать электронику: измерительную, цифровую и силовую. 
• Умею писать встроенное программное обеспечение для ARM, 8051.
• Умею производить спроектированную электронику. Умею организовать контрактное производство электроники штучное, мелко серийное и серийное.
• Умею создавать сети для организаций, на 100-1000 рыл. Вообще-то умею создавать виртуальные транс-национальные корпорации не слезая со стула, вернее проектировать и настраивать аппаратуру и ПО для них. Умею в одиночку обслуживать 40 Кисок (CISCO) и десяток серверов. Хорошо знаком с GNU/Linux, FreeBSD, Oracle Solaris. Могу дописывать и исправлять код каждой из них, если понадобится.
• Знаю толк в информационной безопасности. Могу настроить инфраструктуру организации, включая сервисы e-mail, web, vpn, ip телефонию и многое другое. 
• Программировать люблю, на С. Умею проектировать все, что касается ядра операционной системы. Могу написать свою операционную систему с нуля. Так и сделал. Файловую систему для SSD. Люблю описывать сетевые протоколы, люблю писать интерпретаторы разных языков. Умею проектировать алгоритмы управления роботами. Умею расчитывать траектории движения без использования вещественных чисел с компенсацие ошибок округления. Умею писать реализации криптографических алгоритмов. Умею все это реализовать в аппаратуре программируеых матриц. Могу процессор спроектировать, правда не за один месяц.
• Могу писать пркиладное программное обеспечение, имею богатый опыт разработки ПО для систем реального времени.
• Могу эксперименты ставить по физике полупроводников, знаю толк в террагерцовой оптоэлектронике, правда лет 10 этим не занимался.

Почему-то не умею извлечь прибыль из этого, странно.
Скучаю, ищу работу...

Упражнения с BACnet - журнал операций

Я когда начинал описывать протокол BACnet, сделал такой сетевой "драйвер", который загружает формат Wireshark PCAP и подает на вход службы разбора пакетов APDU. Это позволило написать и отладить разбор форматов ASN кодирования BACnet.
С тех пор прошло много времени. Много кода... Сейчас Мне понадобилось сохранять конфигурацию контроллера в ответ на команду Re-инициализации. Вообще, как себе представляешь можно работать в контроллере с базой данных объектов? И при этом требуется сохранять конфигурацию между сессиями, при перезагрузке. Я придумал, что надо писать во флеш журнал операций модификации данных. Просто так, линейно вперед, пока место не закончится. Когда контроллер просыпается после перезагрузки он должен накатить журнал операций. Журнал, чтобы не плодить разные сущности я решил паковать прямо командами протокола, транзакциям. Например, транзакция понятная - создать объект, дописать свойства к объекту. По сути мы можем сохранять все операции записи типа WriteProperty, прямо из буфера обмена без изменения писать в журнал - во флеш. Вот такая замечательная идея.
Как отладить замечательную идею?
Журнал во флеш, флеш в контроллере. Контроллер перегружается... Не удобно.
Кроме контроллера есть сервер сбора данных, роутер BACnet. Предлагаю сохранить его конфигурацию в файл и поизучать каким-нибудь инструментом. Нужна отладка для журнала файловой системы.
Берем приложение, пригодное для отладки, для рассматривания журналов. Да вот оно - Wireshark. В итоге я пишу функцию, функционал, для записи данных не во флеш, а в файл со структурой PCAP и изучаю, что получилось в итоге. А в итоге имею последовательность команд протокола CreateObject, CreateObject, CreateObject... и на этом выводе примерно за день отлаживаю функцию экспорта объектов из базы данных в журнал.

База данных - это структура типа дерево объектов.
Foreach (DeviceInfo->ObjectTree, (Callback) save_object_cb, to_file);
Т.е. описываем некоторую процедуру обхода дерева с записью каждого элемента в файл/журнал. Так выглядит сериализация дерева.
Каждая запись состоит из шапки::{длина данных, тип записи}, пакета данных, контрольной_суммы. После чего можно записать следующую запись. В шапку можно записать и тип записи, в нашей концепции - это тип сервиса - CreateObject, DeleteObject, WritePropertyMultiple.
По сути нам не нужен этот идентификатор сервиса, потому что CreateObject без параметров можно расценивать как Delete, А операцию повторного создания объекта расценивать, как WritePropertyMultiple. Писать в журнал поток команд протокола подкупает своей простотой.
При ре-инициализации устройства через перезагрузку можно переписать журнал заново. При загрузке просто накатить этот журнал.

Есть проблема: как минимизировать объем данных в журнале, как делать инкрементную запись -- писать только те параметры, которые изменились с прошлой версии базы.

[07.01.2020]

Написал вывод в журнал, с учетом файлов. Получилось две операции: CreateObject, AtomicWriteFile.
Нашел свой старый проект, журналируемая файловая система для флеш. Тогда мы использовали четыре вида полей: APPEND - дописывание файлов, ATTRIB - дописывание/изменение атрибутов, CREATE - создание, DELETE - удаление, SKIP - технологическая операция - заполнение пробелов, для выравнивания на блок записи.
Сейчас думаю использовать старые идеи только расширить понятие файловой системы до объектов. CREATE -- создание объектов CreateObject, ATTRIB - запись атрибутов в формате WritePropertyMultiple, APPEND - дописывание файлов в форме AtomicWriteFile, DELETE - пометить на удаление из базы DeleteObject.

ARMv7m -- операции BFI и BFX и битовые строки

В приложении к криптографии я задался вопросом как "правильно" писать исходный код на языке Си, чтобы компилятор использовал инструкции BFI (Bit Field Insert) и BFX(Bit Field Extract).

Мой алгоритм берет 4 бита и использует их в качестве индекса массива. В примере я описываю "нелинейное биективное преобразование" по таблице подстановок (ГОСТ 34.12-2015 п.5.1.1).
Для отладки алгоритма я использую компиляцию в ассемблер и рассматриваю, какими командами компилятор выражается.
> gcc -march=native -O3 -o - -S magma.c
-- на экран выводится ассемблерный код оптимизированный под мой процессор
Или под целевую платформу ARMv7e-m
$ arm-none-eabi-gcc  -mthumb -mcpu=cortex-m4 -march=armv7e-m -mfloat-abi=hard -mfpu=fpv4-sp-d16 -o - -O3 magma.c -S

Исходно пишу на Си, рассматриваю результат.
Пишу определение, которое соответствует операции BFX

#define BEXTR(x, n, len) (((x) >> (n)) & ((1 << (len))-1))

В исходнике пишу так:

   for (i=0; i<8; i++){
       s = sbox[i][BEXTR(a,(i*4),4)];
       r |= (s & 0xF) <<(i*4);
   }

sbox - это таблица подстановок.
В ассемблерном коде возникает команда
ubfx...
Но мне никак не удается подобрать обратную операцию.
Пишу определение
#define BFI(x, y, n, len) x = ((x) & ~(((1 << (len))-1)<<(n))) | ((y & ((1 << (len))-1))<<(n))
Но в результате компилятор НЕ использует инструкцию bfi, определение НЕ работает.
Смотрю в документацию (Arm C Language Extensions, ACLE Q2 2019) нахожу пояснение, инструкция BFI описывается средствами языка Си, т.е. должно быть соответствующее ей выражение.
Через некоторое время, дошло, что в расширении языка С бывают свои битовые строки. Такое вот описание битовых полей (см код ниже) позволило однозначно задействовать инструкцию извлечения битовой строки и последующее помещение битовой строки обратно в регистр.

uint32_t T(uint32_t a)
{
    register union {
      struct {
        uint32_t u0:4;
        uint32_t u1:4;
        uint32_t u2:4;
        uint32_t u3:4;
        uint32_t u4:4;
        uint32_t u5:4;
        uint32_t u6:4;
        uint32_t u7:4;
      };
      uint32_t x;
    } r;
    r.x  = a;
    r.u0 = sbox[0][r.u0];
    r.u1 = sbox[1][r.u1];
    r.u2 = sbox[2][r.u2];
    r.u3 = sbox[3][r.u3];
    r.u4 = sbox[4][r.u4];
    r.u5 = sbox[5][r.u5];
    r.u6 = sbox[6][r.u6];
    r.u7 = sbox[7][r.u7];
    return r.x;
}

Это и есть биективное преобразование согласно ГОСТ, именно в таком виде оно попало в мою реализацию. -- Магия!!

В результате компилятор создает такой код:

   ubfx r2, r0, #16, #4   -- загрузить битовое поле
   add  r2, r2, r3
   ldrb r2, [r2,#64] @ zero_extendqisi2 -- загрузить байт из таблицы
   bfi  r0, r2, #16, #4   -- вставить битовое поле

CRISP 1.0 -- Протокол защищенного обмена для индустриальных систем

Протокол не важен. Важно что это первый нормальный набор векторов для Magma (ГОСТ 34.12-2015)в режиме генерации имитовставки (IMIT) и в режим гаммирования (CTR) для данных не выровненных. Сами режимы описаны в документе (ГОСТ 34.13-2015).

Протокол защищенного обмена выпущен ТС26 (Технический комитет по стандартизации российской криптографии) в виде методических рекомендаций.
Основная проблема в реализации режимов блочного шифрования - Это, не глядя в чужой код, получить точное соответствие с векторами. Сразу никогда не получается.
Где-то это мои проблемы - читать не умею. Где-то проблемы разработчика - описано плохо. Но чаще всего - это неточности в описании алгоритма, который надо как-то отладить, и ошибок в реализации может быть больше одной на пути рабочей версии.
Почему путь тернист. Берем сам алгоритм блочного шифра, "Магма" (ГОСТ 34.12-2015): описан, есть примеры - тестовые вектора. Взял алгоритм, сделал реализацию, но - не работает. Нигде в стандарте не сказано, в каком представлении даны числа в тестовых векторах. Порядок следования байт может быть Little-Endian или Network (Big-Endian). С ключами шифрования - тоже самое, они могут быть в нормальном порядке, задом наперед по 64 бита, задом наперед/совсем задом наперед - все число от старшего бита к младшему вывернуто. В процессе отладки я подбираю, как правильно записать, в какой последовательности, входные и выходные данные, чтобы числа сошлись. Сходятся.
Потом беру описание режима блочного шифра (ГОСТ 34.13-2015). В тестовых векторах та же путаница не ясно, как должны сходится эти вектора с изменение порядка следования байт или без изменения. Примеры почему то даны только для случая выровненных данных. В случае с усечением или без выравнивания на 64 бита, алгоритм остается не отлаженным.


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

Например, смотрим определения, чем отличается "бинарное представление", "байтовое представление" и числовое. Как разгадать этот ребус. По-русски ведь написано.
x||y -- конкатенация двоичных строк x и y ... при которой левая часть - слева, а правая - справа.
Потом пишут 1||0 - вот это чего за число. или это 0x02? Или это число 0x01? Оказывается это число 0x80.Оказывается оно таким после серии проб и ошибок.
binary() -- это представление символьной строки в виде байтовой строки. -- отпад.
byte() -- это представление числа в виде байтовой строки, при котором соответствующая итоговой байтовой строке двоичная строка .... -- вынос мозга.
вот у меня есть число 0xBAD как его представить в виде байтовой строки: "0B AD" или "AD 0B". Отсюда и берутся такие неоднозначности. Числа можно представить в "сетевом" порядке или в порядке свойственном аппаратуре LE|BE.
На самом деле, если четко договориться по терминам, что такое "двоичная строка", что такое "байтовая строка", что такое "число", что такое "символьная строка", как одно в другое переводится, то проблем нет. Авторы пытаются прикрыть неточность выражения мыслей обилием математических выражений типа Zn||...||Z2||Z1||Z0 - такие выражения никак не облегчаются понимание вопроса по выравниванию данных. Биты которые однозначно слева в "байтовом представлении" могут быть как слева так и справа. Но после прочтения не всегда это ясно.

Вот например каким количеством вариантов можно описать действие
IV = LSB_32(byte(SeqNum, 6)) -- по сути это означает что мы должны взять младшие 32 бита от байтового представления числа, при генерации байтового представления нужно взять 6 байт от числа, т.е 48 бит. Пишу варианты:
uint64_t SeqNum;
uint32_t IV = htonll(SeqNum); (1)
uint32_t IV = htonll(SeqNum<<16); (2)
uint32_t IV = (SeqNum); (3)
uint32_t IV = htonll(SeqNum)>>32; (3)
...
Или вот другая формула SN = 0^5||MSB_35(byte(SeqNum,6))
Если под функцией byte() автор понимает преставление байт из числового представления в сетевое, то это запишется так:
uint64_t SN = (swap64(SeqNum<<16)) & ((1ULL<<35)-1)); (3)
uint64_t SN = ((SeqNum>>(48-35)) & ((1ULL<<35)-1)); (4)
-- я реально перебираю все эти варианты, ориентируясь всего на один признак, результат равен или нет тестовому вектору. После подбора я нахожу, что автор имел ввиду:
SN = старшие_биты(0^5)||число(MSB_35(число(SeqNum,6)). -- сшивание битовых строк/чисел.
При укладке данных SN в "байтовую строку" происходит преобразование "в сетевом" порядке следования байт.
Строка = ...|| byte(SN,5) ||... -- число SN (40бит, 5 байт) разворачивается в порядке Big-Endian, в сетевом представлении. Откуда это следует? -- Я догадался.
Можно вчитываться в документацию, пытаться понять, все ли точно выражено в описании или все таки бинарное представление итоговой байтовой строки от числового представления сдвинутого на пять бит делается как-то иначе. А если в изложении несколько таких неточностей, или я не так понял и в моей реализации неточность или ошибка.
Может автор прав, есть однозначная запись битовых строк, операция MSB применима к битовым строкам, операция byte() преобразует число в битовую строку. Проверяем определения... Я понял, все дело в том, что битовые строки почему-то неявным образом переводятся в байтовые с выравниванием по старшему биту.
При использовании блочного шифра, я почему-то упускаю: надо данные из сетевого представления переводить в числовое перед использованием функции шифрования или надо результат функции шифрования переводить из локального в сетевое представление перед использованием. Такие вот проблемы. 

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

Криптография на эллипсе

Я много почитал статей по эллиптическим кривым, и теперь хочу описать свои впечатления.

Самое сильное из моих впечатлений - это часики. Это понятная аналогия, проясняет сознание.
Возьмем к примеру окружность, часики с круглым циферблатом, уравнение окружности:
x^2+y^2 = 1

Примем за начало отсчета точку O = (0,1) . 12часов 00 минут.
К чему мы клоним, чтобы сразу стало понятно. Мы вводим группу вращения стрелок на циферблате. Время складывается: минуты в часы...
Стараемся думать про одну стрелку. за минуту стрелка отклонилась на (x1, y1). Эту точку можно выразить через синусы и косинусы единичного перемещения 360/60 = 6 градусов. Координаты единичного перемещения обозначим точкой G=(x1,y1)

Утверждаем что это у нас Группа точек на кривой, Группа в математическом смысле.
Свойство Группы:
1) Существование нейтрального элемента, такого что P+O = P, O = (0,1)
2) Существование обратного элемента для каждого члена группы, -P= (-x,y). Перевели стрелки назад. P+(-P)=O
Вводим операцию удвоения точки, с ней можно будет ввести операцию умножения на скаляр через удвоение и сложение.
3) 2P = ... осторожно можно споткнуться... = (cos(2ф), sin(2ф)) = (x^2-y^2, yx+xy)
4) Закон сложения точек (x1, y1)+(x2, y2) = (cos(a+b), sin(a+b)) = (x1x2-y1y2, y1x2+x1y2)

А теперь можно заставить часики ходить...
Алгоритм №1 умножения на скаляр Q = kP, k-раз по минуте.

Q:=O;
for i=.. downto 0 begin
  Q := 2Q;
  if (k_i !=0) Q:= Q+P;
end

Минуты считаются по модулю 60. Число 60 не является простым, его можно на множители раскладывать. Число Р назовем генератором группы, обозначим буквой G чтобы всех запутать.

Алгоритм №2 умножение лесенкой Монтгомери.

Q:=O; P=G;
for i=.. downto 0 begin
  if (k_i !=0){
      Q:= Q+P, P=2P;
  } else {
      P: =P+Q, Q=2Q;
  }
end
return Q
Эти алгоритмы не зависят от того как выглядит операция удвоения и сложения. Алгоритмов умножения можно придумать великое множество: справа налево, слева направо, комбинированные, с окном, со сложением и вычитанием, с разложениями и окнами.
Оба алогоритма можно свести к одной или двум операциям: удвоение точки Q=2Q и Q=2Q+G
Или иными словами мы на каждом шаге алгоритма вычисляем либо удвоение Q_{2n} зная Q_{n}, или Q_{2n+1} зная Q_{n}, Q_{n+1} и Q_1


Я знаю сколько было времени на часах, когда я их запустил - это мой секрет, могу выразить его в минутах d (число минут). Могу рассказать всем, что если умножить Q = dG получится некоторая точка с координатами (Q.x, Q,y) - которая однозначно связана с моим секретом - это будет точка для проверки подписи. Я хочу подписать сообщение. Мне нужно представить сообщение в виде числа m. Тогда подписанное сообщение - это показание часиков R = (m*d)G. Которое можно проверить с использованием открытой точки: R = mQ.

Цифровая подпись, ее неподдельная сущность, держится на том, что никто не может вычислить обратное число d, зная R, m и Q. Или плохо старается.

Все известные алгоритмы нахождения обратного числа держаться на Алгоритм № 3 НОД наибольший общий делитель. На базе него можно получить алгоритм деления или нахождения обратного числа по отношению к операции умножения. Для изготовления понадобится число типа скаляр и операция над точками - уполовинивание. Уполовинивание связано с неопределенностью при операциях с нечетными числами, которую надо как-то разрешать на каждом шаге алгоритма.
...

И тут пришел Монтгомери со своими кривыми алгоритмами и решил все "упросить": проекция x в операции удвоения не зависит от координаты y!
2P = (2x^2-1, 2xy)
Это значит, что мы можем считать удвоение без использования второй координаты. После этого берем паузу и думаем, а как теперь считать сложение точек без использования y- координаты.
x = x_2 x_3 - y_2 y_3 =... надо выразить через X координаты точек P Q и G.
x = 2 x_2 x_3 - x_1
Утверждение такое:
x_{2n} = 2x_{n}^2 -1
x_{2n+1} = 2x_{n} x_{n+1} - x_1
Начальное состояние для вычисления умножения при n=0 (x_{n}, x_{n+1}) = (1, x1).

По сути венец творения Монтгомери - это утверждение, что операцию вычисления x координаты при сложении точек на эллиптической кривой, можно представить в общем виде, как
x_{m+n} = f(x_m, x_n, x_{m-n}) вот и думай теперь над своим алгоритмом.


Откуда взяты идеи с часами и Алгоритмы Монтгомери:
https://eprint.iacr.org/2017/293.pdf -- оттуда

Черный список

У нас есть телефонная станция на базе Asterisk расположена в датацентре. Всякий кому не лень пытается "взломать" пароль SIP и подключится на халяву к нашей телефонной сети. Меня это беспокоит.
Статистика неумолима. до 80млн. запросов в год на подбор пароля. Это невероятно много запросов, которые захламляют лог. Такого чтобы кто-то подключился и получил телефонных услуг на халяву, не замечено.

Виртуальная сеть для контроллера

Настраивал операционную систему. Выбирал таймеры для реализации виртуальной сети.
Нашел несколько полезных вещей, стал применять.

Таймеры на Windows. Необходимым оказалось подкручивать параметр мультимедиа таймера
Использую функцию timeBeginPeriod(wTimerRes);
Я заметил, что на нескольких разных компьютерах под управлением Windows 7 по разному ведет себя процесс отладки. Он то бежит как угорелый, символы быстро пробегают по экрану, то тормозит. Даже на одном компьютере я получил результат, что  после перезагрузки, он стал медленнее обрабатывать протокол. Замерил параметр задержки, оказалось, что вместо usleep(1000)=1мс, квант времени равняется 15мс. А у меня протокол должен работать с разрешением времени минимум 5мс. При 15мс, виртуальная сеть тормозит и вызывает таймауты в работе виртуальных устройств.

Таймеры pthread. Моя сеть должна работать на Windows, Linux и на моей операционке. Основа виртуальной сети - очередь таймерных объектов -- блоков памяти, которые доставляются в строго определенное время.
Исследовал разрешение таймеров.
clock_getres(CLOCK_MONOTONIC) возвращает 370 микросекунд. На разных процессорах эта цифра разная, но меньше 1мс.
clock_getres(CLOCK_REALTIME) возвращает 15.6мс.
Измерять время надо монотонным, иначе цифры округляются до безобразных величин.

Монотонный таймер оказался не очень то монотонным, на этом потерял целый день на отладку работы виртуальной сети и планировщика.
Чтобы сделать из монотонного таймера действительно монотонный применил такой ход:
(uint64_t)(tv_nsec + tv_sec*1000000000); Иногда в tv_nsec встречаются любые числа, неожиданные. Монотонным таймер становится только после такой операции.

 timestamp = osKernelSysTick();
 while ((uint32_t)(timestamp - tr->wait.timestamp) < tr->wait.timeout) {
     interval.tv_nsec = (tr->wait.timeout - (uint32_t)(timestamp - tr->wait.timestamp));
     clock_nanosleep(CLOCK_MONOTONIC, 0, &interval, &diff);
     timestamp = osKernelSysTick();
 }
Применил такой вариант ожидания. Ожидание применяется перед получением пакета данных.

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

void osAsyncQueuePut(osAsyncQueue_t* queue, void* data)
{
    List_t* tr = g_slice_alloc(sizeof(List_t));
    tr->data = data;
    volatile void** ptr = (volatile void**)&queue->tail;
    do {
        tr->next = atomic_pointer_get(ptr);
        atomic_mb();
    } while(!atomic_pointer_compare_and_exchange(ptr, tr->next, tr));
}
Эта операция добавляет элемент в список - вместо верхнего элемента. Список снимается в одно движение со стороны планировщика:
queue->head = atomic_pointer_exchange(&queue->tail, NULL);
Перед разбором список надо перевернуть, чтобы получить нормальный хронологический порядок.

Вот и все искусство.

Модель описания протокола, часть 2

Потратил около недели чтобы реализовать один конечный автомат из стандарта BACnet/MSTP. Хочу поделиться идеями реализации, пожаловаться на сложности. Во первых, надо сказать я его раза три написал. При том что описание получается от 300+ до 600+ строк. Это такой громадный switch() со встроенными проверками и switch() по командам протокола и ветвлением по разным признакам.

Модель описания протокола

Меня беспокоит проблема описания протоколов. Мне приходится сталкиваться с протоколами обмена для электронного документа оборота. Диаграмма состоит из десятка состояний и переходы между состояниями обусловлены получением документа определенного типа. Диаграмма не учитывает ошибки и обработку ошибок. С другой стороны я пытаюсь осилить стандарт BACnet, там все описание стандарта пропитано статусными моделями. Только переходы и условия весьма сложны, их не описать таблицей. Если добавлять в статусную модель обработку ошибок -- она перестает быть наглядной и пугает, такое описание только путает.

CMSIS RTOS osThread vs C11 threads.h vs pthread

Какой стандарт лучше выбрать. Выбор сделать не могу, однако хочу писать программки для контроллеров и для ПК на одном языке.

SVG: типовые ошибки векторной графики

Изучал исходники сairographics нашел интересные места.

Векторные операции и векторная графика

В этой статье суммирую опыт перевода графического приложения на векторные операции SSE/AVX. Приложение двумерное, плоская графика. Для обработки используется упакованный double. Цель статьи обобщить опыт работы с векторными типами данных приментельно к плоской графике и показать, как можно оптимизировать операции над векторами.

Некоторые операции не могут быть эффективно описаны на языке Си, привожу примеры описания базовых операций с использованием встроенных функций и векторных расширений языка.

Развитие концепции CMSIS RTOS

Зачем? Я пишу операционную систему. Надо. Для автоматики здания надо, для управления огромными боевыми человекоподобными роботами надо.

Чего не хватает? CMSIS RTOS -- хороше API для встроенных приложений. Но мне пришлось пересмотреть ее чуть ли не полностью, перетрясти. Я стараюсь ничего не менять. Но надо.

HTTP сервис, что скрывается в облаках

Как правильно сказать: "Добро пожаловать в облачные технологии"?! Наверное в современной терминологии Облаком называется граница знаний, куда вас не пустили, где технологии недоступные для вашего понимая, вам даже не рассказали что это.
Моя глобальная задача -- множество устройств, которые обладают коллективным разумом.  Частная задача -- эффективное взаимодействие множества серверов в пределах видимости, в одной серверной и в одной локальной сети, в демилитаризованной зоне, где важно быстро и четко.