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

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

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

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

UPD:
вот ещё вариант со знаком. Первый сомножитель — var20:var21, число как прочиталось из АЦП со знаком. Второй сомножитель — в памяти, 32бита со знаком.
На выходе в var12:var13 модуль результата, и отдельно в T знак результата, таким образом достигается диапазон выходных значений +-65535

;**********************************************************************************************
; Signed conversion procedure
;**********************************************************************************************
Muladdsh:	ld	var10, y+	
		ld	var11, y+
		ld	var12, y+
		ld	var13, y+
; Проверка будущего знака результата умножения
		eor	var21, var13	; сравниваем знаковые биты
		bst	var21, 7	; запоминаем во флаге T знак будущего частного
		eor	var21, var13	; восстанавливаем var21
; получение модуля множителей
		sbrs	var21, 7	; Проверка знака конвертируемого числа
		rjmp	_mas_v1pos	; если конвертируемое отрицательно
		ldi	r16, 0		; преобразуем его 
		ldi	r17, 0		; в положительное
		sub	r16, var20	; методом вычитания
		sbc	r17, var21	; из ноля
		movw	var20, r16
_mas_v1pos:	sbrs	var13, 7	; Проверка знака коэффициента K
		rjmp	_mas_v2pos	; Конверсия в положительный при необходимости
		ldi	r16, 0
		sub	r16, var10
		mov	var10, r16
		ldi	r16, 0
		sbc	r16, var11
		mov	var11, r16
		ldi	r16,0
		sbc	r16, var12
		mov	var12, r16
		ldi	r16,0
		sbc	r16, var13
		mov	var13, r16
_mas_v2pos:	rcall	Mul24x16uh	; Var1{01234} = var1{012} * var2{01}
		rol	var11
		rcall	round24		; Округление после деления
; Загрузка офсета		
		ld	var20, y+
		ld	var21, y+		
; Ветвление, исходя из знака результата умножения
; Возможны 4 случая:
; 1) положительный+положительный. возможно переполнение сверху при умножении или сложении
; 2) положительный+отрицательный. Возможно переполнение сверху при умножении, или переход через ноль при сложении
; 3) отрицательный+положительный. возможно переполнение снизу на этапе умножения, или переход через ноль при сложении
; 4) отрицательный+отрицательный. переполнение снизу при умножении или сложении.
		ldi	r16,0
		brtc	_mas_repos	; проверка - должен ли результат быть отрицательным?
; здесь - результат умножения отрицательный
		sbrs	var21, 7
		rjmp	_mas_negpos
; здесь - отрицательный результат умножения и отрицательный офсет. Возможно переполнение снизу.
		clr	var10
		clr	var11
		sub	var10, var20	; меняем знак офсета
		sbc	var11, var21
		add	var12, var10
		adc	var13, var11
		adc	var14, r16
		brcs	_sicon_max	; перенос после сложения
_sicnchkovf:	tst	var14
		breq	_sicon_exit	; проверка на переполнение
_sicon_max:	ldi	r16, 0xff	; Ограничение сверху/снизу
		mov	var12, r16
		mov	var13, r16
		rjmp	_sicon_exit
; здесь - отрицательный результат умножения и положительный офсет. Возможно переполнение снизу(сверху по модулю) или изменение знака.
_mas_negpos:	sub	var12, var20
		sbc	var13, var21
		sbc	var14, r16
		brcc	_sicnchkovf	; проверка перехода через 0, смена знака
; здесь - произошел переход через 0 при вычитании. результат имеет противоположный знак.
		clt			; смена знакового флага
_sicnchsgn:	sub	r16, var12
		mov	var12, r16	; смена знака результата
		ldi	r16, 0
		sbc	var13, r16
		mov	var13, r16
		rjmp	_sicon_exit
; Положительный результат умножения
_mas_repos:	sbrs	var21, 7	; проверка знака офсета
		rjmp	_mas_pospos
; Здесь - положительное число после умночения и отрицательный офсет. Возможно переполнение, возможен переход через 0
		clr	var10
		clr	var11
		sub	var10, var20	; меняем знак офсета
		sbc	var11, var21
		sub	var12, var10	; отнимаем от произведения модуль офсета
		sbc	var13, var11
		sbc	var14, r16
		brcc	_sicnchkovf	; нет займа, знак не сменился, проверяем на переполнение
; смена знака результата после сложения. Результат отрицательный, и ложно обозначен как положительный
		set			; меняем знаковый флаг
		rjmp	_sicnchsgn	; меняем знак результата	
_mas_pospos:	
; Здесь - положительное число после умночения и положительный офсет. Возможно только переполнение.
		add	var12, var20
		adc	var13, var21
		adc	var14, r16
		brcs	_sicon_max	; перенос после сложения
		brne	_sicon_max	; переполнение после сложения
_sicon_exit:	ret

И процедура ассиметричного умножения 24*16 бит для неё:

;***************************************************************************
; Var1x{01234} = var1x{012} * var2x{01}
;***************************************************************************
Mul24x16uh:	ldi	r16,0
		mul	var12, var21	; b*x
		mov	var13, r0
		mov	var14, r1

		mul	var12, var20	; b*y
		MOV	var12, R0
		add	var13, R1
		adc	var14, r16	; zero

		mul	var11, var21	; c*x
		add	var12, R0
		adc	var13, R1
		adc	var14, r16	; zero
		
		mul	var11,var20	; c*y
		mov	var11,R0
		add	var12,R1
		adc	var13, r16	; zero
		adc	var14, r16	; zero

		mul	var10,var21	; d*x
		add	var11,R0
		adc	var12,R1
		adc	var13, r16	; zero
		adc	var14, r16	; zero

		mul	var10, var20	; d*y
		mov	var10, R0
		add	var11, R1
round24:	adc	var12, r16	; zero	- и повторный вход для выполнения округления при необходимости
		adc	var13, r16	; zero
		adc	var14, r16	; zero
		ret