четверг, 4 февраля 2021 г.

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

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

Пишу программы под всякую странную аппаратуру. В странной аппаратуре плохо считается операция деления. Например, берем самую совершенную, казалось бы, архитектуру 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

вторник, 1 сентября 2020 г.

Электронное правительство: настройка сервисов

 Мы по роду деятельности занялись "чисткой" в правительстве. Пишем запросы и анализируем ответы из правительства, потом пишем еще запросы и анализируем систему в целом. Систему в целом называем "электронное правительство" - подходящий термин.

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

Структура сети

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

mail.kodeks.vpn (Postfix) with ESMTP -- Постфикс мы уважаем, но по сути это сервер VPN на базе Linux, сервер осуществляет удаленный доступ. mail.kodeks.vpn (unknown [10.30.3.73]) Слово Кодекс в контексте не нравится. Проблема в том что ВСЕ правительство, все комитеты (в SPB) работает через один IP адрес [10.128.66.165] через один сервер [10.30.3.73]. Выход из строя или неправильная работа этого адреса - проблема для всей сети. Пользователи сети никак не отличаются, нельзя выделить источник по IP адресу. Аутентификация пользователей в Postfix по паролю или сертификату безопасности при отсылке сообщений по SMTP не производится. Сервер на нижнем уровне позволяет скрыть источник данных и не проверяет аутентичность сообщений. unknown [10.128.66.165], обратные адреса во внутренней сети не различаются. Соответственно, правила фильтрации трафика в Postfix не настроены.

В системе развернуты два сервиса проверки сообщений на вирусы

X-Antivirus: Dr.Web (R) for Unix mail servers drweb plugin ver.6.0.2.8
X-Antivirus-Code: 0x100000
X-KLMS-Rule-ID: 3
X-KLMS-Message-Action: skipped
X-KLMS-AntiSpam-Status: not scanned, whitelist
X-KLMS-AntiPhishing: not scanned, whitelist
X-KLMS-AntiVirus: Kaspersky Security 8.0 for Linux Mail Server, 
    version 8.0.1.721, not scanned, whitelist

Для проверки трафика на вирусность используется сервер
mail.smolny.uts (mail.cn.uts [10.20.20.3]) Все работает на Postfix один из сервисов работает, как Milter с пересылкой на localhost.
Сервер ничего не проверяет, вообще ничего, все письма проходят по белому списку, список настроен по IP адресу, а IP адрес у всех один и тот же.

mail.gov.spb.ru (mail.gov.spb.ru. [176.97.35.9]) -- это внешний интерфейс системы, через него осуществляется выход наружу из приватных сетей Правительства (СПб) в сеть интернет.

Результат проверки на вшивость исходящих сообщений:
Received-SPF: softfail (google.com: domain of transitioning noreply@letters.gov.spb.ru does not designate 176.97.35.9 as permitted sender) client-ip=176.97.35.9;
Authentication-Results: mx.google.com;
       spf=softfail (google.com: domain of transitioning noreply@letters.gov.spb.ru does not designate 176.97.35.9 as permitted sender) smtp.mailfrom=noreply@letters.gov.spb.ru

Проще говоря все иходящие сообщения классифицирются как убогая попытка навредить простому пользователю, потому что не проходят проверку. Гугл классифицирует их как "опасные", не проверенные, разве что в спам не перекладывает.

Фильтрация на основе SPF записей -- это старый и весьма простой механизм проверки от подмены источника сообщений.

В современных системах корпоративной почты используется механизм аутентификации серверов DKIM + SPF, используется цифровая подпись DNS зоны DNSSEC. Все это поддерживается в связке BIND 9 + Postfix + OpenDKIM.

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

Регламент

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

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

Чего почитать простому админу перед увольнением:

https://ru.wikipedia.org/wiki/Sender_Policy_Framework

https://ru.wikipedia.org/wiki/DomainKeys_Identified_Mail

https://ru.wikipedia.org/wiki/DNSSEC

понедельник, 29 июня 2020 г.

Дезинфекторы

У нас на производстве стали плодиться дезинфекторы. Думаю, если пойдет так дальше, под дезинфекторы можно выделить отдельный раздел.

Мобильный дезинфектор

Это такое устройство с двумя колесиками, которое распыляет адскую смесь (дезинфектор) в воздухе. Для распыления используется компрессор на десять атмосфер. Воздух под давлением подается в форсунку, где смешивается с жидкостью и выходит наружу в виде струи (факела). Близким по принципу работы устройством является - краскопульт для нанесения краски.
Отдадим должное краскопультам. От краскопультов взяли конструкцию форсунки. Форсунки бывают разных видов: HP, HVLP... HP (High Pressure) – высокое давление. HVLP (High Volume Low Pressure) - высокий объем, низкое давление. Мы работали с обоими типами. В серию пошли форсунки высокого давления, потому что нужно было обеспечить дальнобойность. Факел получили длиной 3-4метра. Однако, форсунки низкого давления тоже имеют право на существование, когда необходимо обрабатывать поверхности с расстояния до 1 метра.
По заданию, устройство должно работать в двух или трех режимах, по расходу дезинфектора (хим. раствора): "туман", "увлажнение", "орошение",.. "автомойка". Регулирование расхода химии выполняется с использованием маленького дозирующего насоса: мембранного или перистальтического. Подача воздуха от компрессора - через клапан электромагнитный. Устройство во время работы подсвечивает факел узконаправленным светом, мощным светодиодным светильником встроенным в форсунку - это для удобства, чтобы факел был виден.
Устройство включает сенсорный графический интерфейс (TFT-дисплей с тач-панелью) на контроллере. Рассчитывает время и циклы обработки по объему помещения и норме расхода. Пишет файл журнала.
Устройство может применяться для обработки общественных помещений, контейнеров, грузов, транспортных средств.

Стационарный дезинфектор для мобильного пункта КТ

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

Проходной шлюз для хим. обработки сотрудников

Это другая технология. В основе лежит концепция адиабатических систем увлажнения воздуха с распылением жидкости под высоким давлением в 30 атм. В системе используются специальные мини форсунки расположенные рядами сверху, снизу и по бокам шлюза, создают облако тумана вокруг объекта. Шлюз работает в автоматическом режиме по датчику присутствия и двери. При закрытых дверях включается обработка на заданный интервал. Шлюз имеет систему дренажа, отработанная жидкость (конденсат) собирается в бак. Объем шлюза подсвечивается.
Шлюз - мобильный, собирается как туристическая палатка с тентом на каркасе.
Используется по назначению для обработки персонала в защитных костюмах.

UV-C дезинфектор отработанного воздуха

Это перспективная разработка, выполняется по заказу медиков. Основана на новой технологии светодиодов диапазона UV-C (жесткий ультрафиолет). Пока что у светодиодов UV-C (255нм, 275нм) низкая квантовая эффективность и низкий показатель срока службы в сравнении со светодиодами видимого диапазона. Однако, технология не стоит на месте, мы ожидаем, что уже через год светодиодные систему УФ дезинфекции воздуха будут превосходить аналогичные люминесцентных лампах, как это было со светодиодным освещением. UV-C излучение убивает вирусы с высокой эффективностью. Система которую мы проектируем обрабатывает вдыхаемый/выдыхаемый воздух с эффективностью от 90% до 99%. Заказана система (приставка) для обработки отработанного воздуха в установку ИВЛ. Система монтируется в "выхлопную" трубу (трубку), представляет собой линейку светодиодов с последовательным запуском. Светодиоды включаются в момент выдоха и подсвечивают "бегущей волной" убийственного света выхлоп воздуха внутри трубки.
Задача - убить всех микробов и вирусов на выдохе, с одного прохода воздуха. Никакая система УФ дезинфекции на люминесцентных лампах не может обеспечить компактности. Данная технология позволяет делать носимое оборудование и встраивать UV-C дезинфекторы прямо в респиратор. Сегодняшний уровень технологии UV-C LED (УФ светодиодов) позволяет это сделать. Стоит дороже ватно-марлевой повязки. НО предполагаю, что к концу года цены станут адекватными, чтобы за 30-50у.е можно было выпустить маски со встроенным УФ-дезинфектором.
На сегодня мы запустили в производство опытные экземпляры электроники: регуляторы на 3 зоны облучения и 9 светодиодных линеек с диммированием со стабилизацией тока светодиода, светодиодные модули на UV-C - платы на металлическом основании. Кроме того, подготовили к выпуску модуль светодиодный для UV-C линз LEDIL. Линзы позволяют сфокусировать излучение и использовать УФ-дезинфектор для обработки продуктов на конвейере. Получается УФ-бактерицидная завеса.

Есть потребность? Обращайтесь. Я чувствую, в нашем арсенале не хватает еще одного типа дезинфектора.

четверг, 26 марта 2020 г.

Почему оптимизация memset накрылась...

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

Такую же причуду удалось получить на функции memcpy. GCC чудит и внутрь вставляет рекурсию, функция вызывает себя. Это происходит при использовании оптимизации -O3.

Отключить оптимизацию можно ключом -fno-tree-loop-distribute-patterns.

понедельник, 23 марта 2020 г.

ELF32 для ARM Cortex-M

Мы изучили идеи, которые вложены в загрузчик GNU/Linux. Поигрались с возможностями компилятора (gcc, clang) и линковщика (ld). Готовы описать, как выглядит сборка проекта и как ее упростить. Наша цель встроить в контроллер ARM Cortex-M4 функцию динамической загрузки приложений и загрузку динамических библиотек. Зачем? Чтобы менять программу в процессе работы. Я вынашиваю идеи как делегировать на контроллер обработку данных.

[IHI0044] ELF for the Arm Architecture
ELF format
ELF32 -- бинарный формат для представления объектных и перемещаемых файлов, исполняемых файлов, разделяемых данных и динамических библиотек. Формат состоит из разделов (сегментов). Исполняемый файл может содержать загружаемые сегменты (программы). Сегменты содержат: таблицы символов, строки с именами объектов и функций, а также необходимые данные для сборки кода (таблицы перемещений).
Таблица перемещения - описывает, как нужно модифицировать переменную или вызов функции, с учетом того, что на этапе компиляции не определены абсолютные адреса переменных и адреса функций. Т.е адреса должны подставляться в скомпилированный код на этапе сборки. Фактически код после компилятора представляет из себя шаблон, в который нужно вписать ссылки. Для каждой целевой платформы определяется, как именно вписывается ссылка в сегмент кода и сегмент данных - правила модификации кода. Наша целевая платформа - ARM с системой команд Thumb32. Арм определил ряд правил сборки кода [IHI0044]. Этот список содержит около сотни разных правил. Правила могут быть сложные, сложнее, чем просто записать адрес по указанному смещению. Обычно правила выглядят примерно так: Взять содержимое инструкции, декодировать инструкцию, к аргументу инструкции добавить смещение - адрес сегмента, вычесть адрес самой инструкции, и закодировать инструкцию обратно. Переменные (данные, data) и функции (код, text) могут размещаться в разных сегментах, поэтому при сборке возникают дополнительные операции.
Мы действуем методом обратной разработки. На входе уже есть готовый ELF формат с таблицей перемещений. В формате используются всего четыре разные правила. Заметим, что загрузчик исполняется в составе операционной системы под совершенно определенную архитектуру. Мы будем обрабатывать только те правила, которые возникают на данной архитектуре. Таким образом, число правил существенно уменьшается. Получаем всего 2 правила сборки:
R_ARM_ABS32 -- это правило абсолютного перемещения. Data: (S + A) | T. Сложить содержимое ячейки и адрес сегмента. Для адреса функции не забыть правильно указать флаг системы команд.

R_ARM_THM_CALL -- это правило кодирования инструкции вызова функции Thumb32: ((S + A) | T) - P (маска=0x01FFFFFE) \sa R_ARM_THM_PC22. Сложить содержимое ячейки и адрес сегмента и вычесть адрес самой ячейки. Для ссылки на функцию добавить флаг системы команд Thumb - единичку в младшем разряде. Такое же правило форматирования используется для вызовов R_ARM_THM_JUMP24. Тут мы используем некоторое упрощение- длина вызова - 22 бита вместо 24, два бита не обрабатываем. На вопрос почему, ответ: потом что так обрабатывает линковщик GNU binutils, такое упрощение рекомендует ARM.

Компилятор GCC использует правило R_ARM_TARGET1 для формирования таблицы инициализации - таблицы конструкторов и деструкторов. Это правило может быть таким же, как и R_ARM_ABS32, мы определяем как выглядит это правило, потому что инициализация - это дело разработчика операционной системы, наше дело. Мы используем тот же метод, R_ARM_ABS32.

Компилятор Clang использует еще пару правил. Эти правила возникают из идеи грузить константы 32 бита в два приема - отдельно младшую (lower16), отдельно старшую часть (upper16).

Сделал утилиту, которая работает в точности как readelf (GNU binutils), разбирает таблицы символов и перемещений. Зачем сделал, не очень понимаю. Результатом является возможность реализации динамической загрузки и динамической линковки кода в контроллере, разбор elf формата. Практического приложения пока нет. Пока разбирался с возможностями линковки увидел возможность собирать операционку с использованием статических библиотек. Статические библиотеки можно прикомпиливать целиком или только те фрагменты, которые используются в сборке, на которые ссылается код. Таким образом без ущерба для целостности удалось ужать объем ядра. Если функции не используются они в образ не попадают.

четверг, 20 февраля 2020 г.

DSP на OpenCL -- строительные блоки и блочные диаграммы

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

Блоки в языке OpenCL C - это такие чудные вложенные функции со шляпкой. В GCC поддерживается расширение языка Си, т.н. "nested functions", похожие сущности -- гнездящиеся функции. В языке OpenCL C 2.0+ можно на стороне устройства описывать из блоков асинхронный процесс обработки.

Достаточно крупные операции в DSP надо делать inline функциями, вернее давать возможность компилятору оптимизировать(векторизовать) большие блоки кода. Один из основных методов дающих существенный прирост производительности - векторизация циклов. Это значит что одиночные циклы, для которых число итераций определяется снаружи следует писать в виде inline функций. Мне нравится набор функций, который предлагает ARM CMSIS DSP. Но есть оговорки и вопросы оптимизации под аппаратуру. У нас есть задача переносить коэффициенты и модели из настольного компьютера в контроллер. В компьютере мы используем векторизацию float, в устройстве управления и сбора данных, сервере автоматизации используем ARM Neon, планируем использовать OpenCL. Нам нужен объект и формат данных для переноса "тензоров" и организации "потока тензоров" средствами протокола обмена, между платформами, нам нужен единый API (оптимизированный под конкретную архитектуру) для представления блочных диаграмм. Концептуально, мы говорим про внешний интерфейс для представления функционального блока FunctionalBlockObject, объекты-контейнеры данных TensorObject. А также мы будем рассматривать некоторый набор имен программных блоков, представленных в библиотеке DSP, чтобы ссылаться на их использование в блочной диаграмме.
Говоря про формат данных для внешнего представления блочных диаграмм, мы на сегодня будем ориентироваться на форматы принятые в BACnet: правила кодирования ASN.1, формат Json и CSML. Эти форматы мы выбираем прежде всего, чтобы не плодить сущности, поддержка этих форматов уже присутствует в составе нашего ПО.

Тензоры

Тензор - это массив данных определенного типа: float, double, int, char. Мы говорим про внешнее представление данных. Тнезор обладает свойством выравнивания данных - это его внутреннее свойство, определяется автоматически при создании объекта. Тензор имеет порядок укладки данных при записи и чтении, формат записи. Порядок записи включает: число каналов (feature), пакет данных (batch), пространственную размерность (spatial). Например, изображение с камеры имеет размерность 2d (spatial), организовано как поток кадров (batch), может иметь несколько каналов (feature). Звук имеет размерность =1 (spatial), организовано в поток по времени (batch), и разделено на каналы (features). Тензор - это единица обрабатываемой информации, с точки зрения функционального блока DSP или NN.
Мы будем исходить из того, что передать тензор целиком в одном пакете может оказаться не возможно. Поэтому обращение к тензору может быть как к массиву, ReadRange - запрос на чтение фрагмента массива. Мы предполагаем некоторую процедуру заполнения входных данных, объектов связанных между собой блочной диаграммой, затем запуск преобразования. Результатом преобразования может являться набор тензоров. Результат преобразования может распространяться по подписке CoV. Преобразование запускается по готовности входных данных, т.е. когда все входные данные(тензоры) записаны.
Тензоры могут быть вида TensorOutputObject, TensorInputObject. TensorOutputObject - обладает командным свойством на PresentValue. Когда записывается содержимое тензора, происходит запуск зависимостей по готовности.

Функциональные блоки

Функциональный блок содержит ссылку на исполняемый код(kernel) - DeviceObjectReference, список аргументов, список фиксированный идентификаторов объектов ARRAY_OF(DeviceObjectReference) - зависимости. В процессе исполнения, блок оформляет подписку на соответствующие объекты списка зависимостей. Функциональный блок имеет статус(состояние,контекст) - это некоторая структура данных, в которой могут хранится промежуточные величины. Контекст сохраняется и может быть перезаписан в режиме out-of-service.
Для блока мы определяем фазы: Инициализация, обновление контекста, обработка, пауза/завершение/сериализация контекста. Блочная диаграмма это не сам конвейер, это пакетное задание, которое может выполняться на конвейере с использованием данных. Конвейер это очередь исполнения по готовности.
Функциональный блок может обладать отладочной информацией (опция), для отладки используются штампы времени: queued, submit, start, end.

Программы

Программа - это текст или бинарник. Со стороны хост-приложения возможно определить из бинарника список блоков и список агрументов(входных/выходных портов). Каждый порт описывается 'форматом записи', однако формат записи из программы не определить.

Строительные блоки

Список функций API
dot(a,b)
операция свертки вектора sum_{0}^{N-1}(a[i]·b[i])
conv(a,b,n)
операция частичной свертки, выполняется с разворотом вектора sum_{0}^{N-1}(a[i]·b[n-i-1]). В это определение попадает корреляционная функция
float cascade_biquad_df2T(z,x)
вычисление значения каскада y = z+b0·x
float cascade_biquad_df2T_update(S,x)
обновление регистров каскада (S)
float fir(S,x)
FIR фильтр. Интерполяция.
pid(S,x)
PID регулятор.
lms_update(S,x)
LMS подстройка параметров адаптивного фильтра
IIR фильтры обрабатывают batch составляющую, время.
Для задач мультирейт используются фильтры LPF_decimate (по времени).
Аудио процессор может обрабатывать полосы и сшивать их на выходе. Для этого используются фильтры Allpass. В каких-то случаях могут использоваться фильтры iir_lattice.
Пид регулятор - это фильтр типа biquad. LMS фильтр - подстройка.

PID регулятор

 

Передаточная функция:
H(s) = ( KP + KI·1/s + KD·s ) Замена: 1/s = T/2·(1+z-1)/(1-z-1), s = 1/T·(1-z-1)

Передаточная функция в дискретном пространстве:
H(z) = KP + KI·T/2·(1+z-1)/(1-z-1) + KD/T·(1-z-1)
или 
H(z) = KP + KI·T·1/(1-z-1) + KD/T·(1-z-1)
Ki := KI·T
Kd := KD/T
Может быть представлена биквадратной секцией

void pid_init(S, x){
  S->b0 =  Kp + Ki + Kd;

  S->b1 = -Kp -2Kd
  S->b2 =  Kd;

  S->z1 = 0;
  S->z2 = 0;
}
// первый проход
float pid(S, x){
  return b0*x+z1;
}
// второй проход, обновление регистров
void pid_update(S, x){
  z1 =(b1-b0)*x + (z2 - z1);
  z2 = b2*x;
}
При разработке PID регулятора мы также используем два прохода: первый проход встроен в поток обработки, второй - поток обновления регистров и состояния фильтров. Пока что у нас получается что абсолютно все фильтры представлены выражением: y = g*x+Z, включая и PID регулятор. Далее мы хотим добавить свойство автоматической подстройки каэффициентов FIR фильтров и IIR фильтров. На втром цикле обновления регистров, мы добавляем возможность подстройки коэффциентов по методу LMS (). Есть известное преобразование FIR фильтра в LMS - получается адаптивный фильтр. Мы создаем сеть, в которой каждый функциональный блок должен обладать адаптивностью. Точнее можно сказать, что в процессе работы возникают периоды обучения - подстройки коэффициентов модели. Таким образом, фаза UPDATE может сопровождаться "обучением", если известна величина ошибки выходного сигнала. Для каждого слоя сети, мы можем определить этот признак "обучения".

LMS фильтры

Фильтры адаптивные (LMS) мы будем изготавливать из фильтров типа IIR и FIR. Метод LMS можно применять для обучения (подстройки весовых коэффициентов) нейронной сети. Мы вводим единую методику обучения/распознавания для нейронных сетей (NN) и для структур цифровой обработки сигналов (DSP). Нам известен способ изготовления адаптивного цифрового фильтра из FIR фильтра эту методику мы хотим расширить на остальные задачи.


Для некоторого преобразования c передаточной функцией представленной полиномом H(z) = h0+h1z-1+h2z-2+...+hNz-N , можно подстраивать коэффициенты, если известна ошибка выходного значения. Ошибка может быть пересчитана в поправку каждого конкретного коэффициента.
Ошибку на выходе определим, как квадрат отклонения.
E = (d[n] - y[n])2 = (d[n] - sum(b[i]z-i))2
изменение коэффицента для фильтра можно определить как градиент. Используем метод "градиентного спуска"
w[n+1]=w[n]-μ·∂E/∂w
Здесь μ - это величина шага подстройки. Для стабильного проявления методики, бывает полезно ограничить ошибку в каких то пределах, используя функцию clamp(). Берем частную производную, не сильно увлекаясь частностями. wi[n+1] = wi[n] - 2μ·e[n]·∂e[n]/∂wi
∂E/∂wi = -2(e[n])(x[n-i])

bi[n+1]=bi[n]+2μ (e[n])(x[n-i]).
Это все! Каждый коэффицент подстраиваем исходя из ошибки на выходе. IIR фильтр - это комбинация двух фильтров FIR. Для подстройки коэффициентов A(z) в IIR-фильтре с передаточной функцией H(z)=B(z)/A(z) используем аналогичный вывод:
ai[n+1]=ai[n]+2μ (e[n])(y[n-i]).

{Тут в изложении присутсвует некоторое искусство, с которым надо поразбираться. Работать будет в любом случае, однако математически правильный вариант брать частную производную не от уравнения в конечных разностях. Варианты, которые можно получить исследуя конечные разности и интегралы от дельта-функций - это линии задержки разного вида и фильтры на поток ошибок. К сожалению, быстро ответ не нашел. Нашел, что этим когда-то занимался I.-D. Landau, }
{1976 Feintuch’s algorithm}
{A parallel-type adaptive algorithm is proposed which utilizes the fast least-squares method. Its computational complexity is much less than that of I. Landau's (1978) method, which is based on hyperstability theory Hyperstability requires much computation because it involves matrix operations. The proposed method has nothing to do with hyperstability. It is also shown theoretically and using computer simulation, that the convergence performance of the proposed method is better than that of the series-parallel-type fast least-squares method.}
Линия задержки будет разная для разных форм представления фильтров. Глядя на презентацию, мне почему-то кажется, что тут может быть ошибка в интерпретации, что-то не то с частной производной. Должна быть простая форма. Мне видится простая форма при вторичном использовании регистров z. Предположительно выглядит так, доказать надо на практике:

void lms_dt2t_update(S, x, ue){
  b0+= ue*x;
  b += ue*x;// вектор
  a += ue*y;// вектор
}
void lms_update(S, x, ue){
  b0+= ue*x;
  b += ue*z;// вектор
  a += ue*z;// вектор
}
Если мы применяем LMS к фильтру типа PID должны быть веса разные для коэффициентов. Коэффциенты A не подстравиваются.

суббота, 15 февраля 2020 г.

DSP на OpenCL - Изоморфные преобразования алгоритмов

Хочу поделиться знаниями, вернее акцентировать внимание на известном факте.. Каково место цифровой обработки сигналов в коде. Хотим затолкать цифровую обработку в поток OpenCL, выполнять обработку на потоке данных. При этом не дает покоя одна мысль. Мы хотим получить наивысшую скорость работы и малые задержки. Поэтому мы каждый алгоритм DSP хотим разделить на две части. В первой части рассчитываем результат за минимальное число операций, во второй обновляем значения переменных. Я начну писать, потом видно будет.
В разделе DSP есть два вида фильтров: IIR- и FIR- фильтры. Прежде всего, оговорюсь все это относится к анализу LTI-систем, для которых можно описать передаточную функцию Y(s)=H(s)X(s) в непрерывном виде и далее перейти к дискретному виду передаточной функции H(z) путем замены переменной s←2/T(z-1)/(z+1). Передаточную функцию можно представить в виде отношения полиномов B(z)/A(z). Т.е.в виде дроби (дробно-рациональной функции) - сверху и снизу дроби полином z. Причем буква Z-1 на самом деле представляет из себя оператор задержки в дискретном по времени процессе. Как вообще люди в процессе своего образования доходят до этого - загадка. Думаю все сознательные электронщики и программисты должны это знать, но полагаю мало, кто знает. Наше образование не было полным в том виде, в каком оно может служить инструментом программиста или электронщика. Статья для избранных. На самом деле я уже начинаю сомневаться, что вообще этому учат. Когда-то учили анализу систем, учили передаточным функциям. Предлагаю, сейчас найти учебник на английском (по теме Linear time-invariant system и Digital Signal Processing). Для работы, чтобы получить готовый набор коэффициентов использую GNU Octave или MATLAB с модулем Signal, он умеет готовить матрицы коэффициентов. Использую конспекты лекций по курсу DSP & Digital Filters, потому что мне это близко и понятно, именно в таком виде. Автор курса ссылается на учебник.
[#] "Digital Signal Processing" Sanjit K. Mitra, 4th edition, McGraw Hill, 2011. ISBN 0071289461
Для стационарных систем (LTI) передаточная функция — это дробно-рациональная функция комплексной переменной s. Тернистый путь познания лежит через преобразования Фурье и Лапласа. Однако более короткий путь можно найти через дифференциальные уравнения в конечных разностях, только этот путь неведом выпускникам со специальностью кибернетика, так учат физиков. У меня в образовании был курс ТФКП(теория функции комплексного переменного) и курс математической физики (который символизирует наезд дельта функцией на уравнения), от которого ничего тут использовать не буду.
1) Систему можно представить в виде дробно рациональной функции N-ого порядка в странном пространстве наполненном операторами задержки и умножителями, -- это понятно программисту. Причем порядок числителя не выше чем у знаменателя. 2) Безумную дробь можно разложить на полюса и нули(полюса - это решения полинома в знаменателе, нули - решение полинома в числителе), путем факторизации полинома(еще одно умное слово - разложить на скобки, - так говорят в школе). Длинную дробно рациональную функцию большого порядка можно представить в виде суммы или произведения этих самых дробно-рациональных функций второго порядка. (точка!! - никакой больше теории, только картинки).
Иногда я чувствую, что являюсь хранителем некоего сокровенного невербального знания :). Могу доказать теорему не произнося слова, а рисуя только квадратики, треугольнички и стрелочки. Передаточные функции в дискретном пространстве можно строить в виде графов и представлять математические преобразования над функциям, как тождественные преобразования графов. Граф в точности отражает архитектуру программной или аппаратной реализации алгоритма цифровой обработки. На картинке буквой Z обозначаем оператор задержки, фактически это будет переменная (буфер, в аппаратуре - D-триггер), в которой хранится значение с предыдущего шага обработки.

Рис.1 Биквадратный фильтр. Прямое представление 1-типа DF1 (слева) и прямое представление 2-типа DF2 (справа)

Рисунок #1 получен путем прямой записи выражения B(z)X = (b0+b1·z-1+b2·z-2)X и затем выражения A(z)Y = (1-a1·z-1-a2·z-2)Y, где z-1 - оператор задержки на такт.
В первом случае (DF1)сложность, 5 умножений и 4 сложения: 5m+4a. Во втором (DF2) сложность та же, но для получения хранения промежуточных значений достаточно всего двух регитров Z.

Рис.2 Биквадратный фильтр. Прямое представление 2-типа DF2 (слева) и прямое представление второго типа транспонированное DF2T (справа)

Метод Транспонирования цифрового фильтра:
1) Меняем направление обхода, все стрелочки и треугольники разворачиваем. 
2) Слева будет выход, справа вход. Меняем местами буквы X и Y.
3) Разворачиваем буферы Z.
4) Где были узлы ставим "плюсик", где плюсики - ставим узел.
Цель можно считать достигнутой, мы развернули фильтр таким образом (DF2T), что на первом проходе можно получить в одно действие результат: Y = b0·X+z0 От входа X до выхода Y одно суммирование и одно умножение. Сложность первого прохода: 1m+1a. Этот проход совершается в той же самой математике, что и процесс свертки слоя нейронной сети.
К чему все это?. А, к нейролизации цифровой обработки сигналов в дискурсе для непосвещенных. Все что ниже первой линии будем считать частью процесса обучения, т.е. расчитывать и обновлять значения Z будем на втором проходе сети.
Какие у нас открываются возможности при использовании параллельной архитектуры. Правильно, фильтры могли бы считать паралельно, а не последовательно. Привычное разложение передаточной функции высокого порядка - в каскад (произведение) фильтров второго порядка. Каскад - это последовательная обработка фильтров, предполагает загрузку одного ядра/потока. Возможно представление в виде суммы, можно разложить в сумму фильтров первого и второго порядка. Cмотрю на структуру каскада фильтров. Каскад расчитывается, как sum(z1k+b0kX). Переобозначим переменную, S=sum(z1k) -- сумму можно расчитывать на втором проходе.


Рис.3 Каскад биквадратных фильтров.

Получаем финально такой результат: каскад фильтров можно представить на первом проходе Y =Z+ gX. Весь фильтр в одно действие!
// Алгоритм 
float cascaded_iir_df2t(float x){
// первый проход
   y = g*x+s;
   return y;
}
float cascaded_iir_df2t_update(float x){
// второй проход, вычисляем s
   s=0;
   x=g*x;
   for (i=0; i<len; i++){
     y = x+z1[i];
     z1[i] = b1[i]*x - a1[i]*y + z2[i];
     z2[i] = b2[i]*x - a2[i]*y;
     s+=z1[i];
     x = y;
   }
   return s;
}

Преобразования фильтров


Рис.4 Тождественные преобразования. a) замена переменных, b) объединение идентичных буферов, c) свойство линейности преобразования.

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

На рисунке #6 представлена одиночная операция транспонирования над одним каскадом фильтра FIR.
Рис.6 Преобразования фильтра FIR. Транспонирование каскада.

Векторизация фильтров

Следующий пункт - оптимизация под оборудование. Возможно будет лучше считать в векторном типе с использованием явной векторизации. Исходно мы используем транспонированную форму фильтра.
Рис.7 Прямая и транспонированная форма FIR-фильтра

Ниже пример реализации FIR фильтра. Также разделен на два прохода.
// FIR фильтр с явной векторизацией
float  fir_dft(float x){
  y = g*x+z.s0;
  return y;
}
float4 fir_dft_update(float4 z, float x){
  Z0=(0.0f);
  z = shuffle2(z,Z0, (int4){1,2,3,4});
  z =z + H*x;
  return z;
}
В данном варианте расчета, для обновления регистров фильтра нужно всего две процессорных инструкции shuffle и fma.

четверг, 6 февраля 2020 г.

Функция активации в слое нейронной сети

В статье дано описание функций активации. Выбран класс гладких функций удобный для реализации на CPU с поддержкой инструкций AVX512 и на GPU Intel Gen8-Gen11. При работе с GPU мы ориентируемся на язык OpenCL C. Дается некоторое количество утверждений, основанных на предположениях автора. Мы бездоказательно утверждаем, что использование класса гладких Smoothie-функций дает достаточную точность, для получения адекватных результатов работы сети.

Применяют два варианта функции активации в промежуточных слоях, остальные сводятся к производным. Мы хотим обобщить понятие о функции активации. В обобщенном виде к функции предъявляется ряд требований, таких как гладкость, монотонность и уровни насыщения. Функция задана на области определения и насыщается при определенных значениях. Важен внешний вид функции и важно, чтобы производная функции работала, давала соответствующий градиент в правильном направлении. Величина градиента не сильно важна и сказывается только на времени сходимости при обратном распространении ошибки.
Вариант 1) Сигмоида. Функция типа ступенька с размытыми краями. В обобщенном виде сигмоида может быть заменена на полное отсутствие сигмоиды лишь бы насыщалась на уровне 0 и 1. Т.е. вполне реально активацию можно заменить на операцию умножения с насыщением. Производная сигмоиды: s'(x) = s(x)(1-s(x)) выражается через значение самой функции, т.е. не надо знать ничего, кроме значения на выходе.
Вариант 2) Тангенс гиперболический. Это такая функция, которая монотонно возрастает имеет один перегиб. Имеет область значений [-1,+1].Все что реально нужно знать про эту функцию - она монотонная, нечетная и производная этой функции: (1+t(x))(1-t(x)). Тангенс гиперболический можно выразить через сигмоиду и наоборот. По сути это одна и та же функция.
Надо ли говорить, почему функция сигмоида и тангенс гиперболический не подходит под описание?! Потому что мы говорим про насыщение, а насыщение у сигмоиды не наблюдается. Это значит что нейрон будет "обучаться"(подстраивать коэффициенты) только в области, где производная не равна нулю.
Давайте предположим, что традиция использовать экспоненту в функции обусловлена исключительно требованием гладкости и дифференцируемости функции, которая в свою очередь продиктована методом вычисления градиента. Давайте предположим, что нам удаться сбалансировать систему таким образом, чтобы она не уходила в насыщение в процессе обучения.

tanh(x) = (1-e-2x)/(1+e-2x)
sigm(x) = 1/(1+e-x)
relu(x) = ln(1+ex) -- SmoothReLU, softplus, выпрямитель с гладким переходом.
Следует отметить, что производная от SmoothReLU дает сигмоиду. Функцию tanh можно выразить через сигмоиду и наоборот:

tanh(x) = (1-e-2x)/(1+e-2x) = 2sigm(2x) - 1
sigm(x) = 0.5*(tanh(x/2)+1)
Может правильнее было бы сказать, решением какого дифференциального уравнения являются эти функции. Тут самое время включить мозги, алгоритмы обучения сводятся к минимизации ошибки. Функции активации можно вводить не аналитически, а через уравнение. Если для каждого шага обучения можно оценить величину ошибки, функция ошибки может вводится, как отклонение от уравнения. Задача обучения будет сводится к минимизации этой ошибки.

d t(x)/d(x) = (1+t(x))(1-t(x))
d s(x)/d(x) =    s(x)·(1-s(x))
-- Производная имеет смысл градиента.

Кусочно-полиномиальные приближения функции

Мы не хотим использовать экспоненту категорически, поэтому предлагаем ориентироваться на приближенное вычисление сигмоиды или вообще заменить функцию активации на что-то приблизительно по форме напоминающее сигмоиду - ступеньку с размытыми краями. Все нормальные разработчики используют такой подход. Табличная функция активации выгядит примерно так:
// Кусочно-полиномиальные приближения функций
float activation(float x)
{
const int N=8;
   x = clamp(X_MIN, X_MAX,x);
// округляет в меньшую сторону до целого числа.
   idx = floor(x)-X_MIN;
   t  = x-floor(x);// дробная часть
   v0 = Coeff0[idx];
   v1 = Coeff1[idx];
   v2 = Coeff2[idx];
// квадратная интерполяция 
// v(t) = t*(t*v2+v1)+v0
   v = fma(t, v2,v1);
   v = fma(t, v, v0);
   return v;
}
Примерно в таком виде Intel описывает некоторую оптимизацию в собрании сочинений Intel Architectures Optimization Manual, глава 7.5.2 Fused post GEMM, Sigmoid Approximation with Minimax Polynomials. Которая оптимизация сводится к интерполяции полиномами второй степени. Этот вариант запасной. Он работает, не имеет противопоказаний, это доказано на опыте.
Можно определить функцию активации следующим образом.
На диапазоне [0+eps,1-eps] - происходит линейный рост с производной =1.
(0-eps, 0+eps) и (1-eps, 1+eps) -- сглаженный переход к насыщению - квадратным полиномом.
Назовем функцию, которая осуществляет плавный переход к насыщению Smoothie-выпрямитель. В масштабе 1.0 Smoothie-выпрямитель мы определим на интервале (-1.0,1.0) со значениями производной на конце интервала 0 и 1, соответственно.
Квадратичная интерполяция функции Smoothie-выпрямитель
F (x)  =  a2 x^2 + a1 x + a0
F'(x)  = 2a2 x   + a1
// граничные значения функции
F ( 1) =  a2 + a1 + a0 = 1
F (-1) =  a2 - a1 + a0 = 0
// производная непрерывная
F'( 1) = 2a2 + a1      = 1
F'(-1) =-2a2 + a1      = 0

a0 = 1/4, a1 = 1/2, a2 = 1/4, 

На области определения (-1,1)
SmoothieRect (x) = (x+1)(x+1)·0.25
SmoothieRect'(x) = (x+1)·0.5
В общем виде этот пример можно рассматривать, как метод вычисления коэффициентов полинома на интервале. Мы задаем производную и значение функции в каждой строчке таблицы. Таким образом рассчитываем коэффициенты интерполяции полиномом второй или третьей степени.
Позволю себе небольшое отступление для обоснования расчета коэффициентов табличного метода. Для каждого интервала нам нужно рассчитать коэффициенты полинома. Для этого нам нужно решить систему линейных уравнений полученных подстановкой граничных условий на интервале (0,1). Кто не помнит, надо посмотреть метод Гаусса-Жордана для решения систем линейных уравнений.
// Решение системы линейных уравнений.
Запишем функцию и его производную. Надо найти коэффициенты ai
F (x)  =  a3 x3 + a2 x2 + a1 x + a0
F'(x)  = 3a3 x2 +2a2 x  + a1
Шаг 1. Строим матрицу из множителей коэффициентов полинома ai, 
при заданном x и заданном значении функции bi:
  F(0) = (0, 0, 0, 1 | b0)  
 dF(0) = (0, 0, 1, 0 | b1)  
  F(1) = (1, 1, 1, 1 | b2)  
 dF(1) = (3, 2, 1, 0 | b3)  
Шаг 2. Вычитаем так чтобы получилась диагональная матрица.
b3 = b3-3b2
  F(0) = (0, 0, 0, 1 | b0)  
 dF(0) = (0, 0, 1, 0 | b1)  
  F(1) = (1, 1, 1, 1 | b2)  
 dF(1) = (0,-1,-2,-3 | b3-3b2)  
b2 = -2b2+b3
  F(0) = (0, 0, 0, 1 |  b0)  
 dF(0) = (0, 0, 1, 0 |  b1)  
  F(1) = (1, 0,-1,-2 |-2b2+ b3)  
 dF(1) = (0,-1,-2,-3 |  b3-3b2) 
b3 = b3-3b2+2b1+3b0
b2 = 2b0+b1-2b2+b3
          a3 a2 a1 a0
  F(0) = (0, 0, 0, 1 |  b0 )  
 dF(0) = (0, 0, 1, 0 |  b1 )  
  F(1) = (1, 0, 0, 0 | 2b0+ b1-2b2+b3 )  
 dF(1) = (0,-1, 0, 0 | 3b0+2b1-3b2+b3 ) 

Итого получаем: 
a0 = b0;
a1 = b1;
a3 = 2b0+ b1-2b2+b3
a2 =-3b0-2b1+3b2-b3
Далее мы рассмотрим вариант функции tanh-like активации заданный на интервале (-1,1) полиномом третьей степени.

Гладкие функции активации

Обращаю внимание на существование иных функций, smoothie-функций, которые подпадают под определение функции активации. Например, кубическая интерполяция Эрмита, функция smoothstep, S(x)=3x2 - 2x3. Эта функция введена в язык OpenCL и достойна применения. Функция имеет хорошую производную S'(x) = 6(x - x2) = 6x(1-x). Даже производная похожа на сигмойду, sigmoid-like function.
Продемонстрируем вывод полиномиальной функции активации, представляем функцию активации полином третьей степени. Приравниваем производную вне области определения (-1, 1) нулю и приравниваем значение вне области определения -1 и 1 соответственно, чтобы получить tanh-like функцию.
Вычисляем коэффициенты полинома для граничных условий:
Для полинома и его производной
T (x) = a3 x3 + a2 x2 + a1 x + a0
T'(x) =3a3 x2 +2a2 x  + a1
Составляем систему уравнений для граничных условий
T ( 1) =   a3 +  a2 + a1 + a0 = 1
T (-1) = - a3 +  a2 - a1 + a0 =-1 
T'( 1) =  3a3 + 2a2 + a1      = 0
T'(-1) = -3a3 + 2a2 - a1      = 0
Из выражения для производной получаем:
a1 = -3a3, a2 = 0.
Подставляем в выражение для функции:
T ( 1) = -2a3 + a0 = 1
T (-1) =  2a3 + a0 =-1
Получаем набор коэффициентов:
a0 = 0, a1 = 3/2, a2 = 0, a3 = -1/2
Функция Tanh-like на множестве (-1,1)
T (x)= (3x - x3)1/2 = 2*Smoothstep(x + 1/2)-1
T'(x)= (1  - x2)3/2 = 3/2 (1-x)(1+x)

Интерполяция Sigmoid-like на (-1,1).
S (x)= (T(x)+1)/2 = Smoothstep(x + 1/2)
S'(x)=  T'(x)/2
На языке OpenCL C кубическую интерполяцию функции активации Sigmoid и Tanh можно выразить с использованием базовой функции smoothstep, которая может быть реализована аппаратно или хорошо оптимизирована для векторных типов. В спецификации Khronos указано, что функции реализуются с использованием инструкций fma и mad.
T(x) = 2*smoothstep(-1.0, 1.0, x+0.5)-1.0;
S(x) = smoothstep(-1.0, 1.0, x+0.5);

float smooth_tanh(float x){
  x = clamp(x, -1.0f, 1.0f);
  return (3.0f - x*x)*x*0.5f;
}
float smooth_sigm(float x){
  x = clamp(x, -1.0f, 1.0f);
  return (3.0f - x*x)*x*0.25f+0.5f;
}
float smooth_step(float edge0, float edge1, float x){
  float t;
  t = clamp((x - edge0)/(edge1 - edge0), 0, 1);
  return t * t * (3 - 2 * t );
}
Коль скоро затронули тему Smoothie-функций, можем сюда же определить смузи-функцию выпрямителя, SmoothReLU. Производная от SmoothReLU должна давать функцию SmoothStep, а производная от ReLU - функцию Step, пороговую активацию. Строго говоря, SmoothReLU должна представлять собой полином четвертой степени.

Тестирование

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

Заключение

Я хочу поставить току в обсуждении. Самая простая функция активации - это ее отсуствие. Это и есть самая эффективная реалиазция.

float activation_s(float x){
  return clamp(-1.0f, 1.0f, x);
}

float activation_u(float x){
  return clamp(0.0f, 1.0f, x);
}

float activation_relu(float x){
  return fmax(0.0f, x);
}

Все что нам нужно - ограничивать значение на выходе.
Кстати, чтоб развеять мистический туман вокруг встроенной функции clamp, скажу, что ее можно реализовать двумя функциями: fmax и fmin, которые в свою очередь представлены соответствующими векторными инструкциями.
// Эмуляция функции clamp
float clamp(float x, float minval, float maxval) { 
  return fmin(fmax(x, minval), maxval); 
}
Такой вид функции активации позволяет легко перейти к реализации алгоритмов в целых числах. Активация в целых числах - это просто насыщение.
Пример такой операции оптимизированной под AVX-512:

activation_u8(a,b,c,d){
// добавление константы относится к свертке.
// a = _mm512_adds_epi32(a,a0);
// b = _mm512_adds_epi32(b,b0);
// c = _mm512_adds_epi32(c,c0);
// d = _mm512_adds_epi32(d,d0);
 v0= _mm512_packus_epi32(a, b);
 v1= _mm512_packus_epi32(c, d);
 return _mm512_packus_epi16(v0, v1);
}

activation_i16(a,b){
// a = _mm512_adds_epi32(a,a0);
// b = _mm512_adds_epi32(b,b0);
 return _mm512_packs_epi32(a, b);
}


суббота, 25 января 2020 г.

AVX512-VNNI -- Синтез алгоритма для сверточной нейронной сети

Необходимо сделать самый быстрый алгоритм для свертки (convolution). На свертках делается современное представление о нейронных сетях, так и называется сверточные нейронные сети (CNN). Производители процессоров заявляют: мы добавили искусственный интеллект (AI) в процессоры. По сути это означает только одно - увеличили длину вектора для операции умножения. Операции свертки используются в цифровых фильтрах обработки сигналов(DSP) и в задачах связанных с обработкой звука и изображения. Есть задачи выделения контуров изображений и задачи выделения движущихся объектов на изображении. Есть задачи увеличения четкости изображения, пересчет цветности. Все подобные задачи цифровой обработки изображений требуют матричного умножения NxN элементов. При этом коэффициенты могут быть заданы числами со знаком 8бит. А значения над которыми производится операция - 8 бит без знака. Под этот класс задач компания Intel разработала систему команд AXV512-VNNI. В систему команд VNNI входит всего две инструкции. Одна из них производит покомпонентное умножение 8битных целых без знака на 8битные целые со знаком и последующее сложение с накоплением в 32 битном приемнике.
Математически выражается как-то так. Знаковые и беззнаковые компоненты перемножаются и суммируются группами по четыре, чтобы не нарушать разрядность вектора. При большом количестве таких операций может возникать переполнение разрядов. Для нейронных сетей абсолютно все выводы классификации делаются на переполнении, если есть переполнение то это и означает переход от шаманства с бубнами к логике. В данном случае операция насыщения поверх сложения не сильно оправдана. Надо выполнить более 16 тыс таких операций, чтобы сложились условия для переполнения.
VPDPBUSDS — Multiply and Add Unsigned and Signed Bytes with Saturation
p1 = ZERO_EXTEND(a.byte[4*i+0]) * SIGN_EXTEND(b.byte[4*i+0])
p2 = ZERO_EXTEND(a.byte[4*i+1]) * SIGN_EXTEND(b.byte[4*i+1])
p3 = ZERO_EXTEND(a.byte[4*i+2]) * SIGN_EXTEND(b.byte[4*i+2])
p4 = ZERO_EXTEND(a.byte[4*i+3]) * SIGN_EXTEND(b.byte[4*i+3])
DEST.dword[i] = SIGNED_DWORD_SATURATE(DEST.dword[i] + p1 + p2 + p3 + p4)
Так же выглядит вторая операция - умножение знаковых 16 битных чисел с накоплением результата в 32 битном приемнике.
VPDPWSSDS — Multiply and Add Signed Word Integers with Saturation
p1 = SIGN_EXTEND(a.word[2*i+0]) * SIGN_EXTEND(b.word[2*i+0])
p2 = SIGN_EXTEND(a.word[2*i+1]) * SIGN_EXTEND(b.word[2*i+1])
DEST.dword[i] = SIGNED_DWORD_SATURATE(DEST.dword[i] + p1 + p2)
Вот я это так и вижу. Первая операция - для умножения матриц, для обработки видео. Вторая - для умножения векторов, для обработки звука.
Я буду исследовать производительность Throughput & Latency(CPI - clock per instruction) с использованием LLVM-MCA, инструмента. Подбираю, как должен выглядеть алгоритм на языке высокого уровня, чтобы выжать из этих инструкций максимум производительности.
Под рукой у меня мануал на 5 тысяч страниц с описанием процессорных инструкций накопленных за всю историю Intel [Intel® 64 and IA-32 Architectures Software Developer’s Manual, Volume 2:Instruction Set Reference] и таблица производительности инструкций на ядрах 10th Gen Intel Core, с микро-архитектурой Ice Lake. Кроме того, практикую магию с использованием компиляторов GCC 9.2 и Clang 9.0. Вывожу результат компиляции в ассемблер и анализирую производительность выделенных фрагментов полученного кода инструментом LLVM-MCA, анализатором кода.

$ gcc -march=icelake-client -O3 -S -o nn.s nn.c
$ llvm-mca.exe -mcpu=icelake-client -timeline nn.s
Производительность заявленная в таблице - одна инструкция с разрядностью 512 бит на такт, т.е. 64 шт умножений. Пишу цикл свертки по вектору (по горизонтали изображения) и по вертикали.
#include <x86intrin.h>
static inline
__attribute__((__target__("avx512vnni")))
void mac_vec_512(__m512i *acc, uint8_t * src, const __m512i* weight, const size_t length)
{
    int x=0;
    for (; x<length/64;x++){
        __m512i w = _mm512_load_si512(&weight[x]);
        __m512i v = _mm512_loadu_si512(&src[x*64]);
        acc[0] = _mm512_dpbusd_epi32(acc[0], v,w);
    }
    if (length&63) {
        __mmask64 mask = (1ULL <<(length&63)) -1;
        __m512i w = _mm512_load_si512(&weight[x]);
        __m512i v = _mm512_maskz_loadu_epi8(mask, &src[x*64]);
        acc[0] = _mm512_dpbusd_epi32(acc[0], v,w);
    }
}
При разработке алгоритма я исхожу из положения - компилятор самостоятельно уберет из процедуры проверки констант и выполнит развертывание цикла. В предыдущей статье я уже объяснял основной прием работы с масками, повторяться не буду. Если размерность вектора не кратна размеру регистра, используется маска. В данном случае по маске выполняется загрузка вектора.
Ниже представлен результат анализа производительности цикла, временная диаграмма. Цикл разернут компилятором.

[0,0]     DeeeeeeeeER    .    .    .    .    vmovdqu32   (%rdx), %zmm2
[0,1]     D========eeeeeeeeeeeER   .    .    vpdpbusd    (%r8), %zmm2, %zmm0
[0,2]     D===================eER  .    .    vmovdqa64   %zmm0, (%rcx)
[0,3]     .DeeeeeeeeE-----------R  .    .    vmovdqu32   64(%rdx), %zmm3
[0,4]     .D===========eeeeeeeeeeeER    .    vpdpbusd    64(%r8), %zmm3, %zmm0
[0,5]     .D======================eER   .    vmovdqa64   %zmm0, (%rcx)
[0,6]     . DeeeeeeeeE--------------R   .    vmovdqu32   128(%rdx), %zmm4
[0,7]     . D==============eeeeeeeeeeeER.    vpdpbusd    128(%r8), %zmm4, %zmm0
[0,8]     . D=========================eER    vmovdqa64   %zmm0, (%rcx)
Видно, что ожидаемой производительности не получили, результат инструкции VPDBUSD выходит каждые четыре такта.
То что происходит дальше - шаманство с бубном. Я беру фрагмент ассемблерного кода и подбираю руками комбинацию инструкций, которая дает максимальную производительность. Я знаю что конвейер способен загружать из памяти два регистра одновременно.

[0,0]     DeeeeeeeeER    .    .    .    .    vmovdqu32  (%rdx), %zmm2
[0,1]     DeeeeeeeeER    .    .    .    .    vmovdqa64  (%r8), %zmm3
[0,2]     D========eeeeER.    .    .    .    vpdpbusd   %zmm3, %zmm2, %zmm0
[0,3]     .D===========eER    .    .    .    vmovdqa64  %zmm0, (%rcx)
[0,4]     .DeeeeeeeeE----R    .    .    .    vmovdqu32  64(%rdx), %zmm2
[0,5]     .DeeeeeeeeE----R    .    .    .    vmovdqa64  64(%r8), %zmm3
[0,6]     . D==========eeeeER .    .    .    vpdpbusd   %zmm3, %zmm2, %zmm0
[0,7]     . D==============eER.    .    .    vmovdqa64  %zmm0, (%rcx)
[0,8]     . DeeeeeeeeE-------R.    .    .    vmovdqu32  128(%rdx), %zmm2
[0,9]     .  DeeeeeeeeE------R.    .    .    vmovdqa64  128(%r8), %zmm3
[0,10]    .  D=============eeeeER  .    .    vpdpbusd   %zmm3, %zmm2, %zmm0
[0,11]    .  D=================eER .    .    vmovdqa64  %zmm0, (%rcx)
При этом производительность не изменилась, но существенно скоратилась задержка. Т.е. выигрыш будет на сравнительно небольших векторах. Видно что задержка в 4 такта определяется исключительно зависимостью от выходного результата в регистре ZMM0. В следующем приближении я уменьшаю зависимость от выходного регистра.
$ llvm-mca.exe -mcpu=icelake-client -timeline nn1.s
[0,0]     DeeeeeeeeER    .    .    .    .    vmovdqa64    (%r8), %zmm4
[0,1]     DeeeeeeeeER    .    .    .    .    movdqu32    (%rdx), %zmm5
[0,2]     D========eeeeER.    .    .    .    vpdpbusd     %zmm4, %zmm5, %zmm3
[0,3]     .DeeeeeeeeE---R.    .    .    .    vmovdqa64    64(%r8), %zmm4
[0,4]     .DeeeeeeeeE---R.    .    .    .    vmovdqu32    64(%rdx), %zmm5
[0,5]     .D========eeeeER    .    .    .    vpdpbusd     %zmm4, %zmm5, %zmm2
[0,6]     . DeeeeeeeeE---R    .    .    .    vmovdqa64    128(%r8), %zmm4
[0,7]     . DeeeeeeeeE---R    .    .    .    vmovdqu32    128(%rdx), %zmm5
[0,8]     . D========eeeeER   .    .    .    vpdpbusd     %zmm4, %zmm5, %zmm1
[0,9]     .  DeeeeeeeeE---R   .    .    .    vmovdqa64    192(%r8), %zmm4
[0,10]    .  DeeeeeeeeE---R   .    .    .    vmovdqu32    192(%rdx), %zmm5
[0,11]    .  D========eeeeER  .    .    .    vpdpbusd     %zmm4, %zmm5, %zmm0
Получил пропускную способность 64 умножения на такт. Результат операции выдается на каждый такт. В данном случае производительность ограничена пропускной способностью загрузки данных. А вот дальше происходит невозможное. Оказывается, производительность можно УДВОИТЬ, если коэффициенты для последующих умножений равны или загружены в регистры заранее. Вот пример такой оптимизации:
$ llvm-mca.exe -mcpu=icelake-client -timeline nn1.s
[0,0]     DeeeeeeeeER    .    .    vmovdqu32      (%rdx), %zmm5
[0,1]     D========eeeeER.    .    vpdpbusd       %zmm4, %zmm5, %zmm3
[0,2]     DeeeeeeeeE----R.    .    vmovdqu32      64(%rdx), %zmm5
[0,3]     D========eeeeER.    .    vpdpbusd       %zmm4, %zmm5, %zmm2
[0,4]     .DeeeeeeeeE---R.    .    vmovdqu32      128(%rdx), %zmm5
[0,5]     .D========eeeeER    .    vpdpbusd       %zmm4, %zmm5, %zmm1
[0,6]     .DeeeeeeeeE----R    .    vmovdqu32      192(%rdx), %zmm5
[0,7]     .D========eeeeER    .    vpdpbusd       %zmm4, %zmm5, %zmm0
# Еще один вариант оптимизации удвоенной производительностью
[0,0]     DeeeeeeeeeeeER .    .    vpdpbusd         (%rdx), %zmm9, %zmm0
[0,1]     DeeeeeeeeeeeER .    .    vpdpbusd       64(%rdx), %zmm9, %zmm1
[0,2]     D=eeeeeeeeeeeER.    .    vpdpbusd      128(%rdx), %zmm9, %zmm2
[0,3]     .DeeeeeeeeeeeER.    .    vpdpbusd      192(%rdx), %zmm9, %zmm3
[0,4]     .D=eeeeeeeeeeeER    .    vpdpbusd      256(%rdx), %zmm9, %zmm4
[0,5]     .D=eeeeeeeeeeeER    .    vpdpbusd      320(%rdx), %zmm9, %zmm5
[0,6]     . D=eeeeeeeeeeeER   .    vpdpbusd      384(%rdx), %zmm9, %zmm6
[0,7]     . D=eeeeeeeeeeeER   .    vpdpbusd      448(%rdx), %zmm9, %zmm7
Чтобы получить максимальную производительность надо добиться, чтобы выход компилятора выглядел именно так. Я оговорюсь. Результат может быть недостоверным. Надо проверять на целевой архитектуре. Следует отметить, что при оптимизации небольших циклов практически все векторные инструкции есть смысл отделять от инструкций загрузки регистров, так код получается быстрее, хотя и выглядит немного избыточным, у планировщика появляется больще возможностей загрузить конвейер неоднородными операциями.
#include <x86intrin.h>
static inline
__attribute__((__target__("avx512vnni")))
void mac_vec_512(__m512i *acc, uint8_t * src, const __m512i* weight, const size_t length)
{
    int x=0;
    for (; x<length/64;x++){
        __m512i v = _mm512_loadu_si512(&src[x*64]);
        acc[x&7] = _mm512_dpbusd_epi32(acc[x&7], v,weight[x]);
    }
    if (length&63) {
        __mmask64 mask = (1ULL <<(length&63)) -1;
        __m512i v = _mm512_maskz_loadu_epi8(mask, &src[x*64]);
        acc[x&7] = _mm512_dpbusd_epi32(acc[x&7], v, weight[x]);
    }
}
Тут мы используем восемь аккумуляторных регистров. Потом вектора аккумуляторов надо сложить между собой, выполнить горизонтальное суммирование. Такой красивый результат, как в предыдущей "ручной" опмтимизации, получить с помощью компилятора не удается. В большинстве случаев компилятор выдает что-то не эффективное. Сказывается еще одно аппаратное ограничение - тормозит декодер иснтрукций, у него своя производительность, ограниченная длиной буфера разбора. Не всегда этой длины хватает, чтобы обрабатывать четыре инсрукции в кодировке EVEX-512 на одном такте.
Наилучший результат без ручных оптимизаций __(впишу когда протестирую, хотя бы на симуляторе)__ умножений на такт. Теоретический предел 128 умножений на такт.

Нашел у Intel идею для описания подобной операции VDPDUSD там, где ее нет, на инструкциях AVX2. # VPMADDUBSW υ8×s8 → s16 multiples, VPMADDWD broadcast1 s16 → s32, and VPADDD s32 → s32.
// цикл на операциях AVX2
    for (x=0; x<length/32;x++){
        v = _mm256_loadu_si256((void*)&src[x*32]);
        v = _mm256_maddubs_epi16(v,weight[x]);
        v = _mm256_madd_epi16(v,E);
        acc= _mm256_add_epi32(acc, v);
    }
Для AVX2 компилятор выполняет оптимизацию сносно, алгоритм модифицировать не пришлось. Могут быть расхождения в результате, может возникать насыщение в операции VPMADDUBSW. В этом случае рекомендуется снижать разрядность весовых коэффициентов. Теоретический предел для AVX2 получается 32 умножения на такт, инструкции VPMADDUBSW и VPMADDWD имеют производительность две на такт, 0.5 CPI.

Эта статья об удвоении производительности. Сначала берем хорошую новую инструкцию в 2-4 раза более производительную чем раньше, получаем результат 1/4 от заявленной производительности. Потом начинаем шаманить и водить руками, получаем 100% - теоретический предел. Потом ударяем в бубен и получаем все 200%, о такой возможности производитель умолчал. Потом замеряем реальную производительность в цикле и получаем скромный результат всего __160%__ от заявленной в документации, что тоже не плохо. Еще остается методика оптимизации.

четверг, 23 января 2020 г.

Набор инструкций AVX-512 -- Переосмыслить всё!

-- переписать все алгоритмы на новую систему команд! - смысл такой.
Современные компиляторы никак не умнее программиста. Есть очень небольшая возможность оптимизировать циклы или ветвления за счет использования инструкций AVX-512. Одна из возможностей новой системы команд, на которую упирает разработчик (Intel) - возможность замены циклов и операций ветвления в циклах на инструкции с маской, это когда операция применяется не ко всему вектору, а к отдельным частям этого вектора. Вряд ли компилятор в ближайшие несколько лет научится полиморфить алгоритмы. Сейчас компиляторы работают по шаблону, шаблоны пишут люди, такие же программисты. Шаблонами все многообразие не описать. Изучая вторую неделю эффективность компиляторов версии 9.2+ понимаю, что они очень и очень далеки от искусственного интеллекта. В большинстве случаев мне почему то удается переписать код после оптимизации настолько, что он становится в два раза быстрее. Есть незначительное число вариантов циклов, которые поддаются оптимизации.
Краевые эффекты (хотел сказать "побочные"). Допустим мы можем эффективно умножать матрицы 8х8 элементов. Как умножать матрицы 7х7 элементов. Или 9х9 элементов. Вот тут на помощь приходит новая система команд с масками. Маска позволяет загрузить не кратное число элементов в безумно длинный векторный регистр. А вы знали, что векторная инструкция AVX-512 способна выдавать результат, сложения или умножения на каждом такте. Т.е. это буквально удваивает производительность, но только в том случае, когда алгоритм изначально можно параллелить. Если мы делаем из алгоритма AVX2 новый алгоритм увеличивая разрядность, такого прироста не наблюдается.
Копирование байтовых строк. Вот первая идея.

for(всех частей кратных векторному регистру){
  v = _mm512_loadu_si512(src); src+=64;
  _mm512_storeu_si512(m64, v)
}
if (есть остаток строки) {
// формируем маску
  __mmask64 mask = (~0ULL)>>(64-(len&63));
// загружаем строку произвольной длины меньше длины регистра
  __m512i v = _mm512_maskz_loadu_epi8(mask, src);
// сохраняем обработанную строку, фрагмент строки.
  _mm512_mask_storeu_epi8(dst, mask,v);
}

Для примера рассмотрим реализацию алгоритма BASE64. Чтобы не перегружать мозг, только основной цикл. Каждые 6 бит кодируются печатными буквами из словаря 64 символа. Кодировка применяется при передаче вложений в почтовые сообщения и веб-приложения. Короче, весь интернет забит этой кодировкой. Надо вам переслать скан договора или фотоархив, а он вот так вот кодируется BASE64 [RFC 4648].
// Реализация алгоритма кодирования BASE64
uint8_t* __attribute__((__target__("avx512vbmi")))
base64_enc_axv512(uint8_t *dst, uint8_t *src, int length)
{
const __m512i lookup = _mm512_loadu_si512(base64_table);
/* [a5..a0:b5b4][b3..b0:c5..c2][c1c0:d5..d0] -- расположение бит в памяти
   32-6,32-12,32-18,32-24 - сдвиги в младшем слове [31:0]
   64-6,64-12,64-18,64-24 - сдвиги в старшем слове [63:32]
   Вместе они формируют константу сдвигов.
*/
const __m512i shifts = _mm512_set1_epi64(0x282E343A080E141AULL);
const __m512i revert = (__m512i)(__v64qi){
    -1, 2, 1, 0, -1, 5, 4, 3, -1, 8, 7, 6, -1,11,10, 9,
    -1,14,13,12, -1,17,16,15, -1,20,19,18, -1,23,22,21,
    -1,26,25,24, -1,29,28,27, -1,32,31,30, -1,35,34,33,
    -1,38,37,36, -1,41,40,39, -1,44,43,42, -1,47,46,45
    };
__mmask64 mask = (1ULL<<48)-1;// маска 48 бит
    while (length>=48){// Производительность(CPI) 2 такта  на цикл 64 байта, ускорение 32 раза!!
        __m512i v = _mm512_maskz_loadu_epi8(mask, src);
        src+=48;
        v = _mm512_permutexvar_epi8 (revert, v);// переставить местами байты BSWAP32
        v = _mm512_multishift_epi64_epi8(shifts, v);// 32-6,32-12,32-18,32-24...
        v = _mm512_permutexvar_epi8(v, lookup);// игнорирует старшие 2 бита в байте.
        _mm512_storeu_si512(dst, v); dst+=64;
        length-=48;
    }
    . . .
}
Если реализация данного алгоритма вдохновляет, проследуйте по ссылке.. [0x80]. Для подготовки алгоритма использую справочник псевдо-функций по системе команд, стараюсь не списывать [Intel Intrinsics Guide].
Я понимаю, мир меняется очень серьезно. Компиляторы не дают ожидаемого эффекта, не поспевают за нововведениям Intel. Я вынужден использовать псевдо-функции (intrinsic).
Как вы понимаете, закон Мура в действии. Развитие идет по экспоненте, раньше частота удваивалась, потом удваивалась разрядность, потом ядреность, потом векторность, потом смысловая емкость системы команд. Что-то обязательно должно удваиваться, чтобы вы покупали все новые и новые гаджеты. Надо предлагать сразу несколько алгоритмов, потому что оптимизации подобного рода нужны на платформе ARM Neon, ARM Helium, Intel SSE, на AVX, на AVX2, AVX512 и так далее, каждые три года новая идея способная удвоить производительность. У программиста мозг не резиновый, программисты -- вымирающий вид. Сегодня это векторность в сочетании со смысловой емкостью системы команд.
Надо чтобы программа на этапе загрузки решала, какую оптимизацию использовать в зависимости от поддерживаемой системы команд. Многие куски кода, такие как функции кодирования, я повторно использую в проектах, могу иногда для разминки мозгов применять разные оптимизации. Эти оптимизации нужно оставлять в бестиариуме, прикомпиливать подходящие варианты.
Вернемся к обсуждению алгоритма. Данная реализация ускоряет в 32 раза процесс обработки, в сравнении с алгоритмом, который мне казался вполне быстрым, потому что тратил всего один такт на каждый байт. Данная реализация тратит два такта процессора на 64 байта. Попробую доказать. Ниже привожу результат анализа загрузки ресурсов ядра процессора Ice Lake.

Resource pressure by instruction:(нагрузка на вычислительные блоки)
[0]    [1]    [2]    [3]    [4]    [5]    [6]    [7]    [8]    [9]    Instructions:
 -      -      -      -     0.65   0.35    -     1.00    -      -     vpermb    (%rdx), %zmm2, %zmm0
 -      -      -     0.97    -      -      -      -     0.03    -     subl      $48, %r8d
 -      -     1.00    -      -      -      -      -      -      -     vpmultishiftqb    %zmm0, %zmm3, %zmm0
 -      -      -      -      -      -      -     1.00    -      -     vpermb    %zmm1, %zmm0, %zmm0
 -      -      -      -     0.02   0.33   1.00    -      -     0.65   vmovdqu64 %zmm0, (%rax)
 -      -      -     0.98    -      -      -      -     0.02    -     addq      $48, %rdx
 -      -     0.98    -      -      -      -     0.02    -      -     addq      $64, %rax
 -      -      -     0.02    -      -      -      -     0.98    -     cmpl      $47, %r8d
 -      -     0.04    -      -      -      -      -     0.96    -     jg        .L3


Timeline view:
                    0123456789
Index     0123456789          012

[0,0]     DeeeeeeeeeeER  .    .    vpermb (%rdx), %zmm2, %zmm0
[0,1]     DeE---------R  .    .    subl   $48, %r8d
[0,2]     D==========eER .    .    vpmultishiftqb %zmm0, %zmm3, %zmm0
[0,3]     D===========eeeER   .    vpermb %zmm1, %zmm0, %zmm0
[0,4]     .D=============eER  .    vmovdqu64      %zmm0, (%rax)
[0,5]     .DeE-------------R  .    addq   $48, %rdx
[0,6]     .DeE-------------R  .    addq   $64, %rax
[0,7]     .DeE-------------R  .    cmpl   $47, %r8d
[0,8]     .D=eE------------R  .    jg     .L3

[1,0]     . DeeeeeeeeeeE---R  .    vpermb (%rdx), %zmm2, %zmm0
[1,1]     . DeE------------R  .    subl   $48, %r8d
[1,2]     . D==========eE--R  .    vpmultishiftqb %zmm0, %zmm3, %zmm0
[1,3]     . D===========eeeER .    vpermb %zmm1, %zmm0, %zmm0
[1,4]     .  D=============eER.    vmovdqu64      %zmm0, (%rax)
[1,5]     .  DeE-------------R.    addq   $48, %rdx
[1,6]     .  DeE-------------R.    addq   $64, %rax
[1,7]     .  DeE-------------R.    cmpl   $47, %r8d
[1,8]     .  D=eE------------R.    jg     .L3

[2,0]     .   DeeeeeeeeeeE---R.    vpermb (%rdx), %zmm2, %zmm0
[2,1]     .   DeE------------R.    subl   $48, %r8d
[2,2]     .   D==========eE--R.    vpmultishiftqb %zmm0, %zmm3, %zmm0
[2,3]     .   D===========eeeER    vpermb %zmm1, %zmm0, %zmm0
[2,4]     .    D=============eER   vmovdqu64      %zmm0, (%rax)
[2,5]     .    DeE-------------R   addq   $48, %rdx
[2,6]     .    DeE-------------R   addq   $64, %rax
[2,7]     .    DeE-------------R   cmpl   $47, %r8d
[2,8]     .    D=eE------------R   jg     .L3
Результат получен на анализаторе загрузки ядра, LLVM-MCA для целевой платформы Ice Lake. Видно, что на каждый цикл алгоритма, см. строчки с [1.0] по [1.8] выходит результат, тратится всего два такта. Долго запрягаем, загружаем конвейер команд в работу, потом быстро-быстро грузим. Пропускная способность алгоритма определяется не задержкой, а тем как быстро проходится второй, третий, и последующие циклы. Первая таблица показывает загрузку по портам(вычислительным блокам). Вычислительные блоки не однородные, первые два - это делители, есть четыре унифицированных АЛУ, два канала загрузки из кеш и один канал сохранения в кеш. Векторные инструкции неравномерно распределяются между вычислительными блоками. Некоторые инструкции AVX могут работать параллельно. Если в спецификации на процессор написано два модуля AVX-512 это значит что инструкции AVX-512 могут выходить две за такт, в параллель по разным портам. Одна инструкция может быть представлена последовательностью микроинструкций, которые грузятся сразу на несколько вычислительных блоков. D- дешифрация команд, e - обработка на конвейере, E - исполнение, R - переназначение/сохранение регистров, = - простой планировщика, - простой порта.
В отсутствие возможности пощупать конкретную архитектуру или будущую систему команд, Intel предлагает отлаживать алгоритмы на эмуляторе Intel Software Development Emulator [SDE]. А оптимизацию предлагает тестировать на этом самом Анализаторе LLVM-MCA. Так и сделал. Алгоритм BASE64 отлаживал на эмуляторе SDE.

$ gcc -march=icelake-client -DTEST_BASE64 -O3 -o base64.exe base64.c
$ echo -n 'hello world!' | sde.exe -icx -- ./base64.exe

пятница, 27 декабря 2019 г.

CRC -- реализация через умножение без переносов

CRC - (Циклический избыточный код; англ. Cyclic redundancy check) используется для контроля целостности данных при передаче по сети или записи данных на носитель. Разные протоколы используют различные варианты, которые различаются разрядностью и образующими полиномами. Для рассчета CRC используется арифметика Галуа, в частности операция умножения. Умножение в конечном поле Галуа можно разделить на две операции: умножение без переносов и редуцирование. В большинстве случаев программисты не сильно задумываются над математикой и используют готовые таблицы и алгоритмы.
Допустим в нашем распоряжении есть процессорная инструкция умножения без переносов. Мы хотим оптимизировать и тем самым ускорить в несколько раз расчет контрольных сумм CRC, в десятки раз. Целью статьи является разработка общего подхода к синтезу алгоритмов содержащих умножение в конечном поле. Процесс создания алгоритмов для расчета CRC оказался весьма трудоемким. Хотя все эти алгоритмы подобны друг другу, но требуют отладки и расчета коэффициентов.
{Статья большая, я ее буду медленно дописывать и исправлять... прошу прощения за возможные неточности.}
На платформе Intel операции умножения без переносов соответствует инструкция pclmul. На платформе ARMv8-А также есть аналогичные инструкции, на базе которых можно эффективно реализовать операцию умножения без переносов, так называемая инструкция полиномиального умножения (PMULL) разрядностью операнда 64бит. На платформе ARMv7-A+NEON есть поддержка в виде инструкций vmull.p8, векторное полиномиальное умножение 8бит. На платформе ARMv8-М +MVE Helium ожидаем поддержку в виде векторных инструкций полиномиального умножения p8 бит и p16 бит

Все варианты операции CRC отличаются разрядностью, полиномами и порядком следования бит в числе.

В нашем программном обеспечении используются операции CRC:
CRC8/BACnet
CRC16/BAСnet
CRC16/MODBUS
CRC32K/Koopman
Кроме того, мы используем
CRC32С/Castagnoli
CRC32
CRC32B

Для описания алгоритмов я буду использовать инлайн-функции, которые призваны генерировать векторные инструкции. Базовые операции:
REFLECT128 -- Разворачивание последовательности бит в векторе
BSWAP128 -- Переворачивание порядка следования байт в векторе
LOAD128U -- Загрузка данных в векторный регистр без выравнивания
SLL128U -- логический сдвиг влево числа 128 бит.
Умножение без переноса:
CL_MUL32 () -- умножение без переносов двух чисел 32 бита, результат 64 бита
CL_MUL32H() -- старшие 32 бита произведения двух 32 битных чисел
CL_MUL32L() -- младшие 32 бита произведения
CL_MUL128() -- умножение младшей или старшей части вектора 64 бит на старшую или младшую часть второго вектора. Результат - число 128 бит.

Что такое CRC

По сути вот что это такое, с математической точки зрения,

CRC = (M(x) • xt) mod P(x), 
где t = deg(P(x)) - степень полинома, M(x)- сообщение(число)
-- сообщение M(x) представляют в виде длинного двоичного числа и подвергают операции редуцирования в конечном поле.
Идея оптимизации сводится к накоплению буфера с отложенной операцией редуцирования, операциям сдвига в конечном поле и последующему финальному редуцированию на буфере. Буфером является векторный регистр.
Традиционно обработка CRC ведется на потоке данных, по одному байту. В реальных задачах возникает необходимость расчета CRC при передаче данных по сети или записи на носитель, с расчетом по шапке и телу пакета. Расчет может производится по частям: сначала от шапки пакета, затем от тела пакета или сразу от всего пакета. Когда чтение данных осуществляется блоками, расчет может производится так же по блокам.
Мы хотим разработать скоростной алгоритм, обрабатывать до 16 байт на один цикл, ориентируясь на максимальную разрядность векторного регистра. Планируем использовать операции умножения без переносов. Мы хотим рассчитывать параллельно до 16 контрольных сумм пользуясь тем, что преобразования линейные. Результат параллельного расчета фрагментов сообщения можно складывать. Операция сложения представлена битовой операцией исключающего ИЛИ (⊕). Сложение может выполняться до или после редуцирования.
Редуцирование - это уменьшение разрядности числа после операции умножения без переноса.
Метод ускорения CRC с использованием инструкций умножения без переноса (pclmul) был предложен до нас[2,3], Intel. Однако в литературе представлен только метод для CRC-32. Кроме того, в предложенной методике операция редуцирования выполняется без использования оперций умножения. В данной статье я привожу варианты алгоритмов другой разрядности 8,16,24,32, 64, и нескольких образующих полиномов, привожу методику рассчета сдвиговых коэффициентов, и констант для редуцирования полиномов. Предложен свой вариант алгоритмов для отраженных полиномов, отраженного порядка бит.
Список литературы, рекомендую предварительно ознакомиться:
[1] A painless guide to CRC error detection algorithms
[2] Intel/ Fast CRC Computation Using PCLMULQDQ Instruction
и несколько реализаций на GitHub и в исходниках Apple, в частности BZIP2.
[3] Intel® Carry-Less Multiplication Instruction and its Usage for Computing the GCM Mode

Операция схлопывания

Сообщение М можно разложить на два фрагмента со смещением:
M1(x) • xs ⊕ M0(x) = M1 • (xs mod P) ⊕ M0
Эту операцию называют "схлопывание" (folding), слов не хватает, слова "редуцирование" и "сложение" уже заняты и наполнены смыслом. Сразу предлагаем воспользоваться правилом: редуцирование от редуцирования на результат не влияет. Поэтому мы можем сначала выполнить редуцирование сомножителей, а потом применять операцию умножения.
Константы (xs mod P) для разных показателей степени 32 64 96 128 .. можно рассчитать заранее и использовать для ускорения алгоритма.

Допустим, максимальное число, над которым удобно производить операцию "схлопывания" = 128 бит. Это прежде всего определяется разрядностью операции CL_MUL. В нашем случае операция может выполняться над 64 битными числами и результат имеет разрядность 128бит. Тогда определим операцию схлопывания так, преобразование вектора 256 бит [M4:M3:M2:M1] в 128 бит:
[C2:C1] = (M4 • K3) ⊕ (M3 • K4) ⊕ [M2:M1]
-- результат безрассудного умножения имеет разрядность 128 бит
K4 = x128 mod P
K3 = x192 mod P
K4 и K3 - это сдвиговые константы, возведение в степень - это операция сдвига на соответствующее число разрядов. Таким образом мы заменяем операцию сдвига в поле на операцию умножения без переноса. Соответственно, для последующего преобразования в CRC-32 нам понадобится операция редуцирования из 128-в-96 и 96-в-64 бит и финальное редуцирование 64-в-32 бита. Надо научиться рассчитывать константы для сдвигов и схлопывания.

Редукция Barrett'a

Редукция основана на такой идее, которая изначально родилась для целых чисел. Хотим подобрать такое число Ux, чтобы результат умножения Ux*Px - давал единицу с высокой точностью. Например, хотим поделить целое число в конечном поле 32 бита на 10, а операция деления напрочь отсутствует в нашем арсенале. Тогда мы можем найти такое целое число, чтобы результат операции (A*Ux)>>32 в точности давал бы результат A/10 для всех целых чисел. По сути мы ищем такое целое число, чтобы оно работало как multiplicative inverse (обратное число по отношению к образующему полиному).
Редукцию Барретта можно было бы определить так:

            T1 = floor(R(x)/x^32) * [1/P(x)];   R/P
            T2 = floor(T1/x^32) * P(x);         int(R/P)*P;
            CRC= (R+int(R/P)*P) mod x^32;       R-int(R/P)*P
только под делением мы понимаем операцию уполовинивания в поле, под сложением - битовую операцию XOR и под произведением - умножение без переносов, под делением на x^32 понимаем сдивиг на 32 разряда вправо. Короче, математическое определение редуцирования мало что дает на практике. Можно сказать, зная алгоритм редуцирования для целых чисел я угадал, как он будет выглядеть в конечном поле.
Операция  floor(R(x)/x^32) в конечном поле это просто сдвиг на 32 разряда вправо, а для отраженного полинома - сдвиг влево.

Функция barrett_calc рассчитывает магическое число Баррета, которое в дальнейшем позволяет редуцировать полином.
/*! Алгоритм №1. рассчитывает константу Барретт'а */
uint64_t barrett_calc(uint64_t poly, int bits)
{
    int n = bits;
    uint64_t r = poly;
    uint64_t v = 0;
    while (--n){
        if (r & (1ULL<<n)) {
            r ^= poly>>(bits-n);
            v |= 1ULL<<n;
        }
    }
    if (r) v|=1;// для деления с остатком округляем в большую сторону
    return v;
}
{Я совсем не помню откуда я взял этот алгоритм.}
Операция редуцирования методом Барретта 64-в-32
/*! Алгоритм №2. Редукция Барретт'а*/
uint32_t R (uint32_t v1, uint32_t v0)
{
    uint32_t t = CL_MUL32H(v1, U) ^ v1;// U имеет разрядность 33 бита
    return CL_MUL32L(t, P) ^ v0;
}
Константа U рассчитывается по алгоритму Барретта. Константа U имеет разрядность 33 бита. Поскольку старший бит в число 32 бита не укладывается, к результату умножения CL_MUL32H может добавляться значение операнда, - та самая единица в старшем разряде - это следует не потерять.

Такой же алгоритм можно использовать для обратного порядка обработки CRC, обратного порядка следования бит.
/*! Алгоритм №3. Редукция Барретт'а для отраженных чисел */
uint32_t RR (uint32_t v1, uint32_t v0)
{
    uint32_t t = CL_MUL32L(v1, UR);
    return CL_MUL32H(t, PR) ^ t ^ v0;// PR имеет разрядность 33 бита
}
Во втором случае используется отраженный полином Pr, с отраженным порядком следования бит. И отраженная константа Ur.

Мы только что определили, как выглядит операция поиска обратного (деления с остатком?) не-пойми-чего [1/P] с помощью алгоритма уполовинивания. см barrett_calc
Надо проверить результат barrett_calc. Если все верно, то произведение [1/P]•P должно давать что-то похожее на единицу с достаточной точностью. Проверяем...
(0x104D101DF • 0x104C11DB7) = 1.00000000.490D678D (65 бит);
Результат интересный, на единицу похож смутно, отдаленно, но похож. Если сдвинуть на 64 разряда, то единица получается. А вот в младших 32 битах лежит магическое число, для следующей операции пригодится - обратите внимание остаток равен сдвиговой константе.
Можно проверить ту же операцию на отраженном полиноме.
(0x1F7011641 • 0x1DB710641) = 1.63CD6124.00000001 (65 бит)
Чтобы поверить в магию умножения над отразить все биты в числе, все 65 бит... Короче, это тоже можно считать единицей, с заданной точностью 32 бита. Т.е. полиномы можно умножать задом наперед, но результат будет вывернут.

Правило для умножения без переноса для чисел с отраженным порядком бит выглядит так:
reflect(A) • reflect(B) == reflect(A • B)<<1

Результат нужно сдвинуть на один разряд. Таким образом, мы можем предложить алгоритм для отраженных полиномов, с использованием той же операции умножения без переносов. Алгоритм будет выглядеть очень похожим, с той лишь разницей, что при работе с отраженными полиномами всегда нужно сдвигать результат умножения на разряд влево.
Мой личный опыт строится на историческом осознании роли CRC и развитии CRC в нашем программном обеспечении, А также из чтения статей и изучения чужого опыта. В контроллерах мы применяем табличные методы расчета CRC. Причем, всегда стоит вопрос минимизации исходного кода прошивки, поэтому мы используем маленькие таблицы. Привожу пример табличного алгоритма, 4бит на шаг алгоритма.
// Алгоритм №4 -- табличный метод 4 бит.
uint32_t CRC32_update   (uint32_t crc, uint8_t val){
    crc = (crc<<4) ^ CRC32_Lookup4[((crc>>28) ^ (val>>4)) & 0xF ];
    crc = (crc<<4) ^ CRC32_Lookup4[((crc>>28) ^ (val   )) & 0xF ];
    return crc;
}
// Алгоритм №5 -- табличный метод 8 бит.
uint32_t CRC32_update   (uint32_t crc, uint8_t val){
    crc = (crc<<8) ^ CRC32_Lookup8[((crc>>24) ^ (val   )) & 0xFF ];
    return crc;
}
Математическое обоснование этих действий.

 [M3:M2:M1]•x8 ==  ((M3•x8 ⊕ M2)•x8 ⊕ M1)•x8
В примере сообщение нарезано по байтам. Видно, что алгоритм будет подобен себе для любой нарезки сообщения: на 1 бит, на 2, на4, на 8, на 16...
Для синтеза таблиц умножения использую метод...
// Алгоритмом № Расчет таблицы умножения, прямой порядок
/*! \brief генерация таблицы подстановки CRC
 \param poly полином
 \param bits число бит в полиноме
 \param size число элементов в таблице 16 или 256
 */
void crc_gen_table(uint64_t poly, int bits, int size)
{
    uint64_t table[size];// = {0};
    uint64_t p =poly;
    int i,j;
    table[0] = 0;
    table[1] = p;
    for (i=1;(1<<i)<size;i++)
    {
        if (p&(1ULL<<(bits-1))) {
            p &= ~((~0ULL)<<(bits-1));
            p = (p<<1) ^ poly;
        } else
            p = (p<<1);
        table[1<<i] = p;
        for(j=1; j<(1<<i); j++) {
            table[(1<<i)+j] = p ^ table[j];
        }
    }
// Распечатать таблицу
    printf("POLY=0x%0*"PRIX64"\n", bits/4, poly);
    for(i=0;i<size;i++){
        printf("0x%0*"PRIX64", ", bits/4, table[i]);
        if ((i&0x3)==0x3) printf("\n");
    }
}
Пользуясь Алгоритмом № умею синтезировать любые таблицы от CRC-8 до CRC-64. Для синтеза таблиц отраженного CRC использую аналогичный алгоритм. По сути таблица - это таблица умножения полинома на число 0,1..15, для 4битной таблицы и на число 0..255 для 8 битной таблицы. Таким образом, можно таблицу заменить на операцию полиномиального умножения с редуцированием по образующему полиному. В общем виде эти алгоритмы запишутся так:
Алгоритм № -- метод c умножением 8 бит.
[c4:c3:c2:c1] ::= [c3:c2:c1:0]⊕((c4⊕ Mi)•x8 mod P)


Как пользоваться редуцированием.
Один шаг над восьмибитными числами, обновляет значение CRC по байтам:

uint32_t CRC32_update_8 (uint32_t crc, uint8_t data){
    uint32_t val=data;
    val = BSWAP32(val); // перестановка байт
    crc = crc ^ val;
// Barrett's reduction
    uint32_t t = CL_MUL32H(crc>>24, 0x04D101DF)^(crc>>24);
    uint32_t v = CL_MUL32L(t      , 0x04C11DB7)^(crc<< 8);
    return v;
}
На самом деле многие алгоритмы и теоремы в конечном поле можно вовсе не доказывать. Доказательством может служить простой перебор чисел. Можно проверить утверждение для всех чисел в поле.
Я пишу серию алгоритмов для обработки по 8 бит, по 16 бит, по 32 по 64. И подбираю операции, отлаживаю для каждого шага. Потом из четырех полученных алгоримов создают пятый, который может работать со всеми длинами. Алгоритм можно вывести из математики, но на самом деле это две паралельные реальности, я смотрю на математику и пишу алгоритм, и наоборот проверяю какие-то действия в алгоритме и переписываю математику, которая этому соответствует. Что-то видно из математики что-то видно из алгоритма.

uint32_t CRC32_update_32 (uint32_t crc, uint32_t val){
 val = BSWAP32(val); // перестановка байт
 crc = crc ^ val;
// Barrett's reduction
 uint32_t t = CL_MUL32H(crc, 0x04D101DF)^(crc);
 uint32_t v = CL_MUL32L(t  , 0x04C11DB7);
 return v;
}
Вот этот вариант алгоритма проходит 32 бита (4 байта) за шаг. Ускорились в 4 раза!

// Алгоритм №. Расчет сдвиговых констант (xt mod P):
uint64_t xt_mod_P(uint64_t poly, int t, int bits) { uint64_t v = poly; do { if (v & MSB){ // старший бит проверяем v = (v<<1) ^ poly; } else { v = (v<<1); } } while (--t); return v; }
Отраженные сдвиговые константы можно получить путем отражения и сдвига результата операции.
/* Алгоритм №. Расчет отраженных сдвиговых констант (xt mod P):
  poly - отраженный полином 33 бита.
*/
uint64_t xt_mod_P_reflect(uint64_t poly, int t, int bits)
{
    uint64_t v = poly;
    do {
 if (v & 1){
     v = (v>>1) ^ poly;
 } else {
     v = (v>>1);
 }
    } while (--t);
    return v;
}


Кроме того, хочу предложить алгоритм расчета сдвиговой константы для сдвига вправо, x-t mod P(x). Этот вариант привожу для полноты картины.
// Алгоритм №. расчет сдвиговых констант, сдвиг вправо
uint64_t xt_mod_P_neg(uint64_t poly, int t, int bits){
    uint64_t v=1;
    do {
        if (v&1)
            v = (v>>1) ^ (1ULL<<(bits-1)|poly>>1);
        else
            v = (v>>1);
    } while (--t);
    return v;
}

После этого все операции сдвига я начал оценивать именно с этой позиции, каждое умножение в поле это сдвиги и набор сдвигов можно представить через операции умножения. Каждый байт сообщения (слово, двойное слово) мы подвергаем сдвигу в поле и в конце должны просуммировать все эти значения. Мы должны решить, как именно разбивать сообщение на группы, на байты и слова. В конечном итоге каждый байт исходного сообщения подвергается своему сдвигу и суммируется с результатом.
Еще одно интересное свойство, которе хотелось бы добавить для полноты изложения - можно предложить обратную операцию, отмена действия UNDO, операция CRC обратима. Можно предложить операцию дописывания шапки сообщения, независимо от тела. Сначала расчитывается CRC от тела сообщения, потом досчитывается CRC от шапки, и в итоге можно вычислить CRC от всего сообщения. И наконец, можно предложить обратное преобразование, которое находит ошибочный бит в сообщении -- восстановление сообщения. Ну и совсем хакерская мысль, можно предложить математическое преобразование, которое меняет сообщение, но сохраняет CRC неизменным. Можно переписать сообщение и подогнать значение CRC, добавив магию в соль.

Итак, начнем с голой идеи. Мне сложно сразу написать правильный, работающий алгоритм, практически невозможно отладить, отладить можно только каждое отдельное действие. Поэтому я создаю алгоритм по шагам. На первом шаге я делаю всего одну опрерацию - редуцирование, - отлаживаю ее.
На втором шаге делаю две операции: редуцирование и сдвиг в поле
M1•(x32 mod P) mod P.
Потом, отлаживаю алгоритм редуцирования и два сдвига,
M2•(x64 mod P)) ⊕ M1•(x32 mod P)) mod P.
Затем - четыре сдвига,
M4•(x128 mod P)) ⊕ M3•(x96 mod P) ⊕ M2•(x64 mod P) ⊕ M1•(x32 mod P) mod P.
Потом делаю алгоритм, который обрабатывает все сообщение целиком с циклическим схлопыванием и финальным редуцированием.
 Алгоритм №. Векторный алгоритм рассчета CRC-32
Шаг 1. загрузить вектор с переворотом 
 [c4:c3:c2:c1] := [M4:M3:M2:M1] ⊕ [crc:0:0:0]
Шаг 2. Цикл для каждого фрагмента 128-бит
    [A3:A2:A1] := [c4:c3]•(x192 mod P) ⊕ [c2:c1]•(x128 mod P)
 [c4:c3:c2:c1] := [ 0:A3:A2:A1] ⊕ [M4:M3:M2:M1]
Шаг 3. Сдвиги,
 [c3:c2:c1] := [с4:c3]•(x96 mod P) ⊕ [с2:c1]•(x32 mod P) 
    [c2:c1] := [c2:c1] ⊕ с3•(x64 mod P)
Шаг 4. Редуцирование 
    [T2:T1] := c2•Ux
    [c2:c1] := T2•Px ⊕ [c2:c1]
Возвращает: crc = c1

1. Вычисляем, например, CRC32 по словам 32 бита. Алгоритм №.
2. Пишу преобразование двойным словам 64 биата


uint32_t CRC32_update_64 (uint32_t crc, uint8_t* data){
 c = BSWAP64(data); // перестановка байт
 c = c ^ crc;       
// схлопывание [c2:c1]•x32 mod P(x) (96 бит в 64 бита)
 с = CL_MUL32(c[1], xt_mod_P(Px, 64))
   ^ CL_MUL32(c[0], xt_mod_P(Px, 32));
// Barrett's reduction
 uint32_t t = CL_MUL32H(crc, Ux)^(crc);
 uint32_t v = CL_MUL32L(t  , Px);
 return v;
}

// Алгоритм №. 
uint32_t    CRC32_update_128(uint32_t crc, uint8_t *data, int len) {
   uint64x2_t c;
   c = BSWAP128(data); data+=16;
   c[1] ^= (uint64_t)crc<<32;
   int blocks = (len+15 >> 4);
   while (--blocks>0) {
      c = CL_MUL64(c[1], K3) 
        ^ CL_MUL64(c[0], K4)
        ^ BSWAP128(data); data+=16;
   }
// схлопывание [c4:c3:c2:c1]•x32 mod P(x) (128бит в 64 бита)
   c = CL_MUL32(c[1]>>32, K4) 
     ^ CL_MUL32(c[1]>> 0, K5)
     ^ CL_MUL32(c[0]>>32, K6)
     ^ CL_MUL32(c[0]>> 0, K7);
// Редуцирование 64 бит в 32 бит
   uint32_t t = CL_MUL32H(c[0]>>32, Ux);
   crc = CL_MUL32L(t, Px) ^ c[0];
   return crc;
}
Далее, планируем увеличить разрядность схлопывания до 256 или 512. На этот алгоритм можно обратить внимание, если под рукой появилась векторная инструкция VPCLMULDQD разрядностью 256 или 512 бит.
Алгоритм №. Векторный алгоритм рассчета CRC-32, 256 бит
Шаг 1. загрузить вектор с переворотом 
 [c8:c7:c6:c5:c4:c3:c2:c1] := [M8:M7:M6:M5:M4:M3:M2:M1] ⊕ [crc:0:0:0:0:0:0:0]
Шаг 2. Цикл для каждого фрагмента 256-бит
 [A3:A2:A1] := [c4:c3]•(x256+64 mod P) ⊕ [c2:c1]•(x256 mod P)
 [A7:A6:A5] := [c8:c7]•(x256+64 mod P) ⊕ [c6:c5]•(x256 mod P)
 [c8:c7:c6:c5:c4:c3:c2:c1] := [ 0:A7:A6:A5: 0:A3:A2:A1] ⊕ [M8:M7:M6:M5:M4:M3:M2:M1]
Шаг 3. Редуцирование 256 в 128
 [c4:c3:c2:c1] := [c4:c3:c2:c1]  
    ⊕ [c8:c7]•(x128+64 mod P) ⊕ [c6:c5]•(x128 mod P)
Шаг 4. если есть фрагмент 128-бит,
    [A3:A2:A1] := [c4:c3]•(x192 mod P) ⊕ [c2:c1]•(x128 mod P)
 [c4:c3:c2:c1] := [ 0:A3:A2:A1] ⊕ [M4:M3:M2:M1]
Шаг 5. Сдвиги,
 [c2:c1] := c4•(x128 mod P) ⊕ c3•(x96 mod P)
    ⊕ c2•(x64 mod P) ⊕ c1•(x32 mod P) 
Шаг 6. Редуцирование 
 [T2:T1] := c2•Ux
 [c2:c1] := T2•Px ⊕ [c2:c1]
Возвращает: crc = c1

Теперь хотим развить тему и научиться синтезировать такие же алгоритмы для отраженного порядка бит, отраженных полиномов. Точно также, я начинаю с алгоритма редуцирования и по шагам отлаживаю каждую операцию, каждый сдвиг.
// Алгоритм №. 
uint32_t CRC32K_update_8 (uint32_t crc, uint8_t *data){
 uint32_t val;
 val = data[0];
 crc = crc ^ val;
 poly64x2_t c = {crc};
 c = SLL128U(c,32+24);// сдвиг влево числа 128 бит
// Barrett's reduction
    poly64x2_t t;
    t  = CL_MUL128(c, (poly64x2_t){Ux,Px}, 0x00);
    c ^= CL_MUL128(t, (poly64x2_t){Ux,Px}, 0x10);
    return c[1];
}
// Алгоритм №. 
uint32_t CRC32K_update_16(uint32_t crc, uint8_t *data){
 uint32_t val;
 val = *(uint16_t*)data;
 crc = crc ^ val;
// подбираю связь между величиной сдвига и шагом алгоритма 
 poly64x2_t c = {crc};
 c = SLL128U(c,32+16);
// Barrett's reduction
    poly64x2_t t;
    t  = CL_MUL128(c, (poly64x2_t){Ux, Px}, 0x00);
    c ^= CL_MUL128(t, (poly64x2_t){Ux, Px}, 0x10);
    return c[1];
}
// Алгоритм №. 
uint32_t CRC32K_update_24(uint32_t crc, uint8_t *data){
 uint32_t val=0;
 __builtin_memcpy(&val, data,3);
 crc = crc ^ val;
// подбираю связь между величиной сдвига и шагом алгоритма 
 poly64x2_t c = {crc};
 c = SLL128U(c,32+8);
// Barrett's reduction
    poly64x2_t t;
    t  = CL_MUL128(c, (poly64x2_t){Ux, Px}, 0x00);
    c ^= CL_MUL128(t, (poly64x2_t){Ux, Px}, 0x10);
    return c[1];
}
// Алгоритм №. 
uint32_t CRC32K_update_32(uint32_t crc, uint8_t *data){
 uint32_t val;
 val = *(uint32_t*)data;
 crc = crc ^ val;
//  сдвиги заменяем на умножение 
 poly64x2_t c = {0,(uint64_t)crc<<32};
 c = CL_MUL128(c& 0xFFFFFFFF00000000ULL, (poly64x2_t){K7<<1}, 0x01);
// Barrett's reduction
 poly64x2_t t;
 t  = CL_MUL128(c, (poly64x2_t){Ux, Px}, 0x00);// {Ux, Px}
 c ^= CL_MUL128(t, (poly64x2_t){Ux, Px}, 0x10);
 return c[1];
}
// Алгоритм №. 
uint32_t CRC32K_update_128(uint32_t crc, uint8_t *data){
 poly64x2_t c = (poly64x2_t)LOAD128U(data);
 c[0]^=crc;
// сдвиги
 c = CL_MUL128(c<<32, (poly64x2_t){K4<<1}, 0x00)
   ^ CL_MUL128(c& 0xFFFFFFFF00000000ULL, (poly64x2_t){K5<<1}, 0x00)
   ^ CL_MUL128(c<<32, (poly64x2_t){K6<<1}, 0x01)
   ^ CL_MUL128(c& 0xFFFFFFFF00000000ULL, (poly64x2_t){K7<<1}, 0x01);
// Barrett's reduction
 poly64x2_t t;
 t  = CL_MUL128(c, (poly64x2_t){Ux, Px}, 0x00);// {Ux, Px}
 c ^= CL_MUL128(t, (poly64x2_t){Ux, Px}, 0x10);
 return c[1];
}
Убедившись что все эти алгоритмы №..№ работают, перехожу к синтезу
// Алгоритм №. 
uint32_t CRC32K_update_N(uint32_t crc, uint8_t *data, int len){
 poly64x2_t c = (poly64x2_t)LOAD128U(data);data+=16;
 c[0]^=crc;
// Циклическое схлопывание folding
 int blocks = (len+15 >> 4);
 while (--blocks>0){
  c = CL_MUL128(c, (poly64x2_t){K3<<1}, 0x00)
    ^ CL_MUL128(c, (poly64x2_t){K4<<1}, 0x01);
  c = SLL128U(c, 32) ^ (poly64x2_t)LOAD128U(data); data+=16;
 }
// Сдвиги редуцирование
 c = CL_MUL128(c<<32, (poly64x2_t){K4<<1}, 0x00)
   ^ CL_MUL128(c& 0xFFFFFFFF00000000ULL, (poly64x2_t){K5<<1}, 0x00)
   ^ CL_MUL128(c<<32, (poly64x2_t){K6<<1}, 0x01)
   ^ CL_MUL128(c& 0xFFFFFFFF00000000ULL, (poly64x2_t){K7<<1}, 0x01);
// Barrett's reduction
 poly64x2_t t;
 t  = CL_MUL128(c, (poly64x2_t){Ux, Px}, 0x00);// {Ux, Px}
 c ^= CL_MUL128(t, (poly64x2_t){Ux, Px}, 0x10);
 return c[1];
}
Что хотелось бы доработать в алгоритме. Убрать все сдвиги. Например, если у меня получается что умножается число со сдвигом, то из двух операций можно оставить одну - умножение. Везде где мне вздумалось при выводе алгоритма написать CL_MUL(a, K<<n), это можно решить иначае CL_MUL(a, K') пересчитать константу. Если константа расчиыватеся, как xs mod P(x) Надо рассчитать константу xs-1 mod P(x)
 Алгоритм №. Векторный алгоритм CRC-32B отраженный порядок 
Шаг 1. загрузить вектор M, 32x4 =128 бит
 [c4:c3:c2:c1] := [M4:M3:M2:M1] ⊕ [0:0:0:crc]
Шаг 2. Цикл для каждого фрагмента 128-бит
    [A3:A2:A1] := [c2:c1]•(x191 mod P) ⊕ [c4:c3]•(x127 mod P)
 [c4:c3:c2:c1] := [A3:A2:A1:0] ⊕ [M4:M3:M2:M1]
Шаг 3. Сдвиги,
 [c2:c1] := 
      с1•(x127 mod P)
    ⊕ с2•(x95 mod P)
    ⊕ с3•(x63 mod P)
    ⊕ с4•(x31 mod P)
Шаг 4. Редуцирование 
    [T2:T1] := c1•Ux
    [c2:c1] := T1•Px ⊕ [c2:c1]
Возвращает: crc = c2
В данном Алгоритме № есть еще один сдвиг, лишняя операция, которую можно заменить на умножение. Число A покомпонентно в цикле сдвигается на 32 бита. Этот сдвиг тоже можно заменить путем коррекции коэффициентов K3=X191-32 и K4=X127-32.
 Алгоритм №. Векторный алгоритм CRC-32B отраженный порядок 
Шаг 1. загрузить вектор M, 32x4 =128 бит
 [c4:c3:c2:c1] := [M4:M3:M2:M1] ⊕ [0:0:0:crc]
Шаг 2. Цикл для каждого фрагмента 128-бит
    [A3:A2:A1] := [c2:c1]•(x159 mod P) ⊕ [c4:c3]•(x95 mod P)
 [c4:c3:c2:c1] := [0:A3:A2:A1] ⊕ [M4:M3:M2:M1]
Шаг 3. Сдвиги,
 [c2:c1] := 
      с1•(x127 mod P)
    ⊕ с2•(x95 mod P)
    ⊕ с3•(x63 mod P)
    ⊕ с4•(x31 mod P)
Шаг 4. Редуцирование 
    [T2:T1] := c1•Ux
    [c2:c1] := T1•Px ⊕ [c2:c1]
Возвращает: crc = c2

Привожу результат практический:


CRC-32/BZIP2 poly=0x04c11db7
init=0xffffffff refin=false refout=false xorout=0xffffffff check=0xfc891918
Px = 0x104C11DB7
Ux = 0x104D101DF
K7 = xt_mod_P(Px,  32) = 0x04C11DB7
K6 = xt_mod_P(Px,  64) = 0x490D678D
K5 = xt_mod_P(Px,  96) = 0xF200AA66
K4 = xt_mod_P(Px, 128) = 0xE8A45605
K3 = xt_mod_P(Px, 192) = 0xC5B9CD4C
K2 = xt_mod_P(Px, 512) = 0xE6228B11
K1 = xt_mod_P(Px, 576) = 0x8833794C
Таблица умножения и константы:

CRC-32B poly=0x04c11db7
init=0xffffffff refin=true refout=true xorout=0xffffffff check=0xCBF43926

POLY=0xEDB88320
Таблица умножения:
0x00000000, 0x1DB71064, 0x3B6E20C8, 0x26D930AC,
0x76DC4190, 0x6B6B51F4, 0x4DB26158, 0x5005713C,
0xEDB88320, 0xF00F9344, 0xD6D6A3E8, 0xCB61B38C,
0x9B64C2B0, 0x86D3D2D4, 0xA00AE278, 0xBDBDF21C,

BarrettR u = x^64/P(x) Ur=0x1F7011641 Pr=0x1DB710641
[1/p]*p(x) = 163CD612400000001
Test =CBF43926 ..ok

K7( 32-1) = 00000001
K6( 64-1) = B8BC6765
K5( 96-1) = CCAA009E
K4(128-1) = 9BA54C6F
K3(192-1) = 65673B46

CRC-32C (Castagnoli) poly=0x1edc6f41
init=0xffffffff refin=true refout=true xorout=0xffffffff check=0xE3069283

POLY=0x82F63B78
Таблица умножения:
0x00000000, 0x105EC76F, 0x20BD8EDE, 0x30E349B1,
0x417B1DBC, 0x5125DAD3, 0x61C69362, 0x7198540D,
0x82F63B78, 0x92A8FC17, 0xA24BB5A6, 0xB21572C9,
0xC38D26C4, 0xD3D3E1AB, 0xE330A81A, 0xF36E6F75,

BarrettR u = x^64/P(x) Ur=0xDEA713F1 Pr=0x105EC76F1
[1/p]*p(x) = 0DD45AAB800000001
Test =E3069283 ..ok

K7( 32-1) = 00000001
K6( 64-1) = DD45AAB8
K5( 96-1) = 493C7D27
K4(128-1) = 3171D430
K3(192-1) = 3743F7BD

CRC-32K/BACnet (Koopman) poly= 0x741b8cd7 (0xEB31D82E)
init=0xffffffff refin=true refout=true xorout=0xffffffff check=0x2D3DD0AE

POLY=0xEB31D82E
Таблица умножения:
0x00000000, 0x83CF0F3C, 0xD1FDAE25, 0x5232A119,
0x7598EC17, 0xF657E32B, 0xA4654232, 0x27AA4D0E,
0xEB31D82E, 0x68FED712, 0x3ACC760B, 0xB9037937,
0x9EA93439, 0x1D663B05, 0x4F549A1C, 0xCC9B9520,

BarrettR u = x^64/P(x) Ur=0x17D232CD Pr=0x1D663B05D
[1/p]*p(x) = 018C5564C00000001
Test =2D3DD0AE ..ok

K7( 32-1) = 00000001
K6( 64-1) = 18C5564C
K5( 96-1) = 9D65B2A5
K4(128-1) = 69F48E4D
K3(192-1) = CBB06D55

[31.12.2019]
Написал векторный алогоритм для CRC-16 с отраженным порядком.
Алгоритм №. Векторный алгоритм CRC-16B отраженный порядок
Шаг 1. загрузить вектор M, 16x8 =128 бит
 [c8:c7:c6:c5:c4:c3:c2:c1] := [M8:M7:M6:M5:M4:M3:M2:M1] ⊕ [0:0:0:0:0:0:0:crc]
Шаг 2. Цикл для каждого фрагмента 128-бит
          [A5:A4:A3:A2:A1] := [c4:c3:c2:c1]•(x191 mod P) ⊕ [c8:c7:c6:c5]•(x127 mod P)
 [c8:c7:c6:c5:c4:c3:c2:c1] := [A5:A4:A3:A2:A1:0:0:0] ⊕ [M8:M7:M6:M5:M4:M3:M2:M1]
Шаг 3. Сдвиги,
 [c2:c1] := 
      с1•(x127 mod P)
    ⊕ с2•(x111 mod P)
    ⊕ с3•(x95 mod P)
    ⊕ с4•(x79 mod P)
    ⊕ с5•(x63 mod P)
    ⊕ с6•(x47 mod P)
    ⊕ с7•(x31 mod P)
    ⊕ с8•(x15 mod P)
Шаг 4. Редуцирование 
    [T2:T1] := c1•Ur
    [c2:c1] := T1•Pr ⊕ [c2:c1]
Возвращает: crc = c2
Тут сообщение и вектор (с) разбиты на элементы по 16 бит. Общая идея ясна алгоритм подобен Алгоритму № для CRC-32B.
Привожу наборы констант.
CRC-16/X25 poly=0x1021
init=0xffff refin=true refout=true xorout=0xffff check=0x906e
POLY=0x8408
Таблица умножения:
0x0000, 0x1081, 0x2102, 0x3183,
0x4204, 0x5285, 0x6306, 0x7387,
0x8408, 0x9489, 0xA50A, 0xB58B,
0xC60C, 0xD68D, 0xE70E, 0xF78F,

BarrettR u = x^16/P(x) Ur=0x1911 Pr=0x10811
[1/p]*p = 00000000019D80001
Test =906E ..ok

K9( 15) = 0001
K8( 31) = 19D8
K7( 47) = 1CBB
K6( 63) = 042B
K5( 79) = 81BF
K4( 95) = 2C27
K3(111) = 8555
K2(127) = 7EEA
K1(191) = A95D
CRC-16/MODBUS poly=0x8005
init=0xffff refin=true refout=true xorout=0x0000 check=0x4b37

POLY=0xA001
Таблица умножения:
0x0000, 0xCC01, 0xD801, 0x1400,
0xF001, 0x3C00, 0x2800, 0xE401,
0xA001, 0x6C00, 0x7800, 0xB401,
0x5000, 0x9C01, 0x8801, 0x4400,

BarrettR u = x^16/P(x) Ur=0x1BFFF Pr=0x14003
[1/p]*p = 000000001D0020001
Test =4B37 ..ok

K9( 15) = 0001
K8( 31) = 9001
K7( 47) = FC01
K6( 63) = D101
K5( 79) = CCC1
K4( 95) = C551
K3(111) = C3FD
K2(127) = C100
K1(191) = CCD0
[01.01.2020]
Зачем публиковать, если нет ни одного ценителя. Ау! Ценители отзовитесь. Нужно публиковать алгоритм?
Алгоритм №. Векторный алгоритм CRC-64/XZ отраженный порядок
Шаг 1. загрузить вектор M, 64x2 =128 бит
 [c2:c1] := [M2:M1] ⊕ [crc:0]
Шаг 2. Цикл для каждого фрагмента 128-бит
 [A2:A1] := c1•(x191 mod P) ⊕ c2•(x127 mod P)
 [c2:c1] := [A2:A1] ⊕ [M2:M1]
Шаг 3. Сдвиги,
 [c2:c1] := с1•(x127 mod P) ⊕ с2•(x63 mod P)
Шаг 4. Редуцирование 
    [T2:T1] := c1•Ur
    [c2:c1] := T1•Pr ⊕ [c2:c1]
Возвращает: crc = c2
Таблица параметров для алгоритма:
CRC-64/XZ poly=0x42F0E1EBA9EA3693 отраженный порядок
init=0xFFFFFFFFFFFFFFFF refin=true refout=true xorout=0xFFFFFFFFFFFFFFFF check=0x995DC9BBDF1939FA

POLY=0xC96C5795D7870F42
Таблица умножения:
0x0000000000000000, 0x7D9BA13851336649, 0xFB374270A266CC92, 0x86ACE348F355AADB,
0x64B62BCAEBC387A1, 0x192D8AF2BAF0E1E8, 0x9F8169BA49A54B33, 0xE21AC88218962D7A,
0xC96C5795D7870F42, 0xB4F7F6AD86B4690B, 0x325B15E575E1C3D0, 0x4FC0B4DD24D2A599,
0xADDA7C5F3C4488E3, 0xD041DD676D77EEAA, 0x56ED3E2F9E224471, 0x2B769F17CF112238,

BarrettR u = x^128/P(x) Ur=0x9C3E466C172963D5 Pr=0x192D8AF2BAF0E1E85 (65-бит)
[1/p]*p(x) = 48663A84688941C50000000000000001
Test =995DC9BBDF1939FA ..ok

K5( 63) = x63  mod P = 0000000000000001
K4(127) = x127 mod P = DABE95AFC7875F40
K3(191) = x191 mod P = E05DD497CA393AE4
Алгоритм №. Векторный алгоритм CRC-64/WE прямой порядок
Шаг 1. загрузить вектор M с разворотом байт, 64x2 =128 бит
 [c2:c1] := [M2:M1] ⊕ [0:crc]
Шаг 2. Цикл для каждого фрагмента 128-бит
 [A2:A1] := c2•(x192 mod P) ⊕ c1•(x128 mod P)
 [c2:c1] := [A2:A1] ⊕ [M2:M1]
Шаг 3. Сдвиги,
 [c2:c1] := с2•(x128 mod P) ⊕ с1•(x64 mod P)
Шаг 4. Редуцирование 
    [T2:T1] := c2•Ux
    [c2:c1] := T2•Px ⊕ [c2:c1]
Возвращает: crc = c1
Параметры алгоритма:
CRC-64/WE poly=0x42F0E1EBA9EA3693
init=0xFFFFFFFFFFFFFFFF refin=false refout=false xorout=0xFFFFFFFFFFFFFFFF check=0x62EC59E3F1A4F00A

POLY=0x42F0E1EBA9EA3693
Таблица умножения:
0x0000000000000000, 0x42F0E1EBA9EA3693, 0x85E1C3D753D46D26, 0xC711223CFA3E5BB5,
0x493366450E42ECDF, 0x0BC387AEA7A8DA4C, 0xCCD2A5925D9681F9, 0x8E224479F47CB76A,
0x9266CC8A1C85D9BE, 0xD0962D61B56FEF2D, 0x17870F5D4F51B498, 0x5577EEB6E6BB820B,
0xDB55AACF12C73561, 0x99A54B24BB2D03F2, 0x5EB4691841135847, 0x1C4488F3E8F96ED4,

Barrett u = x^128/P(x) Ux = 0x1578D29D06CC4F872 Px = 0x142F0E1EBA9EA3693
[1/p]*p(x) = 0x1000000000000000005F5C3C7EB52FAB6
Test =62EC59E3F1A4F00A ..ok

K5( 64) = 42F0E1EBA9EA3693
K4(128) = 05F5C3C7EB52FAB6
K3(192) = 4EB938A7D257740E

Вычисление CRC по обному биту

Не знаю какой в этом смылсл, если мы только что представили вычисление CRC по 128 бит за цикл. НО хочется вернуться к истокам. Рассмотрим под микроскопом каждый отдельный бит. Как выглядит операция CRC над битами. Прежде всего можно ввести операцию сдвиг-редуцирование, как основа для построения более сложных выражений:
XT(с) := MSB(c)==1? (с<<1)^Px: (c<<1);
Такая же операция для отраженного полинома выглядит так:
XT(с) := LSB(c)==1? (с>>1)^Pr: (c>>1);
У операции есть свойство группы, которое позволяет выполнять оптимизации:
XT(a ⊕ b) == XT(a) ⊕ XT(b)
Опрерация обновления CRC может быть представлена в виде:
Cn := XT(Cn-1 ⊕ dn<<(N-1)), где N -разрядность CRC
Можно упростить, раскрыв скобки
Cn := XT(Cn-1) ⊕ dn•Px, где N -разрядность CRC
Для двух бит получаем выражение:
Cn+1 := X2•Cn-1 ⊕ [dn+1,dn]•Px
Можно продолжить, так выведен алгоритм с таблицей умножения 4 бит.
Cn+3 := X4•Cn-1 ⊕ [dn+3..dn]•Px
ИЛИ
CRC = (CRC<<4) ^ PMUL(CRC>>(N-4) ^ Data, Px)
где PMUL -- операция полиномиального умножения и редуцирования. Умножение на полином Px - это по сути - сдвиг.
Cn+3 := X4•Cn-1 ⊕ [dn+3..dn]•XN
[05.01.2020]
Я увидел возможность для дальнейшей оптимизации алгоритма расчета CRC за счет увеличения разрядности редуцирования до 64бит. При этом можно сократить число сдвиговых операций(операций умножения) при редуцировании. Для этого необходимо пересчитать константы Барретта с увеличением разрядности до 64.
Алгоритм №. Векторный алгоритм CRC-32 отраженный порядок
Шаг 1. начальное значение
 [с4:c3:c2:c1] := [0:0:0:crc]
Шаг 2. Цикл для каждого фрагмента 128-бит, кроме последнего
 [с4:c3:c2:c1] := [с4:c3:c2:c1] ⊕ [M4:M3:M2:M1]
 [с4:c3:c2:c1] := [c2:c1]•(x31+128 mod P) ⊕ [c4:c3]•(x31+64 mod P)
Шаг 3. Загрузить последний фрагмент
 [с4:c3:c2:c1] := [с4:c3:c2:c1] ⊕ [M4:M3:M2:M1]
    [c3:c2:c1] := [c2:с1]•(x31+64 mod P) ⊕ [с4:c3]•(x31 mod P)
Шаг 4. Редуцирование 
 [T4:T3:T2:T1] := [c2:c1]•Ur
    [c3:c2:c1] := [T2:T1]•Pr ⊕ [c3:c2:c1]
Возвращает: crc = c3
В данном алгоритме используется модифицированная константа Барретт'а (Ur), рассчитанная для разрядности 64 бита. Привожу константы для алгоритмов с отраженным полиномом.

CRC-32B Zlib CRC-32/ISO-HDLC
poly=0x04C11DB7 refin=true refout=true check=0xCBF43926
Ur = 0xB4E5B025F7011641 Pr = 0x1DB710641

CRC-32C (Castagnoli)
poly=0x1EDC6F41 refin=true refout=true check=0xE3069283
Ur = 0x4869EC38DEA713F1 Pr = 0x105EC76F1

CRC-32K (Koopman) BACnet
poly=0x741B8CD7 init=0xffffffff refin=true refout=true xorout=0xffffffff check=0x2D3DD0AE
Ur = 0xC25DD01C17D232CD Pr = 0x1D663B05D
Те же идеи удалось перенести на CRC-16 и CRC-8 с отраженным порядком бит.
Алгоритм №. Векторный алгоритм CRC-16 отраженный порядок
Шаг 1. начальное значение
 [с8...c1] := [0:0:0:0:0:0:0:crc]
Шаг 2. Цикл для каждого фрагмента 128-бит, кроме последнего
 [с8...c1] := [с8...c1] ⊕ [M8...M1]
 [с8...c1] := [c4...c1]•(x15+128 mod P) ⊕ [c8...c5]•(x15+64 mod P)
Шаг 3. Загрузить последний фрагмент
 [с8...c1] := [с8...c1] ⊕ [M8...M1]
 [c5...c1] := [c4...с1]•(x15+64 mod P) ⊕ [с8...c5]•(x15 mod P)
Шаг 4. Редуцирование 
 [T8...T1] := [c4...c1]•Ur
 [c5...c1] := [T4...T1]•Pr ⊕ [c5...c1]
Возвращает: crc = c5
Алгоритм №. Векторный алгоритм CRC-8 отраженный порядок
Шаг 1. начальное значение
 [с16...c1] := [0...0:crc]
Шаг 2. Цикл для каждого фрагмента 128-бит, кроме последнего
 [с16...c1] := [с16...c1] ⊕ [M16...M1]
 [с16...c1] := [c8...c1]•(x7+128 mod P) ⊕ [c16...c9]•(x7+64 mod P)
Шаг 3. Загрузить последний фрагмент
 [с16...c1] := [с16...c1] ⊕ [M16...M1]
  [c9...c1] := [c8...с1]•(x7+64 mod P) ⊕ [с16...c9]•(x7 mod P)
Шаг 4. Редуцирование 
 [T16...T1] := [c8...c1]•Ur
  [c9...c1] := [T8...T1]•Pr ⊕ [c9...c1]
Возвращает: crc = c9
Параметры алгоритма CRC-8 с отражением
CRC-8/BAC POLY=0x81 refin=true refout=true
Ur = 0x808182878899AAFF(64 бита) Pr=0x103 (9 бит)
X7 mod P = 0x01
X71 mod P = 0x81
X135 mod P = 0xC1
Алгоритм №. Векторный алгоритм CRC-64 отраженный порядок
Шаг 1. начальное значение
 [с2:c1] := [0:crc]
Шаг 2. Цикл для каждого фрагмента 128-бит, кроме последнего
 [с2:c1] := [с2:c1] ⊕ [M2:M1]
 [с2:c1] := c1•K1 ⊕ c2•K2
Шаг 3. Загрузить последний фрагмент
 [с2:c1] := [с2:c1] ⊕ [M2:M1]
 [c2:c1] := с1•K3 ⊕ с2•K4
Шаг 4. Редуцирование 
 [T2:T1] := c1•Ur
 [c2:c1] := T1•Pr ⊕ [c2:c1]
Возвращает: crc = c2
Алгоритм CRC-8, CRC-16, CRC-32 в точности такие же как и для разрядности CRC-64, отличие только в значениях сдвиговых констант Kn и {Ur,Pr}, и число разрядов выходного значения (N). Таким образом алгоритм получается универсальный для всех вариантов CRC с отраженным порядком. НО это еще не финал. Дело в том, что представленный алгоритм работает только с входными данными, длиной кратной 128бит. А мы хотим сделать универсальный алгоритм, чтобы он работал с любыми длинами. Идея такая: можно на шаге 3. загружать последний фрагмент не кратный 128 бит и выполнять соответствующий его длине сдвиг. Для этого вносим соответствующее изменение в алгоритм и просчитываем соответствующие сдвиги, получаем таблицы K3 и K4 из 16 сдвиговых констант.
Алгоритм №. Векторный алгоритм CRC-64 отраженный порядок
Шаг 1. начальное значение
 [с2:c1] := [0:crc]
Шаг 2. Цикл для каждого фрагмента 128-бит, кроме последнего
 [с2:c1] := [с2:c1] ⊕ [M2:M1]
 [с2:c1] := c1•K1 ⊕ c2•K2
Шаг 3. Загрузить последний фрагмент k = (length & 15)
 [с2:c1] := [с2:c1] ⊕ [M2:M1]
 [c2:c1] := с1•K3(k) ⊕ с2•K4(k)
Шаг 4. Редуцирование 
 [T2:T1] := c1•Ur
 [c2:c1] := T1•Pr ⊕ [c2:c1]
Возвращает: crc = c2
Вот теперь моя работа закончена. Я изобрел свой собственный способ расчета CRC! Самый быстрый, самый крутой!
Наборы констант для CRC-16/X-25 BACnet с отраженным порядком
poly=0x1021 init=0xffff refin=true refout=true xorout=0xffff check=0x906e

K1, K2:{0x8E10, 0x81BF} // x^{143}, x^{79}
K3, K4:
[ 1] = {0x4EA8, 0x97B7},// x^{-41},x^{-105}
[ 2] = {0x290C, 0xC1A3},// x^{-33},x^{-97}
[ 3] = {0xCA45, 0x9750},// x^{-25},x^{-89}
[ 4] = {0x1563, 0x5212},// x^{-17},x^{-81}
[ 5] = {0x5188, 0x33C1},// x^{-9}, x^{-73}
[ 6] = {0x0811, 0xD7B6},// x^{-1}, x^{-65}
[ 7] = {0x0100, 0xD06A},// x^{ 7}, x^{-57}
[ 8] = {0x0001, 0xCC8C},// x^{15}, x^{-49}
[ 9] = {0x1189, 0x4EA8},// x^{23}, x^{-41}
[10] = {0x19D8, 0x290C},// x^{31}, x^{-33}
[11] = {0x5ADC, 0xCA45},// x^{39}, x^{-25}
[12] = {0x1CBB, 0x1563},// x^{47}, x^{-17}
[13] = {0x0B44, 0x5188},// x^{55}, x^{-9}
[14] = {0x042B, 0x0811},// x^{63}, x^{-1}
[15] = {0x9FD5, 0x0100},// x^{71}, x^{ 7}
[ 0] = {0x81BF, 0x0001},// x^{79}, x^{15}
Ur=0x859B040B1C581911 Pr=0x10811
[24.04.2021]

Нашел способ использовать один алгоритм с разрядностью 64 бита для расчета CRC-16/CRC-24/CRC-32/CRC-64. Фактически можно использовать алгоритм CRC-64 без изменения со своим набором констант для расчета контрольных сумм меньшей разрядности, например CRC-16. При этом используется правило:

 C = A mod P => C•K = A•K mod P•K

В нашем случае К - это свдиг влево на (64-16=48) бит.

Алгоритм №. Векторный алгоритм CRC-64 прямой порядок
Шаг 1. начальное значение
 [с2:c1] := [crc<<(64-N):0]
Шаг 2. Цикл для каждого фрагмента 128-бит, кроме последнего
 [с2:c1] := [с2:c1] ⊕ [M2:M1]
 [с2:c1] := c2•K2 ⊕ c1•K1
Шаг 3. Загрузить последний фрагмент k = (length & 15)
 [с2:c1] := [с2:c1] ⊕ [M2:M1]
 [c2:c1] := с2•K4(k) ⊕ с1•K3(k)
Шаг 4. Редуцирование 
 [T2:T1] := c2•Ux ⊕ [c2:c1]
 [c2:c1] := T2•Px ⊕ [c2:c1]
Возвращает: crc = c1>>(64-N)
Наборы констант для CRC-16
poly=0x1021 init=0xffff refin=false refout=false xorout=0xffff check=0xd64e

K1, K2:{0xAEFC, 0x650B},// x^{128}, x^{192}
K3, K4:
[ 1] = {0xCBF3000000000000, 0x2AE4000000000000},// x^{-104}, x^{-40}
[ 2] = {0x9B27000000000000, 0x6128000000000000},// x^{-96}, x^{-32}
[ 3] = {0x15D2000000000000, 0x5487000000000000},// x^{-88}, x^{-24}
[ 4] = {0x9094000000000000, 0x9D71000000000000},// x^{-80}, x^{-16}
[ 5] = {0x17B9000000000000, 0x2314000000000000},// x^{-72}, x^{-8}
[ 6] = {0xDBD6000000000000, 0x0001000000000000},// x^{-64}, x^{ 0}
[ 7] = {0xAC16000000000000, 0x0100000000000000},// x^{-56}, x^{ 8}
[ 8] = {0x6266000000000000, 0x1021000000000000},// x^{-48}, x^{16}
[ 9] = {0x2AE4000000000000, 0x3331000000000000},// x^{-40}, x^{24}
[10] = {0x6128000000000000, 0x3730000000000000},// x^{-32}, x^{32}
[11] = {0x5487000000000000, 0x76B4000000000000},// x^{-24}, x^{40}
[12] = {0x9D71000000000000, 0xAA51000000000000},// x^{-16}, x^{48}
[13] = {0x2314000000000000, 0x45A0000000000000},// x^{-8}, x^{56}
[14] = {0x0001000000000000, 0xB861000000000000},// x^{ 0}, x^{64}
[15] = {0x0100000000000000, 0x47D3000000000000},// x^{ 8}, x^{72}
[ 0] = {0x1021000000000000, 0xEB23000000000000},// x^{16}, x^{80}
Ur=0x11303471A041B343 Pr=0x1021000000000000
Наборы констант для CRC-24/OpenPGP
poly=0x864cfb init=0xb704ce refin=false refout=false xorout=0x000000 check=0x21cf02

K1, K2:{0x6243DA, 0xB22B31},// x^{128}, x^{192}
K3, K4:
[ 1] = {0x2471670000000000, 0x7190920000000000},// x^{-96}, x^{-32}
[ 2] = {0xEE50080000000000, 0xC6E2490000000000},// x^{-88}, x^{-24}
[ 3] = {0xCE8F4A0000000000, 0xD19E9A0000000000},// x^{-80}, x^{-16}
[ 4] = {0x1D1CA30000000000, 0xF77C040000000000},// x^{-72}, x^{-8}
[ 5] = {0x6DC6AA0000000000, 0x0000010000000000},// x^{-64}, x^{ 0}
[ 6] = {0x67F3180000000000, 0x0001000000000000},// x^{-56}, x^{ 8}
[ 7] = {0x79152C0000000000, 0x0100000000000000},// x^{-48}, x^{16}
[ 8] = {0xE2DD700000000000, 0x864CFB0000000000},// x^{-40}, x^{24}
[ 9] = {0x7190920000000000, 0x668F480000000000},// x^{-32}, x^{32}
[10] = {0xC6E2490000000000, 0x8309D70000000000},// x^{-24}, x^{40}
[11] = {0xD19E9A0000000000, 0x3609520000000000},// x^{-16}, x^{48}
[12] = {0xF77C040000000000, 0xD9FE8C0000000000},// x^{-8}, x^{56}
[13] = {0x0000010000000000, 0x36EB3D0000000000},// x^{ 0}, x^{64}
[14] = {0x0001000000000000, 0x3B918C0000000000},// x^{ 8}, x^{72}
[15] = {0x0100000000000000, 0xF50BAF0000000000},// x^{16}, x^{80}
[ 0] = {0x864CFB0000000000, 0xFD7E0C0000000000},// x^{24}, x^{88}
Ur=0xF845FE2493242DA4 Pr=0x864CFB0000000000
Наборы констант для CRC-32/BZIP2
poly=0x04c11db7 init=0xffffffff refin=false refout=false xorout=0xffffffff check=0xfc891918

K1, K2:{0xE8A45605, 0xC5B9CD4C},// x^{128}, x^{192}
K3, K4:
[ 1] = {0x052B9A0400000000, 0x876D81F800000000},// x^{-88}, x^{-24}
[ 2] = {0x3C5F6F6B00000000, 0x1ACA48EB00000000},// x^{-80}, x^{-16}
[ 3] = {0xBE519DF400000000, 0xA9D3E6A600000000},// x^{-72}, x^{-8}
[ 4] = {0xD02DD97400000000, 0x0000000100000000},// x^{-64}, x^{ 0}
[ 5] = {0x3C423FE900000000, 0x0000010000000000},// x^{-56}, x^{ 8}
[ 6] = {0xA3011FF400000000, 0x0001000000000000},// x^{-48}, x^{16}
[ 7] = {0xFD7384D700000000, 0x0100000000000000},// x^{-40}, x^{24}
[ 8] = {0xCBF1ACDA00000000, 0x04C11DB700000000},// x^{-32}, x^{32}
[ 9] = {0x876D81F800000000, 0xD219C1DC00000000},// x^{-24}, x^{40}
[10] = {0x1ACA48EB00000000, 0x01D8AC8700000000},// x^{-16}, x^{48}
[11] = {0xA9D3E6A600000000, 0xDC6D9AB700000000},// x^{-8}, x^{56}
[12] = {0x0000000100000000, 0x490D678D00000000},// x^{ 0}, x^{64}
[13] = {0x0000010000000000, 0x1B280D7800000000},// x^{ 8}, x^{72}
[14] = {0x0001000000000000, 0x4F57681100000000},// x^{16}, x^{80}
[15] = {0x0100000000000000, 0x5BA1DCCA00000000},// x^{24}, x^{88}
[ 0] = {0x04C11DB700000000, 0xF200AA6600000000},// x^{32}, x^{96}
Ur=0x04D101DF481B4E5A Pr=0x04C11DB700000000