September 14th, 2021

QuatCore

atan(y/x) на двух умножениях!

Чего-то никак меня не отпустит эта тема, всё кажется, что есть очень простой и эффективный метод, надо только его найти!

Сейчас вот такое приближение "придумал". Напомню: у меня на входе два числа, x и y, причём x2+y2=1.
Допустим, что мы перешли к квадрантам I и IV, т.е x≥0, а потом ещё обеспечили условие x≥|y| (поменяв местами аргументы, если необходимо), т.е угол будет в диапазоне от -45 до +45 градусов. Тогда:



Максимальная ошибка составляет 0,041°, или 2,5 угловые минуты, и, возможно, ещё можно коэффициенты чуть оптимизировать.

Многовато вообще для моей задачи: мне ж этот угол ещё на 2 умножать (т.к в кватернионах угол был половинный), это будет 0,082°, при допустимой ошибке в 0,1°. Но может, кому пригодится такое приближение.

Под катом упоротые выкладки, которые привели к этому выражению... Ну как упоротые - вспоминали школьную математику, ну максимум 1-й курс института.

Collapse )

Наверное, эта тема меня так увлекла, потому что всегда интересно было узнать, а как же компьютер вообще считает сложные математические функции? Как в целом работает, складывает и умножает - уже разобрался, когда АЛУ ковырял, а вот эту математику обходили стороной как могли. В math.h не заглянешь, там одни заголовки, а функции уже откомпилированы под конкретную архитектуру. И в процессор не заглянешь (мой-исключение :), как там FPU вкалывает, а там же явно микрокода выше крыши! Кнут дальше арифметики не пошёл, во многих учебниках по программированию не мудрствуя лукаво говорят: ряд Тэйлора. А тут вроде как и "по делу" пришлось :)
QuatCore

Ай да Пафнутий Львович!

Решил ещё немного поковыряться со своим арктангенсом. Хотел применить алгоритм Ремеза, но начал с узлов Чебышёва.

И для начала со своего "линейного приближения", atan(y/x)≈ky. Всеми правдами и неправдами, решив численно нелинейное уравнение (точнее, отдав его на растерзание WolframAlpha), получил значение k=1,00120570556

А вот Чебышёвский подход: мы хотим приблизить эту функцию на углах от -5,625° до +5,625°, чему соответствует y на отрезке [-A;A] = [-0,098; 0,098], и мы просто выбираем, в каких точках хотим ТОЧНЫЙ ответ. Для линейного приближения точный ответ может достигаться в ТРЁХ точках. Но эти точки - не просто края отрезка и его середина, а величины



где n=3 - количество точек, k=1..3 - номер точки. Получаются значения y1 = A·cos(π/6)≈0,084885, y2=0, y3=-A·cos(π/6). Точное равенство в нуле получается "автоматом". А из одной из оставшихся точек найдём k, также потребовав точного равенства:





То есть, 6 значащих разрядов мы получили точно! Если помножить на 32768, получим 32807,48, что мы округлим до 32807 - это будет наше целочисленное представление. А то значение, что мы нашли в прошлый раз, при умножении на 32768 даст 32807,51 и будет округлено до 32808. Ну, почти :))

А может и в более сложной задаче нам Чебышёв поможет?

Collapse )


О как: в пол-пинка, безо всяких интегралов, улучшили наш результат, не шибко сильно, на 5%. Что же будет, если сюда ещё Ремеза натравить? Впрочем, думаю, ещё пару процентов удастся "выгрызть", не более. Ведь график уже вполне себе Чебышёвский: ТРИЖДЫ ошибка практически достигает своего максимума! По одноимённой теореме, оптимальное приближение ровно так и должно себя вести.
QuatCore

atan на ЧЕТЫРЁХ умножениях

Мишка такой человек — ему обязательно надо, чтоб от всего была польза. Когда у него бывают лишние деньги, он идёт в магазин и покупает какую-нибудь полезную книжку. Один раз он купил книгу, которая называется «Обратные тригонометрические функции и полиномы Чебышева». Конечно, он ни слова в этой книжке не понял и решил прочитать её потом, когда немножко поумнеет. С тех пор эта книга лежит у него на полке — ждёт, когда он поумнеет.
Николай Носов, "Весёлая семейка"


Пока не стал пытаться ещё сильнее оптимизировать коэффициенты для арктангенса на ДВУХ умножениях, поскольку подозреваю, что чудес не бывает - какие-то проценты ещё можно выжать, но мне бы хотелось повысить точность НА ПОРЯДОК. Так что, скрепя сердце, добавил ещё одно слагаемое и нашёл коэффициенты с помощью интерполяции по узлам Чебышёва.

Получил такой результат:



Напомню: мы имеем дело уже с отнормированным вектором, т.е x2+y2=1, и это офигительно упрощает дело! Плюс две формулы приведения, сначала к x≥0, затем к x≥|y|, т.е к углам от -45° до +45°.

В таких условиях данная формула даёт ошибку не более 4,9 угловой секунды! Этого более, чем достаточно для наших 16-битных вычислений, где цена младшего разряда углов составляет 12 угловых секунд.

График ошибки:


И в отличие от того нашего алгоритма с поворотами комплексных чисел, тут и вычислительных ошибок может набраться сильно меньше, за счёт 32-битного аккумулятора и "байпаса".

Collapse )

Наверное, так в итоге и сделаю. Ну, ещё сейчас надо глянуть, не слишком ли криво формулы приведения ложатся, надеюсь, что как-нибудь ляжет, может прям через i^k и регистр Inv, они же позволяют поменять два значения местами "виртуально", а также знаки поменять, не прибегая к убийственным ветвлениям... Забавно, как старая реализация, практически дошедшая до железа, устарела ещё до того, как хоть раз заработала. А всё потому, что МИША ПОУМНЕЛ.

Ещё забавно, что на вычислительной математике, на втором курсе, нам рассказывали, что интерполировать иногда лучше не через равные интервалы, а по узлам Чебышёва, но не стали особенно подчёркивать факт, НАСКОЛЬКО ЭТО КРУЧЕ! Даже тот пример, что был там приведён, оказывался каким-то беззубым - там что в лоб, что по лбу - примерно одно и то же получалось. Ну ладно, лучше поздно, чем никогда :)