nabbla (nabbla1) wrote,
nabbla
nabbla1

Categories:

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

Формулу арктангенса на 4 умножениях ещё немножко оптимизировал с помощью алгоритма Ремеза:



Ошибка уменьшилась с 4,9 до 4,65 угловой секунды, и дальше уменьшать особо некуда.

Если эти коэффициенты записать в 16 битах, выйдет так:



точность немного упадёт, до 6 угловых секунд:


Учитывая, что это половина цены младшего разряда - вроде всё хорошо.

Но только начал писать программу на ассемблере QuatCore, так сразу же упёрся в то, как привести угол в интервал от -45° до 45°, в котором данная формула даёт хорошие результаты...


Первым делом я собирался посмотреть знак числа x и поместить его в 1-битный регистр Inv:
Acc	[X+k]
SUB	0
Inv	S


Этот знак нам нужен, чтобы дополнительно "перевернуть" знак y, иначе мы II и III квадранты перепутаем местами.

Далее надо ещё проверить, что |x|>|y|. С этим тоже проблем вроде нет:
ABS	[X+k]
ABSS	[X+1]


Сейчас во флаге S лежит знак разности |x|-|y|. Его тоже можно куда-нибудь положить, например, в регистр k.

И тут-то мы и обнаруживаем: если |x|<|y|, то нужно отдельно рассмотреть ещё два случая: y>0 (углы 45°..90°) и y<0 (углы -45°..90°)!

Выходит, нужно аж 6 случаев рассматривать:

1. x≥0, |x|≥|y|: углы -45°..+45°. Самый простой случай, формула применяется напрямую:


2. x≥0, |x|<|y|, y≥0: углы +45°..+90°. В формуле надо поменять местами x,y, результат взять со знаком "-" и ещё прибавить π/2:



Другая интерпретация: вектор (x;y) заменить на (y;-x), т.е осуществить поворот на 90° по часовой стрелке, из-за чего результат будет на эти же самые 90° меньше, поэтому "ручками" их добавляем назад.

3. x≥0, |x|<|y|, y<0: углы -45°..-90°. В формуле надо поменять местами x,y, не забыть "убрать знак" y (т.к он сейчас стоит на месте x, который в этой формуле полагается положительным) и ещё ВЫЧЕСТЬ π/2:



Другая интерпретация: вектор (x;y) заменяется на (-y;x), т.е идёт поворот на 90° против часовой стрелки, а потом из результата вычитаем эти 90° назад.

И ещё 3 случая, где x<0.

Совершенно очевидно, что мы должны проверить знак x, и что мы должны посмотреть, что больше по модулю, x или y. Но ещё и знак y тоже оказывается важен, на его основе выбирается знак π/2. Т.е нам нужно 3 сравнения и все эти 3 бита где-то хранить и как-то использовать для выбора знаков и множителей...

А ведь формула даёт результаты в диапазоне 90° (-45..+45), поэтому можно было ожидать необходимости уполовинивания лишь ДВАЖДЫ, но никак не ТРИЖДЫ!

Вся проблема, что диапазон очень неудобный. Куда проще было бы, если мы могли бы находить atan(y/x) ТУПО В I квадранте (x≥0, y≥0), это те же 90°, и потом лишь останется подставить правильный знак! Мы бы могли посчитать x·y, это и даст знак результата. А затем x,y брать просто по модулю.

Тут, правда, мне очередная шлея попала под хвост: а как бы так полезно использовать это произведение x·y, не только для знака?

Попытался найти приближение к арктангенсу в I квадранте в виде:


и обнаружил, что этот x·y здесь на этот сам x·y не сдался! Для диапазона -45°..+45° логично было требовать нечётных выражений, а здесь эквивалентным требованием является выполнение arctan(x,y) = π/2-arctan(y,x). И действительно, "линейное приближение" таково:



Максимальная ошибка достигается на краях диапазона, т.е при x=0 либо y=0, и составляет 1,2°. Да, мы это проходили, только тогда диапазон был -45°..+45°, а сейчас 0°..90°.

По сути, единственная чётная компонента формулы - это константа π/4, а все остальные слагаемые должны при замене (x,y) на (y,x) менять знак. Но x·y этого явно делать не будет!

Подумал грешным делом: может ну её, симметрию, как-нибудь Чебышёв на пару с Ремезом найдут, как этот множитель применить с пользой? Попробовал интерполировать по узлам Чебышёва - тут же получил вырожденную систему уравнений. Видимо, не судьба :)

Но может быть, всё-таки разбить наш квадрант на 2 диапазона и, глядишь, в них уже этому x·y найдётся применение, за счёт меньшего диапазона можно будет снизить количество слагаемых, а за счёт симметрии мы получим одни и те же коэффициенты в двух диапазонах, разве что с разным порядком и знаками.

Я стал искать формулу для диапазона 0°..45° в следующем виде:



Провёл интерполяцию по узлам Чебышёва - и получил следующий график ошибки:


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

Именно здесь у меня и возникло жгучее желание опробовать алгоритм Ремеза!

И он отработал как надо, выдав следующую формулу:


И сразу же покажу формулу для диапазона 45°..90°:



Как я и предполагал, коэффициенты по сути остаются те же самые, только меняются местами и знаками, ну и π/2 появляется.

В принципе, эти формулы я наверное мог бы достаточно компактно исполнить на QuatCore. График ошибки на всём интервале выходит таким:


Да, Ремез отработал превосходно: все максимумы и минимумы ошибки одинаковы по абсолютной величине, и ошибка эта составляет 0,0146°, или 52 угловые секунды.

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


В общем, Буриданов осёл нервно курит в сторонке. Столько разных вариантов, не знаю, что выбрать.

- Есть "старый добрый" алгоритм с поворотами комплексного числа. Тупой до безобразия, но вполне рабочий, и фактически уже реализованный на ассемблере. Он ещё и длину вектора умеет находить! (Т.е либо ему скармливаешь вектор единичной длины - и он находит угол и ещё раз подтверждает единичную длину, либо даёшь вектор произвольный - он его поворачивает практически на нулевой угол, после чего довольно точно даёт длину, но угол "плюс-минус лапоть". Но как ни странно, меня это устраивает) С точностью может быть небольшая засада: за 5 умножений комплексных чисел с хранением результатов в памяти (в 16 битах точности) там накопится некоторая ошибка в младших разрядах.

- Можно atan на 4 умножениях "повернуть" на 45° (заменив y на (|y|-|x|) и x на (|y|+|x|)/2), а знак результата положить равным знаку x·y. Это даст уже 5 умножений, но вроде может получиться весьма компактно и НАИБОЛЕЕ ТОЧНО, поскольку все промежуточные результаты не будут покидать АЛУ, где они хранятся в 32 битах.

- А можно вот этот "компромисс" на 3 умножениях: точность хуже на порядок, но в теории это будет самый быстрый и компактный алгоритм.

Боюсь, в конечном итоге всё равно перепробую все варианты, и только тогда соображу, что мне надо, обычно так бывает. Если вдруг не "припрёт". Сейчас просто затишье: в РКК народ в отпуске до 1 октября, у нас тут тоже никого нет толком. Как вернутся - уже такой фигнёй не займёшься, вернусь к родным бумажкам. Надо ж написать программу ЛОИ (Лабораторно-Отработочных Испытаний), а это не такая программа, что для компьютера, а бумажка такая, нудная до безобразия. Потом ещё будет отчёт по ЛОИ. И, если хватит времени между этими двумя бумажками - может и сам ЛОИ удастся провести...

Poll #2113299 Как мне посчитать арктангенс

Столько вариантов "накопал", а выбрать не могу

На последовательных поворотах комплексного числа
0(0.0%)
На 4 умножениях
2(28.6%)
На 3 умножениях
0(0.0%)
Опробовать все и выбрать наиболее понравившийся
5(71.4%)
Tags: ПЛИС, математика, программки, работа, странные девайсы
Subscribe

Recent Posts from This Journal

  • Ремонт лыжных мостиков

    Вернулся с сегодняшнего субботника. Очень продуктивно: отремонтировали все ТРИ мостика! Правда, для этого надо было разделиться, благо народу…

  • Гетто-байк

    В субботу во время Великой Октябрьской резни бензопилой умудрился петуха сломать в велосипеде. По счастью, уже на следующий день удалось купить…

  • А всё-таки есть польза от ковариаций

    Вчера опробовал "сценарий", когда варьируем дальность от 1 метра до 11 метров. Получилось, что грамотное усреднение - это взять с огромными весами…

  • Post a new comment

    Error

    Anonymous comments are disabled in this journal

    default userpic

    Your IP address will be recorded 

  • 0 comments