nabbla (nabbla1) wrote,
nabbla
nabbla1

Categories:

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-й курс института.


По сути, у нас есть синус и косинус искомого угла, скажем, z, и надо найти угол. Причём sin(z) "загибается вниз" относительно прямой линии, а tg(z) - наоборот, устремляется вверх, то, может быть, нам взять их сумму с некоторыми весами, и будет нам счастье??

То есть, ввести приближение


Правда, чтобы получить тангенс, нам нужно деление, а без него хотелось бы обойтись. Зная, что косинус не шибко далеко убежит от единицы, можно взять приближение:



Вот и приходим к выражению


Причём здесь есть "произвол" сразу в обоих коэффициентах! Как их ни выбирай - выражение будет получаться нечётным, как и положено.

Для того, чтобы подобрать коэффициенты, преобразуем выражение следующим образом:


(раскрыли скобки, воспользовались синусом двойного угла)

Чтобы понять, "есть ли здесь какие-то перспективы", можно выбрать коэффициенты из разложения в ряд Тэйлора вблизи единицы:



Хотим, чтобы коэффициент при z в правой части был единичным, а при z3: нулевой. Выходит:




Откуда α=4/3, β=-1/3. Так получается приближение с "красивыми коэффициентами":


И как водится для ряда Тэйлора, это выражение шикарно работает вблизи нуля, но к краям ошибка растёт довольно резко. На 45° (когда x=y≈0,707) ошибка составляет 0,53°. И это не так уж плохо: примитивное приближение atan(y/x)≈y давало ошибку 4,49°, а более хитрое линейное приближение atan(y/x)≈1,0811y: 1,2°. В этом что-то есть :)

Метод наименьших квадратов (МНК)
В прошлый раз я настаивал на минимаксе, но сейчас решил применить старый добрый МНК, тупо потому что он проще :) По крайней мере привычнее. Минимакс от одного параметра - ещё куда не шло, а от двух параметров - чего-то не вполне понимаю, как подступиться.

А тут вводим интеграл от квадрата ошибки на всём интересующем нас диапазоне углов [0;A]:


Надо подобрать такие α, β, чтобы его минимизировать. Но проще всего сначала этот интеграл взять :)

Раскрываем скобки:


И проинтегрируем каждое слагаемое. Первое - мега-халява:

и в то же время абсолютно бесполезное слагаемое, т.к не зависит от α и β, можно было его сразу вычеркнуть.

Следующее заставляет вспомнить интегрирование по частям. Я его так и не запомнил, выписываю сначала производную произведения:


потом соображаю, что здесь нам подойдёт выражение x·cos(x):


Переношу в левую часть нужный нам кусочек -x·sin(x), а в правую всё остальное:


И только тогда вставляю знаки интеграла:


И наконец, добавляю недостающие 2α и границы интегрирования:


Да, мне очень стыдно :)

Следующее слагаемое очень похоже на это, но с двойкой в аргументе косинуса. Не буду торопиться, а то потом замучаюсь ошибку искать. Сначала производную выписываем:


Выражаем -x·sin(2x) и сразу всовываем неопределённый интеграл:


Добавляем множитель β и границы интегрирования:




В следующем слагаемом нужно использовать формулу произведения синусов, но её я тоже никогда не мог запомнить. Поэтому выписываю косинус от суммы, вот её почему-то помню:

(может потому, что если подсунуть α=β, отсюда придём к косинусу двойного угла, который "косинус квадрат минус синус квадрат". Был бы справа знак "+" - вышел бы полный бред, единичка)

Ну и косинус разности, который отсюда следует, т.к косинус чётная функция, а синус нечётная:

(а тут если α=β единица в правой части будет вполне логична :))

Если из верхней формулы вычесть нижнюю и поделить на два, практически устно получается то что нам надо:


И подставим наши x и 2x:


Ну а интегрируется оно в пол-пинка:




Далее идёт слагаемое с sin2x, и тут "рука тянется" сразу помножить на 1/2 и на длину отрезка (как это в рядах Фурье и при вычислении мощности в цепи переменного тока, и много где ещё), но здесь так нельзя, т.к у нас ни полного периода, ни половины периода не набирается! Надо всё честно:





И ещё не забыть помножить на α2.

И последнее слагаемое очень похожее:



И ещё помножить на β2/4.

Всё, собираем вместе:



Вот теперь пора брать частные производные по α и β и приравнивать их к нулю:





Это система линейных уравнений относительно α и β. Выписывать в явном виде выражения для α и β я не решился - завёл эти 6 коэффициентов в эксель (4 элемента матрицы и 2 элемента вектора), потом решение по правилу Крамера.

В итоге, подставив A=π/4, я получил α=1,379206316, β=-0,382389721267396. Они относительно близки к 4/3 и -1/3, полученные из ряда Тэйлора. Но другие, причём их сумма равна 0,997, вот не единица :)

С такими значениями максимальная ошибка получается на угле π/4, в 0,077° График ошибки выглядит вот так:


(по оси X угол в градусах, косинус и синус которого у нас в качестве входных данных. По оси Y ошибка данного метода, также в градусах)

Замечаем характерную особенность метода наименьших квадратов: ему важнее снизить ошибку на длинном участке графика (там ошибка не более 0,032°), а высокий, но узкий пик в конце диапазона не вносит большого вклада в нашу "интегральную" ошибку.

Тэйлор нервно курит в сторонке: уменьшили ошибку в 16 раз! Но если нас интересует минимакс (минимизация максимальной ошибки), то можно ещё поиграться...

Я пока не стал "всё начинать заново" и теперь подходить к этой проблеме со стороны минимакса. Вместо этого просто стал увеличивать значение A, то есть "как бы расширять диапазон". В итоге, до значения A=0,825 максимальная ошибка на диапазоне от 0 до π/4 УМЕНЬШАЛАСЬ, а потом снова начала расти. На A=0,825 мы получили ровно те значения, что приведены в начале поста, а график ошибки выглядит так:


Да, здесь всё как мы любим: ошибка "в плюс" сравнялась с ошибкой "в минус", но главное, уменьшилась по величине почти что вдвое!

Но я не уверен, что "выжали" из этой формулы всё, что смогли, ведь от значений, полученных по "методу наименьших квадратов", мы вели поиск "по одной оси". Возможно, если поставить целью минимакс конкретно на [0;π/4], получится ещё лучше, хотя сомневаюсь, что существенно.

Кстати, я ещё опробовал МНК на простейшем линейном приближении. Скажем, на [0;π/4] мы предлагали формулу atan(y/x)≈1,0811y и получали максимальную ошибку в 1,2 градуса. А МНК предлагает 1,0634y, что даёт максимальную ошибку в 1,92 градуса, но, понятное дело, меньшую среднеквадратичную.

И примерно такое же соотношение для других диапазонов.


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

  • Тестируем atan1 на QuatCore

    Пора уже перебираться на "железо" потихоньку. Решил начать с самого первого алгоритма, поскольку он уже был написан на ассемблере. В программу внёс…

  • Формулы приведения, что б их... (и atan на ТРЁХ умножениях)

    Формулу арктангенса на 4 умножениях ещё немножко оптимизировал с помощью алгоритма Ремеза: Ошибка уменьшилась с 4,9 до 4,65 угловой секунды, и…

  • Алгоритм Ремеза в экселе

    Вот и до него руки дошли, причина станет ясна в следующем посте. Изучать чужие библиотеки было лениво (в том же BOOSTе сам чёрт ногу сломит), писать…

  • Post a new comment

    Error

    Anonymous comments are disabled in this journal

    default userpic

    Your IP address will be recorded 

  • 1 comment