Комментарии 17
Рекомендую юзать Cordic. Просто и надёжно.
Вы имеете ввиду готовую библиотеку, или алгоритм?
Много лет назад, когда работал в АО "Концерн "Созвездие", мы использовали данный алгоритм на FPGA, для синтеза синуса. Не моя часть, но я тоже ее касался, косвенно. Это были еще времена AT91RM9200, IXP465 и PXA270.
Здесь больше речь об оптимизации вычислений, т.к. при разложении в ряд тейлора, использовал все 32 бита, как дробную часть числа. Т.е. двоичная точка стоит прямо перед числом. А результат: 32-битное число. По этому после умножения, можно использовать только старшую часть результата, без сдвига. Но, потом, подумал, и сложил еще и младшую
Какое преимущество в данном случае даст Cordic?
Не сразу понял о чём идёт речь. Думал что используется метод неподвижной точки f(x) = x. И удивился, разве можно найти синус и косинус таким образом? Корень то легко, а вот тригономертические функции как!? Потом понял, что речь идёт про формат записи числа.
Можно подробнее про корень? Я не занимался ранее этим вопросом. Использовал готовые бибилиотеки. Отлаженные, стандартные.
Просто, был у меня в жизни период, когда я прошел так называемую "работу на работе". Это научило меня нескольким вещам.
Во-первых не доверять сторонним реалзациям, а если доверять, то покрывать их тестами, на 100% как BlackBox. Впрочем, лушче бы иметь открытую реализацию. Функции sin(x), cos(x) разложены в ряд тейлора, с корректно подобранными коэффициентами. Точности вычислений я знаю. Потенициальная ошибка может быть только работе процессора. Отчасти по этому выбран Тейлор. К тому же он FixedPoint.
Удалось достигнуть точности порядка 1LSB, при 32-битных вычислениях. Результат то все равно, фактически, 31-бит, т.к. еще есть знак.
Квадратный корень вычислял используя четыре аппроксимации. Первая,- линейная (всегда ниже оригинального значения). Вторая, апроксимация пораболы. Но пораболу следовало сместить, по этому аппроксимировал квадратично ошибку смещения пораболы по оси X. И наконец, апроксимация остаточной ошибки, которая дала точность результа +/-0.15. Последнюю, просто на глаз сделал, по графику. Поленился корни уравнения считать. А затем использовал итеративное приближение. Но это уже 7-12 бит дробной части из ее 16бит. Операцию деления не использую вообще. Ни при вычислении квадратного корня, ни при вычислении корня квадратного.
А что за метод неподвижной точки?
Как известно, каноническая форма любого уравнения это f(x) = 0. Здесь неявно предполагается, что число 0 входит в область опредения функции f(x), что в общем виде не так, более того эта функция вообще может быть определена не на числовом множестве. Тогда можно записать уравнение в другом виде f(x) = x. Это более универсальная форма, т.к. даже если нуля в области определения нет, то хоть какой то x уж точно существует)).
При этом можно заметить, что если f(x) = x, то f(f(x))=x и вообще f(f(f(...(x)...)))= x. Поэтому точка x называется неподвижной точной (fixed point), т.к. применение к ней функции f не изменяет саму x. А значит рекурсивное применение функции f самой к себе автоматически находит решение уравнения f(x)=x (либо не находит - т.к. в численных методах задача либо сходится, либо расходится))).
Как с помощью этого вычислить квадратный корень? По определению квадратный корень числа a это такое x, что x^2=a и x>=0. Запишем x^2=a в виде f(x) = x. f(x) = a/x = x. Стартуем от любой точки, например 1 и рекурсивно применяем f(x) к самой себе и... В результате выходит бесконечный цикл, в котором два значения x1=1 и x2=2 повторяются, прыгая вокруг правильного ответа. Решение расходится, печаль, бывает.
Но очевидно, что правильный ответ находится между x1 и x2, поэтому перепишем f(x) в таком виде f(x)=(x + a/x)/2. Видно, что f(sqrt(a))=sqrt(a), значит такая форма записи тоже удовлетворяет условию задачи. Попробуйте теперь рекурсивно вызывать f(x)=(x + a/x)/2, начиная с любого положительного x - сами увидите, что задача сходится очень быстро.
Я в эксельке накидал для наглядности. Уже на 4й итерации получился ответ с максимально возможной точностью:
К тому же в отличии от формулы Тейлора этот метод не накапливает ошибку вичислений. Т.к. каждый раз на вход f поступает неточное значение неподвижной точки x, и каждый раз f вычисляет её чуть точнее, тем самым автоматически подавляя ошибку.
Не уловил ход мысли. Подумаю еще.
Судя по последней формуле, Вы пришли к классическому вычислению корня квадратного, используя формулу применяемую в newlibc (__ieee754_sqrt(x)), т.е. формулу Герона. Ее еще в Египте применяли для этих целей. Дробление прямоугольника с длинной равной площади и высотой равной единице (начальное условие, в этом случае не 1.5, а 2 должно быть для корня из двух.
Однако.... При вычислении корня, таким образом, требуется операция деления. Там в, знаменателе стоит Y, а не степень двойки, а значит двоичным сдвигом уже не отделаешься. В моем случае, это деление еще и 64-битное. Деление,- достаточно дорогая операция, которая доступна не везде. В процессорах, не везде, а в FPGA и подавно.
Кроме того, Вы изначально очень удачно выбрали начальное условие, это непременно повлияло на количество итераций.
CORDIC- не не слышали?)
Да, Вы не первый кто упоминает Cordic. Обычно, это на слуху у ПЛИС-оводов, и тех кто занимается ЦОС. Сколько времени это займет на STM32?
У меня нет итеративной части при вычислении тригонометрии, она вычисляется всегда за одно и то же фиксированное время. Я понимаю, что быстрая математика на МК обычно использует таблиные методы с линейной апроксимацией. Это если говорить о ЦОС. Если говорить о моей задаче, то мне нужна определенная абсолютная точность. Это для станка CNC (дома стоит такой, фрезерный).
Если прочитать комментарии от первого до четвертого, то CORDIC прозвучит два раза. Ой, уже три.
А почему не использовать табличные значения, сгеренированные на этапе компиляции?
Так же хранить четверть диапазона и оставшиеся вычислять.
Да, я видел подобные методы. Если посмотреть готовые бибилиотеки cmsis, то можно найти быструю математику c числами в формате q31. Я смотрел их перед тем как писать код. Там, как раз таки таблица с линейной интерполяцией.
Собственно, приблизительно так и происходит в моем случае. Только я не интерполирую табличные значения. Эти значения вычисляются каждый раз, используя ряд Тейлора. Т.е. имеется интерполяция sin(x), на интервале 0..PI/4. Имеется интеполяция cos(x), на интервале 0..PI/4. Аргумент приводится к этому интервалу и нужной функции, и вычисляется точное значение для каждого аргумента. Я просчитал. Точность 32-битного числа, составляет 2.32*10^-10. Я просчитал sin(x) и cos(x) приблизительно с этой точностью на указанных интервалах. Это обычный sin(x), cos(x), но просчитанный в фиксированной точке, и смасшатбированный в диапазон 0..2^32 или [-2^32..2^32-1].
Я тоже большой фанат фиксированных вычислений, они предсказуемы и детерминированы на всех платформах. Мне показалось что в моем случай точности pi/512 за глаза поэтому я использовал табличный подход. Для корня я использовал побитовый поиск, быстр и прост. Но потом понял что нативный способ все равно быстрее и оставил так:
#[inline(always)]
fn sqrt(self) -> FixFloat {
let a = (self as f64).sqrt() as i64;
let b = a + 1;
if b * b > self { a } else { b }
}Еще вопрос как вы храните дробную часть? У вас помножено на какую-то константу?
Начальное приближение корня квадратного, в моем случае, вычисляется аналитически. Это требует 9 умножений.
Первое,- это линейная апроксимация на участке [2^n..2^n+1]. Максимальная абсолютная ошибка результата линейной аппроксимации, порядка 500 (если правильно помню). Линейная аппроксимация требует 1 умножение (обычный kx+b). Коэффициенты при этом являются константами в формате n*sqrt(2)^m или n*2^m, в зависимости от того четное или нечетное число бит аргумента поступило на вход.
Остаток от линейной аппроксимации имеет квадратичную форму. Это порабола Однако, парабола смещенная. Смещение параболы квадратичной ошибки тоже имеет квадрвтичный характер по отношению к аргументу. Корректировка аргумента, требует двух умножений. Первый это вычисление пораболы, второй,- коэффициент масштабирования корректировки. Это разница между центральным значением интервала [2^n..2^n+1], и максимумом аппроксимации. Затем еще два умножения. Это вычисление смещенной пораболы аппроксимации квадратичной ошибки, и некий коэффициент,- максимум квадратичной ошибки. Подобный подход позволил добиться абсолютной точности результата, порядка 6, для 32-битного аргумента.
Затем, аппроксимировал остаток ошибки, уже на глаз. Впринципе, можно посчитать производные, найти максимумы и посчитать точнее. Это еще увеличило точность результата до 0.15. Таким образом начальное значение квадратного корня уложилось в абсолютную погрешность 0.15.
А затем, уже использовал бинарный поиск. Вообще, результат функции sqrt от 32 битного аргумента, в случае получения целочисленного результата, имеет только 16-бит. Я вычисляю 32-бита результата. 16-бит,- целая часть. И 16-бит,- дробная. Дробная часть при начальной аппроксимации вычисляется с точностью 2-9 бит. Т.е. остается бинарным поиском найти 13-7бит.
Константы хранятся непосредственно в коде, зачастую, как пары чисел, сдвинутые на степень 2-ки.
У меня на работе лежит - "Digital Signal Processing Applicaions using the ADSP-2100 Family" от Analog Devices. Там масса информации на которой я учился писать подобные функции. Сама фраза про ряд Тейлора как бы немного уже напрягает. В книге почти как бы то же самое, но используются коэффициенты, которые ищутся на минимизации средне-квадратической ошибки на всем интервале аппроксимации.
Синус, косинус, квадратный корень FixedPoint