nabbla (nabbla1) wrote,
nabbla
nabbla1

Categories:

Громко думаем о знаках для FMPM

Это та часть, из-за которой QuatCore (Quaternion core) так называется! Мне хотелось, чтобы умножение кватернионов записывалось довольно просто, без необходимости "в лоб" привести все 16 умножений и 12 сложений, как это делается обычно. А также не хотел делать отдельной процедуры для умножения на сопряжённый кватернион, или того хуже, выделять память под ещё один кватернион, копировать туда исходный, брать сопряжение - и только после этого умножать, поэтому и придумал регистр Inv, который к тому же выражает вполне себе понятную нотацию "плюс-минус" и "минус-плюс".

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

Изначально таблица была получена из следующих выкладок для перемножения кватернионов:

:









Если слагаемые справа расположить в порядке μ0, μx, μy, μz, то и выйдет наша нынешняя таблица:


Но это не единственный вариант...


Ведь мы можем расположить их в порядке λ0, λx, λy, λz и получить:


Эта выглядит более симметричной: в ней поменять все знаки - всё равно, что транспонировать! Да, кажется, с неё я и начинал, а потом "засомневался" - получил неправильный результат, а всего лишь потому, что неправильно операнды взял. Тут вся разница в том, по какому из индексов устроить цикл, а какой брать как XOR между номером строки и столбца.

А вот для арктангенса эта таблица кажется более предпочтительной! Только там бы заменить j на k:


Внешний цикл будет идти по j от 4 до 0 - будут перебираться элементы таблицы AtanTable, где лежат числа exp(i*pi/2), exp(i*pi/4), exp(i*pi/8), exp(i*pi/16) и exp(i*pi/32)*crazy_const. Но там же лежат углы в радианах, которые мы добавляем или отнимаем на каждой итерации: 0, pi/4, pi/8, pi/16 и pi/32. Поскольку на каждую итерацию требуется по 3 числа (действительная и мнимая часть комплексного числа, и ещё угол), то надо использовать адресацию "3j", это вынуждает использовать в точности регистр j, только для него такое доступно:


И теперь наконец-то все части встают на место: по [Z+3j] я буду хранить угол, по [Z+3j+1]: действительную часть и по [Z+3j+2]: мнимую. Получится вот такая процедура:

;atan1(y,x)
;берёт atan(y/x), результат от -pi/2 до pi/2, в формате Q2.12 (т.е представимые углы от -2 до +2)
;также находит длину вектора (x,y). 

;в [X], [X+1] лежат значения x,y, их трогать не будем
;в [Y] после выполнения будет лежать длина вектора, делёная на 2, а в [Y+1]: угол. 

;точность нахождения угла: все 16 бит, если (x,y) нормированы (т.е x^2+y^2=1), в противном случае не хуже 5,625 градусов (5..6 бит значащих) 
;точность нахождения длины вектора: 0,36% (8..9 значащих бит)

atan1 proc

		; [SP++]	X	;будет изменён в процессе
		; [SP++]	Z	;будет указатель на таблицу углов
		; [SP++]	ijk	;все регистры в хвост и в гриву
		; [SP++]	C
					; ;а вот Z оставим нетронутым! Ну и Acc считаем "временным", равно как и флаги
			
;assign ijk = {invreg, kreg, ireg, jreg};
		ijk		0x0804		;Inv=0, k=2, i=0, j=4
		ZAcc		RoundZero
		SUB		[X+i]		;то есть просто [X]
		[SP]		0		;угол инициализируем нулём
		Z		AtanTable
@@j_loop:	Inv		S
		C		[X+i]		;то есть просто [X], здесь всегда i=0
		i		1	
@@i_loop:	ZAcc		RoundZero
		FMPM		[Z+3j+k]	;то есть [Z+3j+2], он же "синус" (мнимая часть). т.к k=2 на всём протяжении
		C		[X+i]
		FMA		[Z+3j+1]	;он же "косинус" (действительная часть)
		[Y+i]		Acc
		iLOOP		@@i_loop	
;вышло 8 строчек, но пока фигня с j (когда не равно единице, со знаками песец выходит)
		X		Y	;возможно, X сохраним в стек, чтобы процедура минимально меняла регистры
		Acc		[SP]
		PM		[Z+3j+i]	;то есть [Z+3j], угол. 
		[SP]		Acc		
		jLOOP		@@j_loop
		ADD		[Y+i]
		[Y+i]		Acc
		
		; C		[--SP]
		; ijk		[--SP]
		; Z		[--SP]
		; X		[--SP]

		JMP		[--SP]
atan1 endp



Наконец-то всё как надо: одной строкой "i 2" стало меньше (таким грубым способом мы пытались обратиться к [Z+3j+2], т.е получить угол) - теперь именно в этом месте у нас i=0, а ровно так нам и надо. k=2 нужен, чтобы получить нужное нам чередование знаков - один "синус" со знаком плюс, второй - со знаком минус, а если Inv=1 - то всё наоборот, в общем, обведённое синим место в таблице.

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

А что там с умножением кватернионов
Удивительная процедура - она была написана самой первой, потом несколько раз переписывалась, но так до сих пор не пошла в ход! Она пригодится в алгоритме сопровождения, а мы пока алгоритм захвата мучали, без которого никакого сопровождения не видать...

Последний вариант был такой:
QuatMultiply proc
		[SP++]		C
		[SP++]		i 			;заносим в стек все индекс. регистры разом
		[SP++]		j
		[SP++]		k
		j 		3			;j-номер компонента, который вычисляем в данный момент
		i		3    		;k-номер слагаемого
@@row_loop:	ZAcc		RoundZero	;обнуляем аккумулятор
@@col_loop:	C		[X+i^j]		;загрузили левый операнд
		FMPM		[Y+i] 		;помножили на правый, с нужным знаком
		iLOOP		@@col_loop
		k		j
		[SP+k]		Acc 		;храним ответы на стеке
;потому что если X=Z или X=Y (обновляем кватернион),
;то занесение сразу в Z исказит последующие вычисл.
		i		3
		jLOOP		@@row_loop
@@out_loop:	[Z+i]		[SP+i] 	;помещаем ответ куда надо
		iLOOP		@@out_loop
		k		[--SP]
		j		[--SP]
		i		[--SP]		;восстанавливаем исх. значения регистров
		C		[--SP]
		JMP		[--SP]		;возврат из процедуры
QuatMultiply endp


Мучительное запихивание в стек всех индексных регистров - да, в тот момент меня жаба душила включить команду ijk... И на ровном месте танцы с бубном, с присваиванием k=j, потому что адресация [SP+k] есть, а [SP+j] - нет, только [SP+2j]. То есть, здесь тоже замена регистра в таблице, с j на k, идёт на пользу!

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

Выходит, код надо исправить на вот такой:
QuatMultiply proc
		[SP++]	C
		[SP++]	ijk 			;заносим в стек все индекс. регистры разом
		;нельзя использовать ijk, т.к мы перезапишем Inv.
		k 		3		;k-номер компонента, который вычисляем в данный момент
		i		3    		;i-номер слагаемого
@@row_loop:	ZAcc		RoundZero	;обнуляем аккумулятор
@@col_loop:	C		[X+i]		;загрузили левый операнд
		FMPM		[Y+i^k] 		;помножили на правый, с нужным знаком
		iLOOP		@@col_loop
		[SP+k]		Acc 		;храним ответы на стеке
;потому что если X=Z или X=Y (обновляем кватернион),
;то занесение сразу в Z исказит последующие вычисл.
		i		3
		jLOOP		@@row_loop
@@out_loop:	[Z+i]		[SP+i] 	;помещаем ответ куда надо
		iLOOP		@@out_loop
		ijk		[--SP]		;восстанавливаем исх. значения регистров
		C		[--SP]
		JMP		[--SP]		;возврат из процедуры
QuatMultiply endp


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

Но была же причина вводить именно i^j?
Да, это я пытался реализовать умножение матрицы 2х2 на матрицу поворота, "на ходу" сформированную из синуса и косинуса, злосчастная процедура FindRoll:
;здесь i неизвестен, j=1, k=Inv=0
;внешний цикл по j, номер строки результата, и строки co/si
@@j_loop:	k		1	;номер столбца результата, и столбца AfTransf
@@k_loop:	i		1	;номер столбца co/si и строки AfTransf
		ZAcc		RoundZero	;обнулить до 1/2 мл. разр
@@i_loop:	C		[Z+i^j]	
		FMPM		[Y+2i+k]
		iLOOP		@@i_loop
		X		AfTransf	;сюда результат складируем				
		[X+2j+k]	Acc
		kLOOP		@@k_loop
		jLOOP		@@j_loop


Именно здесь, чтобы получить "двойную" адресацию результата (строка+столбец), нужно было, чтобы в таблице одним из регистров был j.

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

Мало того, здесь нам вовсе не обязательно, чтобы поворот осуществлялся в одну сторону, если Inv=0, и в другую, если Inv=1. Тут у нас зафиксировано Inv=0, а потому можно использовать вот эту область:


Так что будет два цикла, один, внешний - по j, внутренний - по k, а регистр i на всём протяжении должен оставаться нулевым.


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

  • 0 comments