nabbla (nabbla1) wrote,
nabbla
nabbla1

Category:

atan1 "на четырёх умножениях" на ассемблере

Немножко подправил таблицу результатов для первого алгоритма. В этот раз вместо "правильных углов" 0°, 12°, 24° и т.д (те, которые я ХОТЕЛ взять для тестов) я взял те целочисленные значения X/Y, которые в действительности выступили входными данными для процедуры. Думал, что ошибка немного уменьшится, т.к часть ошибки была вызвана неточным представлением этих углов. На деле она возросла, как ни странно:



Только в 4 случаях из 15 получили все 16 бит верно, ещё в 8 случаях ошиблись на единичку младшего разряда, и в трёх случаях - на две единицы (в прошлый раз, как ни странно, такой случай был единственный).

В целом, поведение ожидаемое, когда производится 5 умножений комплексных чисел с хранением промежуточных чисел в 16-битном виде. Если каждое округление до 16 бит даёт нам случайную ошибку в 0,316 единицы младшего разряда (таково среднеквадратичное отклонение при равномерном распределении чисел перед округлением, один на корень из 10), и ошибки независимы друг от друга, то поскольку у нас 5 раз сохраняется по 2 числа (т.е действительная и мнимая часть), то ожидаем увеличения ошибки ещё в корень из 10 раз, до 1 единицы младшего разряда. Но распределение будет примерно Гауссовым, с дисперсией 1, и там вероятность, что ошибка будет [-0,5;0,5], т.е ни один бит не пострадает: 38%. Вероятность ошибки в 1 единицу младшего разряда (т.е [-1,5;-0,5] и [0,5;1,5]): 48%, в две единицы и выше: 14%. Худо-бедно похоже на то, что мы увидели.

Вообще, такая точность меня вполне устраивает, но "из спортивного интереса" хочу попробовать реализовать более быстрый и точный алгоритм.


Хочу не трогать слишком много регистров "по чём зря", тем более, что сейчас результат ожидается всего один (угол), пущай он "возвращается" прям в аккумуляторе (не нужно его класть куда-то в память). И ни на какую таблицу ссылаться мы не собираемся, Ну разве что индексные регистры (и Inv впридачу) надо сохранить, так что "заготовка" процедуры такова:
;арктангенс через Чебышева
;x лежит в [X], y лежит в [X+1], результат выдаётся в аккумуляторе
;затирается регистр C и флаг знака, остальные остаются "как есть"
atan1 proc
	[SP++]	ijk
	
	
	
	
	ijk		[--SP]
	JMP		[--SP]
atan1	endp


Далее, нам нужен нулевой регистр i или k, чтобы обращаться к [X] (вот такая уж весёлая у нас адресация), пока, бросив монетку, это ответственное дело доверим регистру k:

	k	0


Найдём x·y, исключительно ради знака. Знак этого произведения будет соответствовать знаку результата:
	Acc		[X+1]
	MUL		[X+k]
	Inv		S


Далее, нужно совершить нечто похожее на поворот на 45°. А именно, в [SP] поместим |x|+|y|, причём в беззнаковом выражении. Максимум мог бы получиться, если x=y=1, хотя мы ожидаем вектор единичной длины на входе, и тогда самое большое, что может получиться - это 1,414. Как-то так:
	ABS	[X+k]
	ABSA	[X+1]
	[SP]	UAC	;в [SP] лежит x'

Назовём полученную штуку x'.

Далее, в [SP+1] поместим |y|-|x|. Можно не волноваться о переполнении, т.к каждый модуль по отдельности принимает значения от 0 до 1 (формат Q1.15), а их разность будет лежать в диапазоне -1..+1. Точнее, целочисленное значение -32768 (-1) может превратиться в 32767 (0,999969), но это не так страшно:

	ABS	[X+1]
	ABSS	[X+k]
	[SP+1]	Acc	;в [SP+1] лежит y'

Эту штуку назовём y'.

Здесь мы обязаны были применить Acc, чтобы значение 32768 "насытилось" до 32767 (значение +32768 не уместится в 16 бит, и будет записано как -32768, т.е ошибёмся ДИАМЕТРАЛЬНО). А в первый раз (чтобы найти x') столь же обязаны поставить UAC (Unsigned ACc), чтобы беззнаковые значения выше 32767 не "насыщались", а уходили в память "как есть".

И пора, наконец-то, начать вычисления!

Реализуем "схему Горнера". В первую очередь находим y'2:
	;начинаем вычисления
	SQRD2	Acc	;результат чуть точнее за счёт байпаса


Точнее, это будет y'2/2, и это деление на 2 надо будет учесть в коэффициенте. Мы собирались использовать -5517/32768, но поскольку мы совершили поворот с умножением на корень из двух, то здесь нужно поделить ещё на 2 (т.к y'2), а имеющееся деление на 2 "оставить":
	C		-5517
	MUL		Acc


Дело в том, что результат всего выражения в кои-то веки надо будет представить в формате Q2.14 (от -2 до +2), вместо Q1.15, и для "преобразования форматов" и нужно это дополнительное деление на 2. Снова радуемся байпасу: лишний раз точность не теряется, мы по сути умножаем 32-битный результат на 16-битное число в рамках 32-битного аккумулятора.

Далее нам нужно из результата вычесть 21779/32768x', причём x' здесь беззнаковый, от 0 до 2. Коэффициент уменьшим вдвое, чтобы прийти к формату Q2.14, и ещё в корень из двух, т.к вместо поворота на 45° мы осуществили поворот с умножением на корень из двух:

	C	10890
	UFMS	[SP]	;-5517y'^2/2-10890x'


Напомню: UFMS - это Unsigned Fused Multiply-Subtract, ровно та команда, что сейчас нужна :)

Остаётся прибавить константу, 54554/32768. Но для перехода к формату Q2.14 и она должна быть вдвое меньше, что весьма удачно для нас: специальной команды "загрузить беззнаковое значение" или "прибавить беззнаковое значение" у нас нет. Хотя скорее всего она и не нужна: разница в "скрытом" старшем бите аккумулятора, который можно тупо проигнорировать, если берёшь UAC вместо Acc. Но здесь и этого делать не придётся:

ADD	27277	;-5517y'^2/2-10890x'+27277


Скобку, наконец-то завершили. Теперь умножаем это безобразие на y':

	C	[SP+1]
	MUL	Acc


Правда, тут же соображаем, что все константы внутри скобки нужно ещё раз поделить на корень из двух, чтобы скомпенсировать тот, что здесь.

И ещё, нужно прибавить 45°, или π/4. Ведь если раньше мы брали x,y "как есть". А теперь, если x=y≈0,707, мы получим y'=0, значит результат вычислений тоже нулевой, а должен быть: 45°. Ну да, прибавим, но в формате Q2.14:

	ADD	12868


Вот, собственно, и всё.

Выпишем всю процедуру целиком:

;арктангенс через Чебышева
;x лежит в [X], y лежит в [X+1], результат выдаётся в аккумуляторе
;затирается регистр C и флаг знака, остальные остаются "как есть"
atan1 proc
	[SP++]	ijk
	k	0
	Acc	[X+1]
	MUL	[X+k]
	Inv	S

	ABS	[X+k]
	ABSA	[X+1]
	[SP]	UAC	;в [SP] лежит x'
	ABS	[X+1]
	ABSS	[X+k]
	[SP+1]	Acc	;в [SP+1] лежит y'
	;начинаем вычисления
	SQRD2	Acc	;результат чуть точнее за счёт байпаса
	C	-1951
	MUL	Acc	;-1951y'^2/2
	
	C	5445
	UFMS	[SP]	;-1951y'^2/2-5445x'
	
	ADD	19288	;-1951y'^2/2-5445x'+19288
	
	C	[SP+1]
	MUL	Acc	;y'*(-1951y'^2-5445x'+19288)
	
	ADD	12868	;y'*(-1951y'^2-5445x'+19288)+12868
	
	ijk	[--SP]
	JMP	[--SP]
atan1	endp


Впрочем, мы ещё одну вещь забыли: назначить знак результату! Тут, увы, совсем "чисто" код не получается. Замена последнего MUL на FMPM невозможна, т.к FMPM ПРИБАВЛЯЕТ ЛИБО ВЫЧИТАЕТ результат умножения из ЗНАЧЕНИЯ В АККУМУЛЯТОРЕ. А нам бы хотелось перед FMPM инициализировать его нулём, но тогда его текущее значение надо "эвакуировать". Например, можно содержимое Acc перенести в 16-битный C, затем обнулить Acc и вызвать, наконец, FMPM:

;	C		[SP+1]
;	MUL		Acc	;y'*(-1951y'^2-5445x'+19288)
	C		Acc
	ZAcc		RoundZero
	FMPM		[SP+1]


Из хорошего: теперь у нас корректно выполнится округление, за счёт того, что к финальному результату прибавим 1/2 от единицы младшего разряда. Из плохого: младшие биты будут утеряны, т.е часть точности у нас уйдёт в трубу.

(потом выбрать, прибавлять ли π/4 или отнимать - не проблема, заменяем ADD на PM).

Более точный вариант: дойти до финального результата (как в приведённом коде процедуры), положить его в тот же C, обнулить аккумулятор и потом выполнить команду PM C:

	C		Acc
	Acc		0
	PM		C


Так длиннее, зато точности терять не должны, т.к финальный результат всё равно будет обрезан до 16 бит. И ещё одну вещь, наверное, надо сделать. Чтобы и в таком варианте округление было правильным, нужно вместо ADD 12868 написать DIV2A 25737 (прибавить значение, поделённое на 2). Итого, вот так:

;арктангенс через Чебышева
;x лежит в [X], y лежит в [X+1], результат выдаётся в аккумуляторе
;затирается регистр C и флаг знака, остальные остаются "как есть"
atan1 proc
	[SP++]	ijk
	k	0
	Acc	[X+1]
	MUL	[X+k]
	Inv	S

	ABS	[X+k]
	ABSA	[X+1]
	[SP]	UAC	;в [SP] лежит x'
	ABS	[X+1]
	ABSS	[X+k]
	[SP+1]	Acc	;в [SP+1] лежит y'
	;начинаем вычисления
	SQRD2	Acc	;результат чуть точнее за счёт байпаса
	C	-1951
	MUL	Acc	;-1951y'^2/2
	
	C	5445
	UFMS	[SP]	;-1951y'^2/2-5445x'
	
	ADD	19288	;-1951y'^2/2-5445x'+19288
	
	C	[SP+1]
	MUL	Acc	;y'*(-1951y'^2-5445x'+19288)
	; C	Acc
	; ZAcc	RoundZero
	; FMPM	[SP+1]
	
	;ADD	12868	;y'*(-1951y'^2-5445x'+19288)+12868
	DIV2A	25737	;ещё 1/2 младшего разряда для правильного округления
	
	C	Acc
	Acc	0
	PM	C
	
	ijk	[--SP]
	JMP	[--SP]
atan1	endp



Итого, 26 слов кода (1 лишнее из-за Hazard на одновременное чтение и запись памяти, транслятор вставляет NOP), или 52 байта. Не используется сегмент данных (только 3 значения хранятся на стеке). То есть, вышло чуточку компактнее предыдущей программы, которая занимала 29 слов кода и 15 слов данных.

Компилируется без проблем, вообще редчайшая вещь у нас - процедура совсем без ветвления! А испытывать начнём в следующем посте.
Tags: ПЛИС, математика, программки, работа, странные девайсы
Subscribe

  • Великая Октябрьская резня бензопилой

    Сегодня прокатился прочистить Абрамцевскую просеку. Как обычно, с приключениями. Выезжал на велосипеде, а вернулся на самокате. Первый раз по этим…

  • Очередная несуразность в единицах измерения

    Когда-то я написал программу PhysUnitCalc - калькулятор, умеющий работать с размерностями. Мне казалось, что я наступил уже на все грабли, которые…

  • Big Data, чтоб их... (3)

    "В предыдущих сериях": мой прибор выдаёт 6 значений: 3 координаты и 3 угла, т.е все 6 степеней свободы твёрдого тела. Причём ошибки измерения этих 6…

  • Post a new comment

    Error

    Anonymous comments are disabled in this journal

    default userpic

    Your IP address will be recorded 

  • 0 comments