Попалась мне недавно статья Синус, косинус, квадратный корень FixedPoint. Автор размышляет как можно не затратно рассчитывать координаты и углы в микроконтроллере. Попробовал я подсказать автору пару аппроксимаций, но он оказался общителен только на тему "упадка автоматизации в РФ", а по делу как то не сложился диалог. Посмотрел, такие статьи с такой тематикой не редкость. Например, очень хорошая статья Как посчитать синус быстрее всех на Xабре. В общем разгрузил себе голову на майских праздниках от главного хобби - геометрической алгебры.
В процессе изучения всего этого, возник у меня вопрос - а зачем вообще нужно аппроксимировать sin,cos, arctan и еще и в привязке к числу в двоичной системе, если есть декартовы координаты?
Из ответа на этот вопрос родилась идея этой статьи. Будет длинно, но если на примере подробно разбираться с работой машинного эпсилон и автоматическим дифференцированием, короче не получится. Следите за мыслью по ходу изложения. Начну с главного тезиса, и разверну по шагам как это работает на примере операций с единичной окружностью. На математический трактат не замахиваюсь, поэтому пишу только основное и не строго.
Автоматическим дифференцированием можно назвать любую конечную разность, например dy=(y(x+ε)-y(x-ε))/(2*ε). Разность здесь взята центральная, так как она дает меньшую погрешность. Параметр ε это машинный ноль. За счет округления до младшего бита его главное свойство: ε^2=0.
Эта статья по сути не более, чем описание основных моментов идеи. И если у кого то появится желание поставить эту идею на строгие математические рельсы, с удовольствием готов поучаствовать. Кто в этом случае опубликует финальную версию мне искренне не важно.
Дискретные синус и косинус как координаты на единичной окружности
(1) Зададим координаты на окружности по школьным формулам.
Преобразуем их для случая единичной окружности r=1.
Если y определяем через x, то сделать это можем только с точностью до отражения относительно оси {x}. В общем такое задание координат задает окружность в виде двух ветвей.
(2) Теперь это надо как то привязать это к числу в двоичной системе: x`=1+i*ε.
Привяжусь к числу fixed point (мантиссе числа в двоичной системе). i - порядковый номер числа.
Координаты x,y принадлежат отрезкам [-1,1]. Преобразуем данный отрезок в формат мантиссы binary18 (как я понял, автор первой упомянутой статьи хотел ограничиться 18-ю битами, чтобы экономить место в микроконтроллере). А мне так понравилось сделать, чтобы оставить в резерве часть битов мантиссы float32, что полезно для исключения ошибок типа overflow и underflow. Хорошая статья про форматы чисел, вдруг кому будет полезным.
В общем далее по умолчанию числа округляются до искусственного машинного ноля ε=2^(-18) операцией: round(x*ε^(-1))*ε.
Мантисса принадлежит отрезку [1,2). Квадратная скобка означает, что граница включена в диапазон. Круглая скобка означает, что граница в диапазон не включена. Пусть число мантиссы 18 бит, тогда число имеет точность ε~5*10^(-6) или в двоичной записи ε=2^(-18).
(3) Для начала найдем коэффициенты преобразования числа из диапазона
x[-1,1] --> x`[1,2). И учтем, что начало отсчета углов в математике обычно от точки x=x`=1.
Тогда можем записать линейную систему уравнений и получить коэффициенты преобразования:
Подставим вычисленные коэффициенты в выражение для x` и сразу выведем обратное преобразование:
(4) Квадраты чисел c упрощением за счет свойства ε^2=0.
(5) Мы вводили, что r^2=x^+y^2=1. Подставляя выражения для x,y получаем это выражение:
Решая уравнение относительно y`, и снова используя ε^2=0 получаем уравнение для сдвинутой окружности.
(6) Количество чисел в каждой координате 2^18=262144
(7) Итак у функции y две ветви. Округления я применил в ключевых местах, где это сделает компьютер. При этом я рассматриваю неплохой случай, в котором программист пользовался округлением до ближайшего. С округлением до меньшего, погрешности, от рассматриваемых здесь, возрастут в 100 раз. Поэтому использую всюду "round".
Функции верхней и нижней ветви и их график:
(8) Оценим погрешности округления. График относительной погрешности приведен ниже.
Погрешности примерно от x=1.02 до x=1.98 в пределах исходных 5*10^(-6). Дальше погрешность возрастает до 10^(-5). В пике погрешность возрастает до 2*10^(-5). Итого погрешность к краям окружности вырастает примерно в 10 раз.
Лекарство от увеличения погрешностей простое:
При углах единичного вектора, (от -пи/4 до +пи/4)~(1.85<x<2) и (от 3пи/4 до 5пи/4)~ (1<x<1.15), выражать x`=f(y`)
Во всем остальном диапазоне выражать y`=f(x`).
Ну или как я можно данный эффект проигнорировать, чтобы не плодить математику.
(9) Теперь зная с заданной точностью x,y на единичной окружности можно расcчитать x,y, которые определяют дискретные cos,sin.
Проверим как будет влиять устранение ε из знаменателя. Оно повышает погрешность в среднем в 2 раза. Упрощать себе жизнь или нет на ваш выбор. Я не упрощаю, и как исключение ε влияет на дальнейшее не проверял.
(10) Ремарка. Переходя к истории про определение углов, надо уточнить следующее.
Косинус и синус кратных углов задают полиномы Чебышева. Формулы из Википедии:
Зная значение cos(φ),sin(φ) для какого то малого угла для перекрытия окружности с кратностью n, по этим формулам не сложно расcчитать cos(nφ),sin(nφ). Множественное возведение в степень пусть вас не пугает. Так как есть хорошая статья на тему быстрого возведения в целую степень, ну и на Википедии на эту тему много есть. Написал для альтернативы, не все знают про эти ряды.
Определение углов через координаты
(11) Подумайте перед тем как пользоваться обратными тригонометрическими функциями, нужны ли они вам вообще. С ними не все так просто. Ряды для обратных тригонометрических функций сходятся медленно, или если простым языком, требуют МНОГО членов ряда для точного вычисления. Для реализации на микроконтроллерах это не вариант. Можно определить угол через длину дуги, что как увидите проще, и позволяет проиллюстрировать ограничения для применения операций с машинным нулем. До этого момента кому то могло показаться, что данный прием панацея..
Штрих так же как и ранее обозначает преобразованную величину, а не производную. Интегрировать будем верхнюю ветвь окружности.
Для единичной окружности угол равен длине дуги окружности. i - как и ранее порядковый номер координаты на окружности.
(12) Дифференциал y определим центральной разностью, так как нецентральная дает большую погрешность.
Математически данная производная - это дискретный котангенс, так как она равна отношению производной синуса к производной косинуса. На графике действительно кривая котангенса, ноль находится между точками i=262142/2 и i=262144/2, именно там, где находится угол пи/2. Произведение с дискретным тангенсом (y/x) добавил на график, чтобы видеть, что их произведение действительно =1.
(13) Безболезненно можно сокращать только элементы типа iε^2, а i^2ε^2 уже нельзя, почему думаю понятно. Возведем cot.d в квадрат и упростим.
Сравним, что получилось, сравнив погрешность для полной формы и сокращенной. Выясняется, что чем ближе к пи/2 (в самой точке пи/2 погрешность 400%), тем метод хуже работает, аналогично тому как это происходило по краям аппроксимации sin. Чем функция ближе к нулю, тем хуже работает такой способ дифференцирования. Ну и, конечно, оно от типа функции зависит.
(14) Посмотрим сколько погрешности накопится при вычислении угла пи. Во-первых видим, что число получилось комплексным, что требует при программировании исключать точки с мнимым результатом. Во-вторых, казалось бы, накопилась не маленькая погрешность определения угла (4-3,14)/3.14=+27.5%
Но посмотрим, что будет если рассчитывать угол для всех точек окружности. Считать 262144 суммы из 262144 элементов, долго даже для мощных компьютеров, поэтому в этой части снова перейду от суммы к интегралу, на суть это не особо влияет.
Оказывается, что угол почти линейно зависит от номера точки, то есть, по сути, имеет место почти линейное отображение множества углов на множество номеров. Дело в каком-то почти постоянном множителе, для каждого дифференциального компонента в сумме.
Вторая функция на графике это прямая для угла от 0 до пи в зависимости от номера.
(15) Если разделить первую функцию на вторую, зависимость оказывается почти линейна, изменения линейного коэффициента при углах > пи/64 происходят в пределах 0,1%. При малых углах погрешность очень большая. Но если такая точность определения угла устраивает, то на этом все. Если не устраивает, то более точно можно, например, методом Ньютона посчитать, с начальным приближением из описанного подхода.
Но можно немного по другому. "1.27326" это коэффициент пропорциональности прямой φ(i) и интегрального угла. Нелинейный множитель пропорционален (пи-φ/2)/(1-cos(φ/2)), что можно использовать, для развития описанного метода.
Чтобы ликвидировать этот нюанс, я воспользовался методом наименьших квадратов в виде решения основанного на псевдо-обратной матрице. После применения интерполяции полиномом 2-й степени от аргумента (пи-φ/2)/(1-cos(φ/2)), остаточная погрешность получилась 10^(-6).
Подробно не расписываю, так как уже много написал.. Кто понимает в этой математике, думаю сразу поймут что я сделал, ну или вопросы зададут.
Итого, аппроксимация для угла с аргументом в виде номера точки оказывается полиномом второй степени от аргумента X, с нулевой константой C0. Зависимость для угла получилась вполне вычислимой даже на микроконтроллере.
Как быстро рассчитывать корни есть много методов. Не решенным остались три вопроса на перспективу:
Как эффективно использовать экспоненциальную часть float32 для расчета. 10^(-10) это за пределами точности представления числа только 18-ти битной мантиссой.
Как эффективно рассчитывать дроби.
Все таки потребовалось и косинусы вычислять от дискретных углов.
Выводы
Описанное выше означает, что зависимость углов от номера на окружности как то можно заложить в логику микроконтроллера, и относительно точно аппроксимировать углы без интегралов и обратных тригонометрических функций.
Нюанс есть всегда, в данном случае он в том, что такая математика потребует коэффициентов для калибровки в контроллере под каждый конкретный метод реализации. Но на мой взгляд это небольшая цена. Когда закончу с геометрической алгеброй, думаю обобщить этот метод на три измерения. Спойлер - там свойства матричного машинного эпсилон еще интереснее, чем здесь. И так же будет интересно зашить в число float32 три мнимых единицы, это так же одна из причин написать здесь про число 18 бит, чтобы оставить свободными необходимые три бита и 2 бита оставить в резерве.