Библиотека математических операций с фиксированной точкой

[menzoda]

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

Пока эти функции реализованы только для архитектуры ARMv7-M (Cortex-M3). В принципе, код можно компилировать и для ARMv7E-M (Cortex-M4), но это будет не совсем оптимально, так как под Cortex-M4 некоторые моменты можно реализовать более эффективно. В дальнейшем хочу добавить еще несколько функций, сейчас вот приступлю к делению и квадратному корню. Так же собираюсь добавить реализацию для других архитектур, может быть даже не ARM если кто-нибудь заинтересуется и захочет помочь.

Циферка 4 в названии означает, что библиотека предполагает использование 32-битных (4 байта) целых для представления чисел с фиксированной точкой. Библиотека поддерживает работу с числами в формате от Q1 до Q31, т.е. каждая функция (за некоторыми исключениями) имеет свою версию для каждого конкретного формата.

Использование 32-битных чисел приводит к необходимости вычисления их 64-битного произведения, т.е. очень желательно, чтобы целевое устройство имело такой аппаратный умножитель, иначе это все надо будет реализовывать ручками, и будет это все работать очень медленно. Поэтому я хочу сделать вторую версию библиотеки, работающую с 16-битными числами. В большинстве случаев их точности будет достаточно, а вычислительных мощностей нужно будет гораздо меньше. Для начала думаю реализовать ее для Cortex-M0, в нем же, если не ошибаюсь, как раз нету умножителя с 64-битным результатом.

Ах да, для ARM это все написано на ARM-овском ассемблере, который armasm, который в Keil. К библиотеке нужно будет добавить какой-нибудь make чтобы её собирать (прям сейчас займусь), кое-то по мелочи допилить, прокомментировать и можно будет уже использовать. Пока же вы можете просто подсмотреть что и как реализовано (если интересно), или выдернуть себе в проект нужную функцию, они все проверены и должны работать. По точности и скорости каждой функции спрашивайте тут, у меня, GitHub вики со всей этой информацией пока еще в планах.

Собственно проект на GitHub

Итого на текущий момент реализовано следующее:

ARMv7-M (Cortex-M3)

clip (ограничивает величину верхним и нижним пределом, чтобы не писать несколько if)
int (возвращает целую часть)
frac (возвращает дробную часть)
round (округляет до ближайшего целого)
sgn (вычисляет функцию signum)
abs (возвращает абсолютное значение)
mul (умножает два числа)
sin (вычисляет синус)
cos (соответственно косинус)

[golf2109]

полезно было бы добавить arctg & sqrt(a`2+b~2) для дискретного фурье

[menzoda]

Хорошо, приму к сведению. Собственно я сам после деления и простого корня хотел еще чего-нибудь добавить, вроде тех же дополнительных тригонометрических функций. К сожалению, это не быстро делается. Вон с теми же простейшими синусами я очень долго возился. Сначала сам алгоритм подбирал, понимал, переваривал, смотрел как у других сделано. Затем пытался наиболее оптимально его реализовать, да чтобы погрешность была минимальная. Куча итераций было: вроде уже всё выжал, а потом бац! Приходит идея в голову, которая еще на пару тактов ускоряет функцию.

[drakon]

menzoda, изобретаете велосипед. Я на cortex M3,4 использую IQmath от TI. Она написана на asm, хорошо документирована и хорошо оформлена. У меня есть скомпилированные библиотеки + исходник для M4.

[menzoda]

Насчет велосипеда то поосторожнее.

Во-первых, там нету некоторых функций, которые уже есть у меня, и которые я собираюсь реализовать.

Во-вторых, я ради примера посмотрел на их реализацию sin (PU). В документах написано 50 тактов, в исходниках аж до 80 тактов. Не знаю где правда, но моя реализация работает максимум за 42 такта, плюс еще два (вроде), если считать выход из функции. Правда, тут надо уточнить как они измеряли свое время, это можно по разному делать, тот же выход из функции можно считать а можно и нет. Посмотрю их исходники поподробнее, тогда можно будет сказать точно, но скорее всего все верно, и моя реализация быстрее.

В-третьих, они используют таблицу на 512 значений. 512 значений, Карл! Каждое по четыре байта, всего 2 Кб памяти. У меня же используется всего 64 значения, а вся таблица занимает 256 байт. Чувствуешь разницу? А если бы я взял такую таблицу как у них, то тогда бы функция работала не за 40 тактов, а за 20-30.

В-четвертых, я хочу сделать уменьшенную версию библиотеки для Cortex-M0 и может каких-нибудь восьмибитников. Такой у Техаса точно нету.

[drakon]

У меня оптимизация проектов в основном по скорости, 2к это совсем ничего.
Иногда приходится с фрискейлами M0+ работать, но там ПЗУ хватает, только с ОЗУ напряг.
Если будет версия для M0 то Ваша библиотека имеет смысл для меня
Если не закрытый проект то выкладывайте результат, будем тестировать.
У Вас оформление функций свое или как у техаса?

[menzoda]

Дело в том, что таблицы и код функций обычно кладут в RAM, у которого 0 wait states, потому что flash на всех микроконтроллерах очень медленная. Я как-то интересовался, так вот у 90% МК на полной тактовой частоте доступ к flash это примерно 3-5 тактов. У некоторых редких экземпляров есть "акселераторы" (например ART у ST), которые ускоряют доступ к памяти, но 1 такт тоже не гарантируют. Так что маленький размер кода и таблиц вполне актуален. Хотя, можно оставить таблицы во flash, а код перенести в RAM. Возможно, при использовании больших таблиц во flash увеличение времени доступа к ним нивелируется ускорением самих функций. В общем надо подумать.

Проект не закрытый. Ссылка на GitHub в первом сообщении. Оформление функций свое. Как я говорил выше все имеющиеся на данный момент функции мной проверены и работают. Если вы их дополнительно где-нибудь протестируете, то будет совсем хорошо. Однако проект в целом пока еще не готов. Я тут застрял на создании makefile, смотрю что люди используют, cmake, make, waf, scons. Вариантов много, все корявые, а я перфекционист... А еще надо вменяемую справку сделать, чтобы со временем работы каждой функции, ее точностью... В общем, сейчас можно просто взять код и вставить в свой проект. Хотя бы просто чтобы я мог получить какой-нибудь полезный feedback.

[drakon]

У взрослых МК и FPU есть, хоть он и медленнее. Я iq математику применяю в МК с частотой 24 МГц.
Удобно когда есть:

Код: Выделить всё

#define GLOBAL_Q                24

а потом:

Код: Выделить всё

#if GLOBAL_Q == 29
#define _IQsin(A)               _IQ29sin(A)
#endif
#if GLOBAL_Q == 28
#define _IQsin(A)               _IQ28sin(A)
#endif
#if GLOBAL_Q == 27
#define _IQsin(A)               _IQ27sin(A)

Можно сразу писать:

Код: Выделить всё

a = _IQsin(b);

Еще бы имена функций или обертки к ним чтоб заменить техасовские Вашими (заменить только src).
Хотелось бы увидеть конфигурирование: куда класть таблицу для каждой ф-и 0) ОЗУ 1) ПЗУ.
cmake и прочее удобно (наверное), но, а если я захочу скомпилировать под IAR? Я так понимаю написано в GCC asm?

[menzoda]

У взрослых МК и FPU есть, хоть он и медленнее. Я iq математику применяю в МК с частотой 24 МГц.

FPU это уже другая тема. Кроме того на данный момент большинство МК (даже с частотой более 100 МГц) его не имеют.

Удобно когда есть GLOBAL_Q

Такую функцию реализовывать пока не собираюсь, потому что работать с одним общим масштабом считаю неудобным (неоптимальным): где-то нужно больше целых разрядов, где-то дробных. Где-то это нужно ради оптимизации, чтобы попусту не приводить промежуточные результаты к одному масштабу. К тому же - Explicit is better than implicit.

Еще бы имена функций или обертки к ним чтоб заменить техасовские Вашими (заменить только src).

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

Хотелось бы увидеть конфигурирование: куда класть таблицу для каждой ф-и 0) ОЗУ 1) ПЗУ.

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

Я так понимаю написано в GCC asm?

Нет, с использованием ARM Compiler toolchain, который по умолчанию идет вместе с Keil MDK-ARM. Мне с ним удобнее работать.

а если я захочу скомпилировать под IAR?

Под IAR скомпилировать не получится, но если IAR использует тот же Application Binary Interface и Procedure Call Standard, то он должен слинковаться и работать с библиотекой скомпилированной в Keil. GNU ARM toolchain это поддерживает, правда я не пробовал, нужно проверить.

[Immortal_Buka]

есть же в CMSIS DSP, чем не устраивает то?

[menzoda]

1. CMSIS DSP поддерживает только Q32.31 и Q16.15 форматы.
2. CMSIS DSP почему-то ориентирована на векторные операции.
3. CMSIS DSP реализует очень маленький набор функций.
4. CMSIS DSP написана на C и существенно более медленна.
5. CMSIS DSP написана на C и занимает больше памяти.
6. Не нашел информации по точности функций CMSIS DSP, но судя по тому, что для синуса она использует таблицу и линейную аппроксимацию между точками таблицы, то как минимум точность тригонометрических функций вызывает большие сомнения.

[Alexey_N]

На гитхабе Page not found ...

[menzoda]

Я тут просто разбираюсь с самим гитхабом, сравниваю его с гитлабом. То там попробовал, то там. Походу дела там кое-что переименовал и ссылка отвалилась. Вот новая. Больше экспериментировать не буду. А вообще, в данный момент закопался с операцией деления. Используемый алгоритм дает на выходе погрешность где-то в полтора бита. Чтобы её исправить нужно будет использовать весьма дорогостоящий набор операций, а так неохота...

[Alexey_N]

А какой алгоритм используете для деления?

[menzoda]

Сначала вычисляю инверсию знаменателя (1/B). Первое приближение беру из таблицы 128x8bit или 256x8bit. Далее две итерации метода Ньютона-Рафсона. После этого умножаю на числитель A. На выходе получается частное:

Код: Выделить всё

round(A/B*2^Q) + E

где Q - степень соответствующая текущему формату с фиксированной точкой, E - погрешность в диапазоне от -1 до 1 включительно, A = a*2^Q - числитель в формате с фиксированной точкой, B = b*2^Q - знаменатель в формате с фиксированной точкой.

[menzoda]

Наверное сделаю две версии функции. Одну с упором на скорость, все-равно погрешность не такая уж большая. Вторую - с максимальной точностью. Да, пожалуй так будет лучше всего.

[Alexey_N]

А простой алгоритм (мантиссы делятся как целые числа, порядки вычитаются) дает бОльшую погрешность чем этот?

[menzoda]

У меня числа не с плавающей точкой, а с фиксированной. То есть мантиссы как бы нет, она в уме (в википедии лучше написано), работа происходит только с основанием. Так что я тоже просто делю целые числа. Собственно, в аппаратных FPU происходит тоже самое, просто там делятся не 32-битные целые, а 24-битные (в случае чисел одинарной точности). Ну и алгоритмы немного другие используются, оптимизированные для аппаратной реализации, какой-нибудь алгоритм Гольдшмидта например.

[Alexey_N]

Оу, что-то я попутал форматы, извиняюсь

[Alex_Megavolt_79]

Я тут просто разбираюсь с самим гитхабом, сравниваю его с гитлабом. То там попробовал, то там. Походу дела там кое-что переименовал и ссылка отвалилась. Вот новая. Больше экспериментировать не буду. А вообще, в данный момент закопался с операцией деления. Используемый алгоритм дает на выходе погрешность где-то в полтора бита. Чтобы её исправить нужно будет использовать весьма дорогостоящий набор операций, а так неохота...

ССыль битая. С первого поста тож не открывается.