Как посчитать арккосинус на ПЛИС, не привлекая внимания санитаров

nabbla1 — 09.11.2025 С косинусом частично разобрались, хотя ещё будет запоздалое хэллоуинское (пока можно начать думать, при чём тут вообще хэллоуин?) продолжение. А теперь - обратная функция.

Подозревал, что с ней будет хуже: косинус всё-таки очень хорошая функция, с производной, никогда не выходящей за диапазон [-1;1]. А вот арккосинус совсем другой - у него в двух крайних точках производная устремляется к бесконечности. Поэтому аппроксимировать его полиномами - занятие печальное...

Честно говоря, у меня нет задачи во что бы то ни стало посчитать его "быстро". Текущая задача - посчитать косинус и арккосинус "малой кровью". Так-то, зная, что нужно вычислять "в обе стороны", хорошим выходом был бы CORDIC. Но с моим "перфекционизмом" я его буду неделю минимум "собирать", придумывать управляющую логику, весь datapath. И опять блоки ЦОС (DSP) окажутся "не при делах", что тоже немного странно - полезная же штука!

Хотелось решить задачу "радиотехнически"...


Но прежде чем сразу же отбросить затею с полиномами, я попробовал всё-таки их подобрать.

В линейном приближении максимальная ошибка составила 1724 единицы младшего разряда, или 9,5°.
Квадратичное: 995 единиц, или 5,5°.
Куб: 674 единицы, или 3,7°.
4я степень: 511 единиц, или 2,8°.
5я степень: 410 единиц, или 2,2°.
6я степень: 342 единицы, или 1,9°.
7я степень: 294 единицы, или 1,6°.
8я степень: 257 единиц, или 1,4°.

И при этом вообще нет уверенности, что при столь длинных вычислениях у меня не накопится нефиговых таких ошибок. Там, как обычно, феномен Рунге, когда полиномы высокой степени зверски "звенят", чуть ошибёшься с вычитанием близких значений - и привет.

И да, я-то хочу точности в 0,005°, в одну единицу младшего разряда, или даже лучше. И как видно, сходимость чрезвычайно медленная, совсем не похожая на ту, что мы наблюдали в косинусе (каждое новое слагаемое улучшает точность на порядок).




На графике приведена величина, обратная к ошибке, в зависимости от степени полинома. Т.е Степень полинома * ошибка - величина практически постоянная! Чтобы получить ошибку в 26 единиц, понадобится полином 80-й степени! А до 3 единиц - и вовсе 800й. Это абсурд, и так мы делать не будем!

Один перспективный вариант - интерполяция кубическими сплайнами, т.е вместо полинома, приближающего функцию ЦЕЛИКОМ, мы делим кривую на несколько отрезков, на каждом отрезке своя кубическая парабола, причём они гладко "склеены" между собой. Возможно, ещё его сделаем. Но тоже уверенности нет. Я помню, что "колокола" Гаусса сплайнами интерполировались очень хреново, очень им хотелось начать "звенеть" вблизи нуля.

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



Операционный усилитель потому так и назвали, что с его помощью можно выполнять всевозможные математические операции, чем активно пользовались в аналоговых компьютерах и не только. Тут он попытается уравнять свои входы, начнёт тащить выходное напряжение в нужную сторону, и в конце концов у него это удастся, так что на выходе мы получим квадратный корень.

Хотя можно было аналоговые схемы и не вспоминать, вычислительная математика у нас с той же задачи началась, найти корень уравнения, то есть такой x, при котором f(x) = 0, притом полагается, что считать f(x) мы уже умеем.

Но, зараза, производная, варьирующаяся от 1 до бесконечности, делает большинство методов, которые мы учили, совершенно бесполезными!

Один из самых простых как раз напоминал чем-то работу операционника. Вот мы взяли "нулевое приближение", что y₀ = π/2-x ≈ acos(x). Затем нашли cos(y₀). Если мы "угадали", это будет в точности x. Вычитая x, мы находим "невязку", её мы с каким-то заранее подобранным коэффициентом прибавляем к y₀, и это приближение будет уже получше. Снова от него находим cos(), вносим следующую коррекцию, и так, пока не придём куда надо.

Что ж, посчитаем "на пальцах" acos(1/2). Нулевое приближение: π/2-1/2 ≈ 1,0708 (или 61°). Берём от него косинус, получаем 0,479. Вычитаем x, то есть 1/2, получаем -0,0205. Прибавляем к нашему приближению, выходит 1,050222 (или 60,17°), косинус: 0,497, вычитаем x, получаем -0,00262. Прибавляем к нашему приближению, выходит 1,0476, или 60,023°.

Выглядит очень неплохо! Но это лишь потому что мы ещё не подошли к краю. А попробуем найти acos(1). Поначалу даже ощущается какой-то прогресс. Но допустим, мы мучительно добрались до приближения 0,01. Считаем от него cos(), получаем 0,99995. Ну да: когда производная обратной функции стремится к бесконечности, у "прямой" функции она фактически нулевая. Поэтому наша оценка близости через косинус - очень кривая. Нам говорят: вы уже почти у цели, сдвиньтесь пожалуйста на 50 пикометров влево. Хм, и ещё на 50. И ещё на 50, вы на верном пути!

Но, если подумать, есть самый дурацкий способ, который здесь сработает очень неплохо! Это старое доброе "деление пополам". Задача как будто специально создана для него. Мы хотим решить уравнение cos(x) = A, при этом x находится в интервале от 0 до π, A - заведомо от -1 до +1. Причём мы заранее можем сказать, что cos(0)=1, как минимум не ниже A, тогда как cos(π)=-1, как минимум, не выше.

Для "деления отрезка пополам" не нужно вообще ничего помнить, кроме текущего разряда двоичного числа.

Т.е на нулевой итерации у нас полный интервал 0000_0000_0000_0000..0111_1111_1111_1111. (мы продолжаем использовать масштаб, когда на 16-битные числа, хоть со знаком, хоть без знака, ложится один полный оборот, например, углы от -π до π радиан)
Чтобы поделить его на два, мы берём косинус от 0100_0000_0000_0000 = 0x4000, что соответствует углу π/2.
Если результат получился выше, чем A, по текущему разряду устанавливаем единичку, иначе - нолик. Допустим, мы снова хотим найти acos(1/2). Первым вы взяли cos(π/2) = 0, это ниже, чем 1/2, поэтому устанавливаем нолик.

Условно, мы сейчас в интервале 0000_0000_0000_0000..0011_1111_1111_1111, хотя всё это запоминать вообще не требуется! Нам ничего толком не нужно, кроме регистра на 16 бит, куда мы будем один за другим вставлять правильные биты. Т.е на первом шаге у нас было 0000_0000_0000_0000, жирным выделен УЖЕ ИЗВЕСТНЫЙ бит. Стало 0000_0000_0000_0000.

Поехали дальше, ставим единицу в следующий неизвестный бит: 0010_0000_0000_0000 = 0x2000, это угол π/4. Косинус получатеся 0,707, это выше, чем 1/2, поэтому устанавливаем единицу: 0010_0000_0000_0000.

Ставим единицу в следующий неизвестный бит: 0011_0000_0000_0000 = 0x3000, это 67,5°. Косинус: 0,383, ниже, чем 1/2, поэтому устанавливаем ноль: 0010_0000_0000_0000.

Ставим единицу: 0010_1000_0000_0000 = 0x2800, это 56,25°, косинус: 0,555, выше, чем 1/2, поэтому оставляем единицу: 0010_1000_0000_0000. И так далее!

Да, с нахождением acos(1/2) ранее показанный итеративный метод справлялся лучше, потому что там мы довольно хорошо "угадали" производную, которая у нас была "зашита" в качестве константы, на которую домножаем невязку. Зато метод "деления на 2" ЖЕЛЕЗНО СРАБОТАЕТ ЗА 16 итераций ВООБЩЕ ВСЕГДА.

Пора реализовать.

Довольно симпатичный модуль, всего 28 строк кода, да и те довольно "рыхлые":



Примерно так он выглядит "схемотехнически":



Регистры довольно просты:
- rD - наиболее тривиальный, он просто сохраняет входное значение, потому что я пока не уверен, что все 16 тактов оно будет здесь держаться (скорее всего в нас "кинут" значением и займутся другими делами, чтобы потом проведать, как мы там поживаем),
- curBit - сдвиговый регистр, у которого всего два режима. Либо он сдвигается вправо, заполняясь слева нулями, либо (по сигналу req) в старший бит заносится единица, во все остальные - нули,
- SAR - "регистр последовательного приближения". У него приоритетный управляющий вход sset, также подключённый к req: по нему регистр устанавливается в "нулевое приближение". Для арккосинуса, так уж получилось, что это НОЛЬ. Если же sset=0, то если ENA=1, он по фронту защёлкнет значение, пришедшее на его вход D, иначе продолжит хранить старое значение.

Ещё у нас стоит сумматор, за ним - ранее сделанный "косинус", и ещё компаратор. Тут я снова пошлю лучи ненависти Альтере: их модуль lpm_compare, СПЕЦИАЛЬНО ДЛЯ ЭТОГО СДЕЛАННЫЙ, раз за разом синтезируется очень толсто. Когда числа беззнаковые, я всегда вместо него ставлю lpm_add_sub, т.е "вычитатель", причём волнует меня исключительно бит переноса cout. Сейчас хотел его же поставить, но числа со знаком, в дополнительном коде, там я продолжаю путаться (отчасти и от реализации lpm_add_sub: кажется, если ему указать, что числа со знаком, то cout вообще работать перестаёт, только overflow) и желать поскорее перевести весь мир на ур. тр. систему, где подобной фигни в принципе нет.

А чтобы понять, как оно работает, посмотрим результаты симуляции:



Мы "затребовали" (request) вычислить acos(1/2). Регистр последовательного приближения SAR (он же и наш выход) сбрасывается в ноль (хотя он и так был в нуле). Сдвиговый регистр устанавливается в 0x4000. Их сумма, 0x4000, подаётся одновременно на вход D регистра SAR, и на вход cos. Выход cos и входное значение D, "на всякий случай" сохранённое внутри, сравнивается между собой. Если Dsign = 1, значит cos(Q) вышел больше чем D.

На первом такте мы посчитали cos(0x4000), т.е от π/2. Правильный ответ 0, но у нас косинус чуть кривоватый пока, он дал -5, в смысле -5/32768≈-0,00015. В любом случае, это меньше, поэтому Dsign = 0, новое значение в SAR не защёлкивается.

На втором такте на сдвиговом регистре 0x2000, оно складывается SAR, что даёт 0x2000, или π/4. Косинус от этого значения: 23171/32768≈0,707122 (правильное значение: 0,707106, или 23170, т.е ошибка на единичку). Теперь "зажёгся" Dsign=1, и новое значение защёлкнулось в SAR.

Так продолжается некоторое время. Единичка в правой части сдвигового регистра "вдвигается" в регистр ready (на "схеме" я изобразил его частью сдвигового регистра, так оно и есть по сути), а после этого остаются сплошные нули, что естественным образом прекращает работу без необходимости в явных счётчиках, конечных автоматах и пр. Понятно, всё есть конечный автомат, и это тоже, нодо сих пор мне больше по душе декомпозиция "в пространстве" (на отдельные составные части, достаточно слабо связанные друг с другом), чем "во времени" (в состоянии 1 делай то-то и тото, перейди в состояние 2, а в состоянии 2 поступай вот так!).

Результат работы: 10922/32768 * π радиан, или 10922/32768 * 180 градусов, то есть 59,996°. Следующее значение, 10923 - это 60,0018°, так что НЕДУРНО, я считаю!

Посмотрим крайние значения. Сначала попытаемся найти acos(32767/32768):



Довольно забавный результат: НОЛЬ. Это, вообще, хорошо. Да, правильный ответ вроде как 0,45°, или 81 единица, но вся проблема, что у нас нет единицы, это самое близкое к ней. И получилось, что поскольку и cos(0) выдал свой "максимум" в 32767, то при нахождении обратной величины мы "как и положено" получили ноль.

Может, как-нибудь расскажу, зачем всё это мне понадобилось, по сути сейчас нужно выступить "адвокатом дьявола" и опробовать альтернативный кватернионам алгоритм работы. И знаете, можно не любить кватернионы, но когда вместо них попытаешься использовать углы Эйлера, непрерывно применять сферические теоремы косинусов для пересчётов, вдруг понимаешь, что кватернионы - это меньшее из зол!

Теперь посмотрим acos(32766/32768):



Ответ получился 90/32768 * 180°, или 0,49°. Правильный ответ: 0,63°. Это всё ещё следствие странно "сплющенного" косинуса вблизи нуля.

И ещё надо отрицательные числа проверить, acos(-1):


Ответ: 32767/32768 * 180°, или 179,9945°. Правильный ответ: 180°, но теперь на выходе мы упёрлись в потолок, это самое большое представимое число.

И напоследок acos(-32767/32768):



Ответ 32667/32768 * 180°, или 179,445°. Правильный ответ: 179,552°

Что ж, для моих нужд этого пока достаточно. Когда я чуть улучшу работу cos(), "автоматом" улучшится и acos().

Оставить комментарий

Популярные посты:

• Ранее
• Архив