Сложить и приумножить

Дискламер: Тут будет описано несколько боянистых боянов. Поэтому, если вы решите что боян, можно не читать, скажете что я разрешил. А на случай если кто-то не знает, пусть это лежит здесь.

В противополжность персонажам известного анекдота, сегодня мы не станем отнимать и делить. Напротив, мы будем умножать и складывать, чтобы получить в результате всем известное уравнение прямой y=kx+b.
В свою очередь, уравнение прямой нам ценно не само по себе, а поскольку оно позволяет строить калибровочные характеристики, например чтобы перевести коды АЦП в микроамперы. Или наоборот — формула одна и та же.

Начнём с умножения. Существует как минимум два способа умножать. Первый, под кодовым наименованием «сдвинуть и сложить» описан в AVR200, и в учебниках по программированию. Достоинство этого способа в простоте и лёгкой масштабируемости. Если захочется умножить два 512-битных числа, то этот метод для вас. Если захочется умножить эти же два числа на процессоре типа mcs48, или i8080, или на пике — добро пожаловать! Работает на всём, что умеет в сложение и сдвиг. Единственный недостаток — работает сравнительно медленно.
Другой вариант работает только там, где есть аппаратный умножитель. Довольно косноязычно этот способ описан в AVR201, но я сейчас расскажу немного подробнее.
Для начала пара слов о размерности сомножителей — обычно я использую коэффициент крутизны длиной 32 бита. 16 бит на целую часть, 16 бит на дробную, это позволяет оперировать с 16ти битными данными ЦАП или АЦП с минимальными потерями точности на округление. И получить такое число проще простого — если у вас есть рассчитанный, скажем, в экселе коэффициент 0.1488 милливольта на бит АЦП — надо просто умножить его на 65536 и округлить. Получится 9752.
Итак, умножаем шестнадцатибитное число на тридцатидвухбитное — и получаем в итоге 48 бит. Из них нам потребуются всего 16, но остальные тоже будут использованы.
Всем тем, кто в школе пропустил тот урок, когда учили умножать в столбик, или кто с тех пор благополучно всё забыл (как я например), я подготовил гифку-залипалку. Красными стрелочками показано умножение, а зелёными — сложение.


Метод хорош, но требует различных регистров для результата и сомножителей, тогда как после умножения один из сомножителей обычно не нужен. Поэтому вот вариант, оптимизированный на эконмию регистров ценой снижения быстродействия.


Теперь отливаем в камне:

.def var10=r2
.def var11=r3
.def var12=r4
.def var13=r5
.def var14=r6
.def var15=r7
.def var16=r8
.def var17=r9
.def var20=r10
.def var21=r11
.def var22=r12
.def var23=r13

; ***************************************************************************
; multiplication - addiction unsigned uses hadware
; ассиметричные сомножители: var2x = 16bit
; array[Y] = 32 бита множитель + 16 бит слагаемое
; результат var12 var13
; ***************************************************************************
Muladduh:	clr	var16		;zero
		ldd	var10,y+3
		mul	var10,var21	; A*X
		movw	var14,r0
		ldd	var10,y+2
		mul	var10,var21	; B*X
		mov	var13,r0
		add	var14,r1
		adc	var15,var16	;+ zero
		ldd	var10,y+3
		mul	var10,var20	; A*Y
		add	var13,r0
		adc	var14,r1
		adc	var15,var16	;+ zero
		ldd	var10,y+1
		mul	var10,var21	; C*X
		mov	var12,r0
		add	var13,r1
		adc	var14,var16	;+ zero
		adc	var15,var16	;+ zero
		ldd	var10,y+2
		mul	var10,var20	; B*Y
		add	var12,r0
		adc	var13,r1
		adc	var14,var16	;+ zero
		adc	var15,var16	;+ zero
		ldd	var10,y+0
		mul	var10,var21	; D*X
		mov	var11,r0
		add	var12,r1
		adc	var13,var16	;+ zero
		adc	var14,var16	;+ zero
		adc	var15,var16	;+ zero
		ldd	var10,y+1
		mul	var10,var20	; C*Y
		add	var11,r0
		add	var12,r1
		adc	var13,var16	;+ zero
		adc	var14,var16	;+ zero
		adc	var15,var16	;+ zero
		ldd	var10,y+0
		mul	var10,var20	; D*Y
		mov	var10,r0
		add	var11,r1
		adc	var12,var16	;+ zero
		adc	var13,var16	;+ zero
		adc	var14,var16	;+ zero
		adc	var15,var16	;+ zero

Вы ещё помните, что в первом сомножителе у нас были 16 бит дробной части? После умножения теперь и результат содержит 16 бит дробной части. Если просто отбросить их, то мы ухудшим точность вычислений. Поэтому будем округлять к ближайшему целому. Как вы наверное помните, при десятичном округлении, пять и выше в отбрасываемом разряде увеличивает остающееся число на единицу. Возможно, вы будете удивлены, но с двоичными числами поступают точно так же. Если старший бит из отбрасываемой группы — единица, то она добавляется к округляемому числу.

; округление при делении на 65536
		lsl	var11; выталкиваем старший бит в перенос
		adc	var12,var16	;+ zero
		adc	var13,var16	;+ zero
		adc	var14,var16	;+ zero
		adc	var15,var16	;+ zero

Начиная с этого момента, младшие 16 бит нам больше не нужны, они уже сыграли свою роль.

Теперь складываем со вторым коэффициентом прямой — со смещением (тот самый b в уравнении). В качестве смещения я обычно использую 16 бит целое со знаком — это разумно большое и достаточное число.
Здесь нам потребуется ветвление, чтобы правильно ограничить результат сверху и снизу. Именно для этого момента мы и хранили старшие два байта результата — если сделать ограничение сверху сразу после умножения — то при отрицательном смещении будет потеряна часть диапазона результата. Поэтому ограничивать надо только после завершения всех арифметических операций. Если смещение положительно, то просто прибавляем его к результату умножения, и проверяем на ограничение сверху.
А если смещение отрицательное — вычитаем модуль смещения из результата умножения, проверяем на ограничение снизу, а затем дополнительно проверяем на ограничение сверху. Только таким образом у нас получится максимально использовать весь доступный диапазон выходного числа.

; Сложение/вычитание 
		ldd	var10,y+4	; Загружаем офсет
		ldd	var11,y+5
		sbrc	var11,7		; Проверяем знак оффсета
		rjmp	_mahso
; ветвь - офсет положительный, складываем
		add	var12, var10
		adc	var13, var11
		adc	var14, var16	;+ zero
		adc	var15, var16	;+ zero
_mahchkm:	or	var14, var15	; проверка переполнения
		breq	_mahexit
		ldi	r16,0xff	; при переполнении заполняем результат 
		mov	var12,r16	; максимальным значением
		mov	var13,r16	
_mahexit:	ret
; ветвь - офсет отрицательный
_mahso:		clr	var20
		clr	var21
		sub	var20,var10
		sbc	var21,var11
		sub	var12,var20	; вычитаем офсет из результата
		sbc	var13,var21
		sbc	var14,var16	; с учётм следующего extra байта
		brcc	_mahchkm	; если не было займа - то проверяем на возможное переполнение сверху.
		clr	var12		; если был займ после экстра-байта, значит имеем значение ниже нуля
		clr	var13
		rjmp	_mahexit

Спасибо за внимание.