Сборочное производство

Собрали линию на новой площадке. Принимаем заказы на автоматический монтаж.

В состав линии входят: трафаретный принтер, чип шутер с приставкой для паллет с процессорами, длинная печка на десять зон и станция контроля.

Хеш плата для Майнинга

Делаем хеш-плату для биткоинов.. с жидкостным охлаждением. Киловатт с квадратного децеметра.

Я поставил рекорд по плотности монтажа, все сигналы и сила разведены по одной стороне. Искусство схемотехники!

Плата должна выдерживать токи в 40-60А через кристалл размером 8х8мм. Дорожки усилены медной полосой. Питание радается каскадом: ток через один каскад чипов питает следующий каскад.

Размер платы всего 210x90мм. Успейте купить!

Цифровая подпись ГОСТ -- быстрое редуцирование по модулю простого числа

Настал момент поделиться ноу-хау, как ускорять криптографию ГОСТ. Вернее одна маленькая деталь, но самая важная. Как быстро выполнять редуцирование - модуль простого числа. Как вообще строится цифровая подпись. Есть математика на бублике, в поле целых чисел по модулю простого числа. На этом строится и RSA и ECDSA и ГОСТ. Можно сказать модуль числа - это самая критическая операция - остаток от деления.

Любая операция на бублике: это умножение большого числа (256 или 512 бит) или сложение большого числа сопровождается последующем редуцированием(взятием по модулю).
Простые числа в криптографии часто выбирваются исходя из простоты восприятия человеком и с возможностью оптимизации вычислений, например в ГОСТ представлены простые числа вида P = 2^(n-1)+p0 и P = 2^n - p0 где n=256 или n=512 бит.

Все операции в поле выполняются по модулю. Например, умножение по модулю - это умножение и модуль числа.

// Умножение с редуцированием по модулю простого числа.
void bn_mulm(bn* r, bn* a, bn* b)
{
	int cy = bn_mul(r, a, b);
	bn_mod(r, cy, Prime);
}

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

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

Что в таком случае редуцирование?

Представим результат операции - как перенос (сx- число, которое вылезло за пределы 512 бит) и число X, которое осталось в пределах 512 бит {cx, x}. При этом число x имеет разрядность 512 бита, а перенос - 32 или 64 бита, зависит от архитектуры процессора, ориентируемся на 64 бита.

Наша задача сводится к тому, чтобы вычесть из числа {cx,x} простое число несколько раз, чтобы перенос стал нулевым. При этом остаток помещается в 512 бит, но может требовать дальнейших операций по редуцированию. Нам важно только одно - результат редуцирования помещается в заданные 512 бит. Операцию можно выполнить в два этапа.

Алгоритм 1. Быстрое редуцирование в поле P = 2^N - p0
1. {cy,x} : = {cx, x} - (2^512-p0)*cx;
В результате cy может быть равен нулю или единице.
2. if cy>0 then {cy,x} : = {cy, x} - P;

Чтобы это работало на P надо наложить некоторые ограничения p0 < 2^(n-m), где m - разрядность переноса.

Алгоритм можно упростить:

Алгоритм 1.1. Быстрое редуцирование P = 2^N - p0
1. {cy,x} : = x + p0*cx;
2. if cy>0 then x : = x + p0;

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

Теперь рассмотрим второй случай, когда простое число представленно суммой P=2^(N-1) + p0, где N -- разрядность поля, 256 или 512.

Алгоритм 2. Быстрое редуцирование P = 2^(N-1) + p0
1. {cy,x} : = {cx, x} - 2*(2^511+p0)*cx = x - 2*p0;
В результате cy может быть равен нулю или (-1). 
2. if cy<0 then {cy,x} : = {cy, x} + P;

При значении cy == (-1) x находится в интервале [-2*p0; -1], в старшем бите единица (1).

Полагаю равносильно заменить P на втором шаге на 2P. Cуть операции редуцирования не меняется.

Алгоритм 2.1. Быстрое редуцирование P = 2^(N-1) + p0
1. {cy,x} : = {cx, x} - 2*(2^511+p0)*cx = x - 2*p0*cx;
В результате cy может быть равен нулю или (-1).
2. if cy<0 then {cy,x} : = {cy, x} + 2*P = x + 2*p0;

Отдельно рассматриваем редуцирование при сложении и вычитании

Алгоритм 3. Сложение с редуцированием
1. {cy, x} := a+b
В результате cy принимает значения 0 или 1.
2. if (cy>0) {cy,x} := {cy,x} - (2^N - p0) = x + p0
Может понадобится третий шаг, очень мало вероятно:
3. if (cy>0) {cy,x} := {cy,x} - (2^N - p0) = x + p0

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

// Алгоритм 4. Умножение с накоплением (со сложением)
uint64_t bn_mla_512(uint64_t* r, uint64_t *a, uint64_t d)
{
    unsigned __int128 ac = 0;
    int i;
    for (i=0;i<8; i++) {
        ac +=  (unsigned __int128)d*a[i] + (unsigned __int128)r[i];
        r[i] = (uint64_t)ac;
        ac = ac>>64;
    }
    return (uint64_t)ac;
}

Операция "Умножение с накоплением" используется в реализации Алгоритма 1.1, быстрого редуцирования.

// Алгоритм 5. Умножение с вычитанием
 int64_t bn_mls_512(uint64_t* r, uint64_t *a, uint64_t d)
{
    __int128 ac = 0;
    int i;
    for (i=0; i<8; i++) {
        ac += (unsigned __int128)r[i] - (unsigned __int128)a[i]*d;
        r[i] = (uint64_t)ac;
        ac = ac>>64;
    }
    return (int64_t)ac;
}

Операция "Умножение с вычитанием" используется в реализации Алгоритма 2.1, быстрого редуцирования.

Эффективность описания операции (Alg.5) сомнительна, в рабочем проекте я использую ассемблер. Данное описание мне понадобилось для отладки под GCC 10 на платформе Intel x86_64.

Сжатие методом Хаффмана (Huffman Codes)

На картинке изображен метод кодирования Хаффмана, предложенный еще в 1952 году, 70 лет назад. Глядя на эту схему, приходит в голову сортированный список или стек(приоритеная очередь) с сортировкой по весовому коэффциенту.
Алгоритм построения кодов состоит из двух этапов. На первом этапе мы составляем и сортируем очередь указанным на картинке способом. На каждом шаге создаем композитный элемент очереди, включающий ссылки на два последних элемента и тем самым формируем из списка дерево -- еще одна динамическая структура. На втором этапе совершаем обход дерева, степень вложенности элементов в дерево будет давать длину кода. Данный метод используется для вычисления длин кодов, одновременно с эти можно вычислять сами коды, приняв направление ветвления, право - лево за ноль или единицу. Расчет кодов предлагается выполнять методом обхода дерева, а это -- рекурсия.

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

No |wght|bits|code
 0 | 20 |  2 |00
 1 | 18 |  3 |010
 2 | 10 |  3 |011
 3 | 10 |  3 |100
 4 | 10 |  3 |101
 5 |  6 |  4 |1100
 6 |  6 |  5 |11010
 7 |  4 |  5 |11011
 8 |  4 |  5 |11100
 9 |  4 |  5 |11101
10 |  4 |  5 |11110
11 |  3 |  6 |111110
12 |  1 |  6 |111111

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

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

Как избавиться от приоритетной очереди. От нее не избавиться, но мы могли бы использовать сортировку массива или создавать динамические данные из массива. Сортировка массива дает слишком много циклов и много операций копирования. В общем случае мы используем объект типа список из массива элементов и методы: "вставить с приоритетом"(list_insert_sorted), "получить следующий", и "выделить новый элемент из массива" (list_new). Используем процедуру обхода списка, в одном направлении -- односвязный список.

Это мы сейчас обсуждаем генерацию кодов, которая востребована на этапе упаковки данных, сжатия.

Декомпрессия выполнятся иначе. Для декомпрессии на вход подается таблица битовых длинн кодов, только третья колонка в нашем примере, вторая и четвертая колонки таблицы: частота использования и коды -- не заполнены.

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

Метод декодирования с использованием дерева ветвления - тоже сравнительно медленный. Чтобы понять насколько медленно скажем так: символы основного алфавита кодируются длинами от 1 до 15 бит. Если разброс длин таблицы кодирования большой, то может понадобиться 15 циклов, чтобы найти нужный код. Однако, тут надо иметь ввиду, что длинные коды встречаются реже, чем короткие. Так что, берем оценку по среднему - восемь циклов на декодирование каждого символа.

Метод декодирования с использованием кодовых таблиц в теории будет работать быстрее, потому что разбор ведется сразу по несколько бит, в среднем по 8 бит для литералов и длин, и по пять бит для смещения. Т.е. ожидаем производительность табличных методов в несколько раз выше. Однако, для декодирования кодов с максимальной длиной 15 бит нужны таблицы декодирования размером в 32-128 кБайт, а их генерировать возможно дольше, чем разбирать поток по одному биту. По этой причине мы рассматриваем методы, которые используют частичное декодирование с помощью таблиц, а затем вычисляют остаток.

Ставим задачу: получить производительность декодирования - один цикл на символ, в среднем. Как это сделать. Для начала надо научиться по префиксу кода, находить индекс в таблице. Готовим коды с выравниванием по левому краю, как указано в таблице, чтобы можно было выполнять сравнение чисел.

 0 | 20 |  2 |00XXX
 1 | 18 |  3 |010XX
 2 | 10 |  3 |011XX
 3 | 10 |  3 |100XX
 4 | 10 |  3 |101XX
 5 |  6 |  4 |1100X
 6 |  6 |  5 |11010
 . . .

В таком виде можно вычислять индекс по кодам. Привожу алгоритм описанный в обзорной статье [ACM Comput. Surv., Vol. 1, No. 1, Article 1. Publication date: June 2019].

    1) Найти такое ℓ чтобы first_code_l[ℓ] ≤ code < first_code_l[ℓ + 1];
-- Коды упорядочены по убыванию частоты использования, найти какому диапазону принадлежит код.
    2) пусть offset ← (code − first_code_l[ℓ]) >> (L − ℓ);
-- находим смещение в таблице - номер символа. 
    3) Находим сам символ s ← first_symbol[ℓ] + offset;

Для поиска в таком режиме применяем следующий цикл, именно он будет определять производительность метода:

// Алгоритм поиска кода по таблице.
    set ℓ ← search_start[code >> (L − t)];
    while code ≥ first_code_l[ℓ + 1] do
        set ℓ ← ℓ + 1.

Я бы сказал так, выражаясь шаблонами: мы используем хеш таблицу для поиска кода.

Сначала из таблицы некоторой фиксированной длины мы выбираем величину минимальной длины кода (t). Затем подбираем код по таблице. Казалось бы выигрыша особенного нет, потому что получаем цикл от минимальной до максимальной длины кода.

Для отладки и исследования алгоритма я ввел еще одну таблицу - search_end - максимальная длина кода с заданным префиксом, которая ограничивает глубину поиска кода.

// Предложенный алгоритм поиска
    set ℓ ← search_start[code >> (L − t)];
    set u ← search_end  [code >> (L − t)];
    while (ℓ < u) && code ≥ first_code[ℓ + 1] do
        set ℓ ← ℓ + 1.

Тут надо принять во внимание, что значения в таблице search_strart (минимальная длина кода) и search_end (максимальная длина кода) для данного префикса чаще всего совпадает, и мы статистически за одну проверку условия (ℓ < u) получаем битовую длину и сам код, без множества итераций. Поиск по таблице можно было бы оптимизировать с использованием бинарного поиска, но это не сработает, потому что коды статистически так распределены, что их проверять лучше последовательно.

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

И тут пришла в голову мысль, а не посмотреть ли как другие решают эту задачу. Посмотрел. Опять нашел код с исходниками Марка Адлера и отказом от права контролировать копии, с 1996 года в составе GNU, gzip -- этот код есть на каждом сервере. Код тридцатилетней давности, 30 лет, в коде упомянаются проблемы компиляции под MSDOS и OS/2, с ипсользованием Turbo C. В файле ТОДО рекомендуется использовать оптимизации CRC от Росса Вильямса, который уже пятнадцать лет, как на пенсии.

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

// Алгоритм декодирования
    do{
        t ← t.data + (взять_из_потока(t.n_bits) & маска_выборки(t.extra)) 
    } while нужно_продолжить(t.extra);
    сдвинуть_поток(t.n_bits);
    выполнить_обработку_данных(t.data, t.extra).

Примерно такой шаблон кода должен решать задачу декодирования. Более того, если этот фрагмент выполняется в цикле, таблицами (t) и переходами внутри таблицы (t.data) можно организовать контектстно-зависимый разбор. Маска выборки может быть линейная или не линейная. Допустим нам нужно декодировать программу из системы команд одного процессора в систему команд другого процессора. Или скажем, нужно написать итератор по системе команд Intel x86, т.е. по префиксам выделять длину инструкции -- тот же алгоритм декодирования, тот же шаблон. Или понадобилось разбирать синтаксис и грамматику в интерпретаторе команд, - тоже самое.

Конкретно, нам нужно выбрать инструкцию, на которой будет все работать. Для пионеров можем предложить реализацию на инструкции Intel AVX512_BITALG:

Synopsis
   __mmask16 _mm_bitshuffle_epi64_mask (__m128i b, __m128i c)

Instruction: vpshufbitqmb k, xmm, xmm
CPUID Flags: AVX512_BITALG + AVX512VL

Т.е на одной такой инструкции можно декодировать входной поток с использованием перестановки бит. Комбинация бит во входном потоке приводит к переходу по таблице декодирования. Результат - маска.

Мы выбираем простой способ формирования масок ((1<<(n))-1), который можно реализовать в одну инструкцию Intel BMI или BMI2;

В зависимости от целевой платформы, мы можем с таким алгоритмом достичь производительности в один такт(цикл) на кодовое слово. И тут происходит основное расстройство - мы живем в золотом веке программирования - этот подход уже реализован в gzip. С той лишь разницей, что реализован он средствами 30 летней давности, в инструкциях ANSI С. Более современная интерпретация алгоритма пристуствует в проекте zlib-ng.

Это сильно расстраивает, но наша цель -- предложить самый быстрый и компактный алгоритм декодирования.

Я потратил выходные на то чтобы научиться составлять таблицы декодирования и подобрать функцию декодера. И вот что получилось.

// Алгоритм декодирования
#define BIT_MASK(n) ((1<<(n))-1)

    t = t + (stream & BIT_MASK(N));
    if (t->extra) {
        t = te + t->n + ((stream>>N) & BIT_MASK(t->extra));
    }
    stream >>= t->bits;
    return t->n;

Функция использует две таблицы: (t) - основаная с разрядностью индекса N и дополнительная (te) - для декодирования кодов с длиной больше N. Декодирование выполняется не более чем в два этапа. На выходе алгоритм выдает декодированное значение. Параметр N выбирается разный для таблиц декодирования литералов (N=8) и для декодирования дистанций (N=5). При данных значениях параметра вторая таблица используется редко. Операции извлечения бит из потока могут выполняться с использованием инструкций Intel BEXTR (BMI1). Полагаю данная реализация является оптимальной с точки зрения производительности.Если при этом удасться сохранить лаконичность кода, то и самая быстрая в своем роде (быстрее gzip, zlib и zlib-ng).

Суммараная производительность алгоритма будет складываться из времени подготовки таблиц, времени работы алгоритма декодирования (inflate), и времени подгрузки данных.