Wurzel ziehen bei Mikrocontrollern

Eine schnelle Quadratwurzel wird recht oft benötigt:

• Effektivwertberechnung (Spannung, Strom) nach Abtastung und Summation der Quadrate
• Entfernungsberechnung zwischen zwei (auch räumlichen) Punkten
• Statistik (Standardabweichung u.ä.)
• Grafikfunktionen (Kreise, Bögen und Ellipsen)

Grundlagen

Folgende Lösungsmöglichkeiten sind üblich und gängig:

Hinweis: Keines dieser Verfahren benötigt die Multiplikation und ist deshalb auch auf weniger leistungsstarken Mikrocontrollern verwendbar.

Summe der ungeraden Zahlen

Eine Tabelle soll jedoch noch einmal die Langsamkeit dieser Routine verdeutlichen:

Tabellenzugriff

2ⁿ * 2n/8 Bytes

Aber - für popelige Anzeigen gibt es ein besseres Verfahren, s.u.

QuadTab:	;!UNGETESTET!
.set i = 0
	.rept 256

	.dw	i*i
.set i = i+1
.endm

;Löwenfangmethode: Wüste halbieren...
	ldiw	zh,zl,(QuadTab+128)*2	;Zeiger auf ersten Testwert (Tabellen-Mitte)
	ldi	r16,128			;Distanz (zum nächsten Testwert)
halbieren:		;8 Runden
	lpm	r0,z+
	lpm	r1,z
	dec	zl

	cp	r0,r3
	cpc	r1,r4
	
	brcc	h1

	sub	zl,r16
	subi	zh,0
	rjmp	h2
h1:
	add	zl,r16
	adci	zh,0
h2:
	srl	r16			;Distanz halbieren
	brne	halbieren		;Distanz=0 wenn fertig
	srl	zh
	ror	zl		;Wortadresse
	subi	zl,LOW(QuadTab)	;jetzt steht in ZL der ganzzahlige WurzelWert
				;und in R1:R0 (immer noch) das Quadrat dazu
	or	r16,zl

Lösung wie schriftlich, aber binär

Die folgende Routine für AT90, ATmega8 ist für 24-bit-Radikand und getestet. Sie ist sogar schneller als der A/D-Wandler:-)

Man beachte für 8- oder 16-bit-Ergebnisse, dass die Doppelwurzel ein Bit mehr erfordert. Der Akku braucht sogar zwei Bit mehr, dies kann mit dem Übertragsflag abgefangen werden!
Daher ist diese Routine für genau 8 oder 16 bit Ergebnisbreite aufzublasen.

.def RL = r2	;Radikand (24 bit), Datenquelle
.def RM = r3
.def RH = r4
.def AL = r12	;Akku (14 bit), Rest (13bit) am Ende
.def AH = r13
.def WL = r18	;Ungerade Doppelwurzel (13bit), Wurzel (12bit) am Ende
.def WH = r19
.def LoopCounter = r17	;Schleifenzähler

;ca. 280 Takte (17,5 µs @ 16 MHz)
Wurzel:
;Doppelwurzel = 0 (LSB wird sowieso gesetzt)
	clr	WL
	clr	WH
;Akku löschen
	clr	AL
	clr	AH

	ldi	LoopCounter,12	;Für jedes Ergebnis-Bit eine Runde
zieh:
;Doppelwurzel "verdoppeln"
	sec
	rol	WL	;[xxxxx1 -> xxxxx01]
	rol	WH
	cbr	WL,2
;Zwei Bits aus Radikand in Akku schieben
	lsl	RL
	rol	RM
	rol	RH
	rol	AL
	rol	AH
	lsl	RL
	rol	RM
	rol	RH
	rol	AL
	rol	AH
;	brcs	setbit	;auskommentieren für 22/30-bit-Radikand
;Akku mit Doppelwurzel vergleichen
	cp	AL,WL
	cpc	AH,WH
	brcs	okay
setbit:
;Doppelwurzel vom Akku abziehen
	sub	AL,WL
	sbc	AH,WH
;Doppelwurzel+=2
	sbr	WL,2	;[xxxxx01 -> xxxxx11]
okay:
;nächste Runde
	dec	LoopCounter
	brne	zieh
;Doppelwurzel halbieren
	lsr	WH
	ror	WL
;fertig
	ret

;PE: BC:IY = Radikand 30 bit
;PA: DE = Wurzel 15 bit
;    HL = Rest 16 bit
;VR: IY=BC=A=0, HL, DE
;Ta: ca. 3300 (1,9 ms @ 1,75 MHz)
Wurzel:
;Doppelwurzel = 0 (LSB wird sowieso gesetzt)
	xor	a	;(4)
	ld	e,a	;(4)
	ld	d,a	;(4)
;Akku löschen
	ld	l,a	;(4)
	ld	h,a	;(4)
	ld	a,16	;(7)
zieh:
;Doppelwurzel "verdoppeln"
	scf		;(4)
	rl	e	;(8) [xxxxx1 -> xxxxx01]
	rl	d	;(8)
	res	e,1	;(8)
;Zwei Bits aus Radikand in Akku schieben
	add	iy,iy	;(15)
	rl	c	;(8)
	rl	b	;(8)
	adc	hl,hl	;(15)
	add	iy,iy	;(15)
	rl	c	;(8)
	rl	b	;(8)
	adc	hl,hl	;(15)
	jc	setbit	;(12/7)
;Akku mit Doppelwurzel vergleichen
	sbc	hl,de	;(15)
	add	hl,de	;(11)
	jc	okay	;(12/7)
setbit:
;Doppelwurzel vom Akku abziehen
	or	a	;(4)
	sbc	hl,de	;(15)
;Doppelwurzel += 2
	set	e,1	;(8) [xxxxx01 -> xxxxx11]
okay:
;nächste Runde
	dec	a	;(4)
	jnz	zieh	;(12/7) => (203) *16 = (3248)
;Doppelwurzel halbieren
	srl	h	;(8)
	rr	l	;(8)
;fertig
	ret		;(10)

Mischtechnologien

Wegen der Effektivität der „schriftlichen“ Routine ist diese Möglichkeit ins Hintertreffen geraten.

;Deshalb hier eine Mischtechnologie für den 24-bit-Input:
;* Das QH wird genommen und die Wurzel bestimmt (0..15)
;* Ausgehend vom Quadrat der Wurzel * 256 sowie der Wurzel * 16
;  wird die Wurzel aus QH:QM bestimmt, das sind wieviel Iterationen max.??
;Beispiel:
;  FF² = FE01
;  sqrt(FE) = F (Wurzel aus High-Byte ziehen)
;  F²  = E1 (Nebenprodukt des Wurzelziehens)
;  F0² = E100 (zusammensetzen, weil 16² = 256)
;  F1² = E2E1 (addiere 2E1 = 2*F0+1)
;  das Beispiel zeigt: max. 16 weitere Iterationen, allerdings 16-bit-Strichrechnung
;Weiter geht's dann genauso mit 24-bit-Strichrechnung

Andere Quellen

• Halbierung des Radikanden, damit auch verringerter Aufwand beim Addieren (?)
• Mathematische Rundung des Ergebnisses
• Vermeidung des Wurzelziehens beim Berechnen des Betrags einer komplexen Zahl
• Eine schnelle Näherungsroutine wäre das Heron-Verfahren.