nabbla (nabbla1) wrote,
nabbla
nabbla1

Categories:

Таблица коэффициентов для atan1

Пора уже добить этот несчастный арктангенс, что-то слишком долго его ковыряю. Мы упоминали, что в оперативной памяти будет лежать таблица в 15 слов: 5 углов и 5 комплексных чисел, им соответствующих. Вопрос был, в каком порядке всё это дело расположить, ну и сами константы "придумать". Пожалуй, вот так:

;таблица поворотов для процедуры atan1
AtanTable:
;завершающая итерация
;пока длина вектора составляет 1,0000032 от исходной (в теории, т.е если коэффициенты заданы ровно эти, но вычисления с абсолютной точностью)
Atan4a	Int16 1608		;pi/32 в формате Q2.14
;если cos(pi/32) и sin(pi/32), это было бы 32610 и 3212.
;если cos(1608/16384) и sin(1608/16384), то 32610 и 3211
;но мы ещё должны поделить на 2, т.к формат 2.14, а вектор был в 1.15
;(и это деление проведём, заменив ADD на DIV2A в конце процедуры)
;а ещё умножить на k=1,00120570556/1,0000032
;в итоге возьмём cos(1608/16384)*k и sin(1608/16384)*k
Atan4x	Int16	32650
Atan4y	Int16	3215


;к этому моменту длина вектора составляет 1,000013 от исходной
Atan3a	Int16 3217		;pi/16 в формате Q2.14
Atan3x	Int16 32138		;exp(i*pi/16) = 0,9808 + 0,1951i
Atan3y	Int16 6393

;к этому моменту длина вектора составляет 1,00000107 от исходной (в теории, это sqrt(23170^2+23171^2)/32768)
Atan2a	Int16 6434		;pi/8 в формате Q2.14
Atan2x	Int16 30274		;exp(i*pi/8) = 0,9239 + 0,3827i
Atan2y	Int16 12540

;к этому моменту длина вектора не менялась
Atan1a	Int16 12868		;pi/4, выраженные в формате Q2.14 (в который влезают значения от -4 до 4)
Atan1x	Int16 23170		;exp(i*pi/4) = 0,7071 + 0,7071i
Atan1y	Int16 23171		;а точнее, exp(i * 12868 / 16384), чтобы осуществить именно тот поворот, что указан в Atan1a

Atan0a	Int16 0
Atan0x	Int16 0		;первый поворот, на +- 90 градусов (0+1i)
Atan0y	Int16 -32768	;взяли знак "-" потому что первый раз мы берём обратный знак (из нуля вычитаем "x")
;длина вектора остаётся прежней


Как всегда, здесь всё задом наперёд!

Кажется, во все времена на ассемблере было удобнее организовывать циклы ДО НУЛЯ, а не от нуля куда-то вверх. Поэтому таблицу эту мы проходим снизу вверх.

На первой итерации мы вместо exp(i*pi/2) поставили exp(-i*pi/2), поскольку с точки зрения реализации оказалось перед входом в цикл загрузить в аккумулятор ноль и вычесть из него значение x, отчего знак выходит обратным. Ничего страшного.

Очень примечательное значение Atan0a = 0. Именно оно отличает atan1 от atan2 :) Поставь сюда π/2 - и метод будет отличать друг от друга все 4 квадранта, а так он "путает" первый с третьим и второй с четвёртым. Напомню: ТАК И ЗАДУМАНО, чтобы при переходе от половинного угла в кватернионе к целому углу не получить абсурдный диапазон -2π..2π и потом героически его ужимать. Зачем...

Далее, присутствует непонятная на первый взгляд асимметрия в значении exp(i*pi/4). Это поворот на 45 градусов, должны получиться два одинаковых значения, но не тут-то было! Дело в том, записать точно угол pi/4 мы не смогли (значение Atan1a), ближайшее, что было: 12868/16384. И по логике вещей, раз уж мы "заявили" этот угол (прибавили его к результату), то и поворот надо осуществить именно на этот угол! Это оказывается угол 45,00013 градусов. Косинус от него, умноженный на 32768, равен 23170,423, а синус: 23170,527. Вот мы и округляем одного вниз, а второго: вверх!

Значения для 3-й и 4-й итерации вполне нормальные, "без особенностей". Что мы брали бы exp(i*pi/8) и exp(i*pi/16), что exp(i*12540/16384) и exp(i*6393/16384), результаты одни и те же.

И наконец, на 5-й итерации мы решили вместе с поворотом провести и небольшое масштабирование, которое существенно повышает точность. В прошлый раз я пришёл к коэффициенту 1,0012 просто "методом научного тыка", а тут попытался найти "точное выражение".

Мы пытаемся заменить asin(x) на простейшее линейное выражение kx. Так оно выглядит в интересующем нас диапазоне, от 0 до 5,625 градусов (от -5,625 до 0 всё то же самое, только знак обратный, обе функции нечётные):


Да, здесь вообще ничего не понятно, графики сливаются. Для наглядности возьму куда более широкий диапазон, до 45 градусов:


Арксинус "загибается вверх", а хорошее линейное приближение должно так его пересечь, чтобы максимальная ошибка на первом участке (где арксинус проходит ниже, чем надо) равнялась максимальной ошибке в самом конце интересующего нас диапазона (где арксинус рванул вверх), два отрезка, нарисованные чёрным цветом.

Это минимаксная задача - минимизировать максимальную ошибку. Может захотеться минимизировать среднеквадратичную ошибку, а ещё можно потребовать несмещённого результата, в некоторых случаях это страшно полезно, как в случае округления. Но этот арктангенс у меня будет в самом конце алгоритма, просто, чтобы показать результаты в понятном виде. Тут мне кажется, что минимакс наиболее обоснован.

Выразим ошибку приближения:


На первом отрезке это выражение уйдёт в минус, затем, к границе интересующего диапазона - в плюс.

Сначала найдём минимум, для чего берём производную и приравниваем её к нулю:



Мы вполне можем объявить неизвестной величиной именно xm (положение минимума), поскольку коэффициент k, который выражается через неё.

Абсолютная величина ошибки в этом минимуме составит:


И оно должно равняться ошибке на границе нашего диапазона, когда x=sin(π/32):



Приравняем их:



Сгруппируем:


Сделаем замену: x=sin(α). Тогда корень в знаменателе заменяется на cos(α), а арксинус - просто на α:



Проще, наверное, уже не получится...

Аналитического решения не имеет. Подумал было α выразить через всё остальное:


Надеялся, это окажется сжимающее отображение - и тогда я эту формулу восприму как итерации по уточнению α - и быстренько приду к ответу. ХРЕНУШКИ: сходиться оно никуда не хочет, плавно куда-то ползёт... Может, меня кто просветит, как "в пол-пинка" его модифицировать, чтобы итеративно сходилось. А пока я задолбался - и зашёл на WolframAlpha:


Притом с мобильного: на последней версии Оперы под Windows XP сайту совсем поплохело...

Да, чуда не случилось, Вольфрам тоже ничего лучше не нашёл, как численно решить. Нас интересует значение α=0,04908147042...

А коэффициент, который мы всё это время искали - это cosec(α)=1/cos(α)=1,00120570556...

Потом я ещё немножко заморочался - решил посмотреть, что происходило с длиной вектора во время каждого из поворотов. Насчитал, что она помножилась на 1,0000032. Не знаю, есть ли в этом какой-то смысл: если помножить это число на 32768, получим примерно 32768,105, т.е такое увеличение длины непредставимо. Как-нибудь поиграюсь, посмотрю - подозреваю, что если прямо все возможные векторы взять, увеличение длины всё-таки получится увидеть, т.е в нескольких местах, где "пропадающая" дробная часть колебалась вблизи 0,5, мы ошибочно склонились не в ту сторону. Не знаю пока...

Но всё же я этот коэффициент на всякий пожарный поделил на 1,0000032, потом помножил на cos(π/32) и sin(π/32) - и таким образом получил последние две константы для таблицы. Самое смешное, что эта добавочка возымела-таки эффект: без неё мы получали 32649,5314388705, что надо было округлить до 32650. А с ней: 32649,4269607042, с округлением до 32649 :)


Наверное, в очередной раз заморочался куда сильнее, чем следовало бы. Очень странная форма перфекционизма - придумать приближённое решение, но потом попытаться выжать из него все-все соки.
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 

  • 4 comments