atan на ЧЕТЫРЁХ умножениях
Мишка такой человек - ему обязательно надо, чтоб от всего была польза. Когда у него бывают лишние деньги, он идёт в магазин и покупает какую-нибудь полезную книжку. Один раз он купил книгу, которая называется «Обратные тригонометрические функции и полиномы Чебышева». Конечно, он ни слова в этой книжке не понял и решил прочитать её потом, когда немножко поумнеет. С тех пор эта книга лежит у него на полке - ждёт, когда он поумнеет.
Николай Носов, "Весёлая семейка"
Пока не стал пытаться ещё сильнее оптимизировать коэффициенты для арктангенса на ДВУХ умножениях, поскольку подозреваю, что чудес не бывает - какие-то проценты ещё можно выжать, но мне бы хотелось повысить точность НА ПОРЯДОК. Так что, скрепя сердце, добавил ещё одно слагаемое и нашёл коэффициенты с помощью интерполяции по узлам Чебышёва.
Получил такой результат:
\approx{y}(1{,}665436193-0{,}665222351x-0{,}168724504y^2))
Напомню: мы имеем дело уже с отнормированным вектором, т.е x2+y2=1, и это офигительно упрощает дело! Плюс две формулы приведения, сначала к x≥0, затем к x≥|y|, т.е к углам от -45° до +45°.
В таких условиях данная формула даёт ошибку не более 4,9 угловой секунды! Этого более, чем достаточно для наших 16-битных вычислений, где цена младшего разряда углов составляет 12 угловых секунд.
График ошибки:
И в отличие от того нашего алгоритма с поворотами комплексных чисел, тут и вычислительных ошибок может набраться сильно меньше, за счёт 32-битного аккумулятора и "байпаса".
Код, непосредственно реализующий эту формулу, выглядит примерно так:
SQRD2 [X+1]
C AtanGamma
MUL Acc
C AtanBeta
FMA [X+k]
ADD AtanAlpha
C [X+1]
MUL UAC
В первой строчке y возводится в квадрат, и 32-битный результат лежит в аккумуляторе.
В строке "MUL Acc" срабатывает "байпас", дополнительная шина, соединяющая младшие биты входа и выхода АЛУ. Благодаря нему, в умножении задействованы практически все 32 бита результата, что повышает точность.
Строка "FMA [X+k]" также хороша: это Fused Multiply-Add. Не зря я назвал её так, а не MAC (Multiply-Accumulate). Разница как раз в том, что FMA способна использовать аккумулятор большой ширины и продолжать накапливать в нём выражение с полной точностью.
И в самом конце, "MUL UAC", байпас срабатывает ещё разок. Вместо Acc тут используется UAC (Unsigned Acc), поскольку в скобке значение будет больше единицы, и ни в коем случае нельзя здесь вызывать "насыщение".
Наверное, так в итоге и сделаю. Ну, ещё сейчас надо глянуть, не слишком ли криво формулы приведения ложатся, надеюсь, что как-нибудь ляжет, может прям через i^k и регистр Inv, они же позволяют поменять два значения местами "виртуально", а также знаки поменять, не прибегая к убийственным ветвлениям... Забавно, как старая реализация, практически дошедшая до железа, устарела ещё до того, как хоть раз заработала. А всё потому, что МИША ПОУМНЕЛ.
Ещё забавно, что на вычислительной математике, на втором курсе, нам рассказывали, что интерполировать иногда лучше не через равные интервалы, а по узлам Чебышёва, но не стали особенно подчёркивать факт, НАСКОЛЬКО ЭТО КРУЧЕ! Даже тот пример, что был там приведён, оказывался каким-то беззубым - там что в лоб, что по лбу - примерно одно и то же получалось. Ну ладно, лучше поздно, чем никогда :)
Николай Носов, "Весёлая семейка"
Пока не стал пытаться ещё сильнее оптимизировать коэффициенты для арктангенса на ДВУХ умножениях, поскольку подозреваю, что чудес не бывает - какие-то проценты ещё можно выжать, но мне бы хотелось повысить точность НА ПОРЯДОК. Так что, скрепя сердце, добавил ещё одно слагаемое и нашёл коэффициенты с помощью интерполяции по узлам Чебышёва.
Получил такой результат:
Напомню: мы имеем дело уже с отнормированным вектором, т.е x2+y2=1, и это офигительно упрощает дело! Плюс две формулы приведения, сначала к x≥0, затем к x≥|y|, т.е к углам от -45° до +45°.
В таких условиях данная формула даёт ошибку не более 4,9 угловой секунды! Этого более, чем достаточно для наших 16-битных вычислений, где цена младшего разряда углов составляет 12 угловых секунд.
График ошибки:
И в отличие от того нашего алгоритма с поворотами комплексных чисел, тут и вычислительных ошибок может набраться сильно меньше, за счёт 32-битного аккумулятора и "байпаса".
Код, непосредственно реализующий эту формулу, выглядит примерно так:
SQRD2 [X+1]
C AtanGamma
MUL Acc
C AtanBeta
FMA [X+k]
ADD AtanAlpha
C [X+1]
MUL UAC
В первой строчке y возводится в квадрат, и 32-битный результат лежит в аккумуляторе.
В строке "MUL Acc" срабатывает "байпас", дополнительная шина, соединяющая младшие биты входа и выхода АЛУ. Благодаря нему, в умножении задействованы практически все 32 бита результата, что повышает точность.
Строка "FMA [X+k]" также хороша: это Fused Multiply-Add. Не зря я назвал её так, а не MAC (Multiply-Accumulate). Разница как раз в том, что FMA способна использовать аккумулятор большой ширины и продолжать накапливать в нём выражение с полной точностью.
И в самом конце, "MUL UAC", байпас срабатывает ещё разок. Вместо Acc тут используется UAC (Unsigned Acc), поскольку в скобке значение будет больше единицы, и ни в коем случае нельзя здесь вызывать "насыщение".
Наверное, так в итоге и сделаю. Ну, ещё сейчас надо глянуть, не слишком ли криво формулы приведения ложатся, надеюсь, что как-нибудь ляжет, может прям через i^k и регистр Inv, они же позволяют поменять два значения местами "виртуально", а также знаки поменять, не прибегая к убийственным ветвлениям... Забавно, как старая реализация, практически дошедшая до железа, устарела ещё до того, как хоть раз заработала. А всё потому, что МИША ПОУМНЕЛ.
Ещё забавно, что на вычислительной математике, на втором курсе, нам рассказывали, что интерполировать иногда лучше не через равные интервалы, а по узлам Чебышёва, но не стали особенно подчёркивать факт, НАСКОЛЬКО ЭТО КРУЧЕ! Даже тот пример, что был там приведён, оказывался каким-то беззубым - там что в лоб, что по лбу - примерно одно и то же получалось. Ну ладно, лучше поздно, чем никогда :)