Atmel Quadratwurzel

5.2.4 Quadratwurzel


Achtung alle folgenden Assembler Routinen sind eine Abbildung der Delphi Routinen. Sie wurden noch nicht optimiert.

5.2.4.1 Version 1

Für einige Funktionen wird manchmal die Quadratwurzel einer Zahl benötigt. Da der Kontroller in der Regel nicht über einen entsprechende Befehl verfügt muß diese Funktion nachgebildet werden. Das Iterationsverfahren nach Heron kommt hier nicht in Frage, da eine Division erforderlich ist. Hier wird ein Verfahren gezeigt, das nur einfache Subtraktion benutzt. Eine sehr schöne Beschreibung dieses Verfahren wird unter http://www.diaware.de/html/wurzel.html beschrieben und soll hier nicht näher erläutert werden. Im rechten Bild ist nur ein Beispiel beschrieben, wie es im hexadezimalen Bereich aussieht.

Dieses Bild zeigt den relativen Fehler in % im Vergleich zwischen der Float und der Integer Operation.

Je größer die Quadratwurzel, je kleiner der Fehler. Wird das Ergebnis der Float Operation gerundet und mit dem Ergebnis der Integer Operation verglichen, ist der Unterschied maximal 1 Bit.

Delphi Routine

Definitionsbereich = 0 ... 2147483647 = MaxInt
Sqrt(2147483647) = 46340,95
Function IntSqrt(w1:Integer):Integer;
Var i0,i1,i2,i3,i4,i5,k0,k1:Integer;
Begin
  k1 := 0;
  i1 := 0;
  i5 := 32;
  For i0 := 1 to 4 do Begin
    i5 := i5 - 8;
    i4 := (w1 shr i5) and $000000FF;
    i1 := (i1 shl 8) or i4;
    i2 := (k1 shl 5) + 1;
    k0 := 0;
    Repeat
      i3 := i1 - i2;
      If i3 < 0 Then Begin
        Break;
      End Else Begin
        i1 := i3;
        i2 := i2 + 2;
      End;
      k0 := k0 + 1;
    Until false;
    k1 := (k1 shl 4) + k0;
  End;
  Result := k1;
End;

Assembler Routine für AVR

Definitionsbereich = 0 ... 4294967295
IntSqrt(4294967295) = 65535

Laufzeiten:
SQRT($00000000) = $0000 :   740 Taktzyklen (187µs bei 4MHz)
SQRT($48D0C840) = $8888 : 1508 Taktzyklen (377µs bei 4MHz)
SQRT($FFFFFFFF) = $FFFF : 2180 Taktzyklen (545µs bei 4MHz)

Codelänge: 168 Bytes

;****************************************************************
; Quadratwurzel
;  Input: r16,r17 = Low,High Byte 1
;         r18,r19 = Low,High Byte 2
; Output: r16,r17 = Low,High Byte Quadratwurzel
; Return: -
; Modify: Alle Register
SQRT:       clr   k1a        ;k1 := 0
            clr   k1b
            clr   i1a        ;i1 := 0
            clr   i1b
            clr   i1c
            clr   i1d
            ldi   i5,32      ;i5 := 32
            ldi   i0,4       ;For i0 := 1 to 4
sqrt01:     subi  i5,8       ;i5 := i5 - 8
            mov   i4a,w1a    ;(w1 shr i5) and $000000FF
            mov   i4b,w1b
            mov   i4c,w1c
            mov   i4d,w1d
            mov   ix,i5
            tst   i5
            breq  sqrt02a    ;kein Shift wenn i5 = 0
sqrt02:     lsr   i4d
            ror   i4c
            ror   i4b
            ror   i4a
            dec   ix
            brne  sqrt02
sqrt02a:    mov   i1d,i1c    ;i1 := (i1 shl 8) or i4
            mov   i1c,i1b
            mov   i1b,i1a
            mov   i1a,i4a
            clr   i2d        ;i2 := (k1 shl 5) + 1;
            clr   i2c
            mov   i2b,k1b
            mov   i2a,k1a
            ldi   ix,5
sqrt03:     lsl   i2a
            rol   i2b
            rol   i2c
            rol   i2d
            dec   ix
            brne  sqrt03
            ldi   ixa,1
            clr   ixb
            clr   ixc
            clr   ixd
            add   i2a,ixa
            adc   i2b,ixb
            adc   i2c,ixc
            adc   i2d,ixd
            clr   k0         ;k0 := 0;
            nop              ;Repeat
sqrt04:     mov   i3a,i1a    ;i3 := i1 - i2
            mov   i3b,i1b
            mov   i3c,i1c
            mov   i3d,i1d
            sub   i3a,i2a
            sbc   i3b,i2b
            sbc   i3c,i2c
            sbc   i3d,i2d
            brmi  sqrt05     ;Exit Repeat If i3 < 0
            mov   i1a,i3a    ;i1 := i3
            mov   i1b,i3b
            mov   i1c,i3c
            mov   i1d,i3d
            ldi   ixa,2
            clr   ixb
            clr   ixc
            clr   ixd
            add   i2a,ixa    ;i2 := i2 + 2
            adc   i2b,ixb
            adc   i2c,ixc
            adc   i2d,ixd
            inc   k0         ;k0 := k0 + 1
            rjmp  sqrt04
sqrt05:     ldi   ix,4       ;k1 := (k1 shl 4) + k0
sqrt06:     lsl   k1a
            rol   k1b
            dec   ix
            brne  sqrt06
            add   k1a,k0
            adc   k1b,ixb
            dec   i0
            breq  sqrt07
            rjmp  sqrt01
sqrt07:     clr   w1d        ;Result := k1 
            clr   w1c
            mov   w1b,k1b
            mov   w1a,k1a
            ret

5.2.4.2 Version 2

Es gibt noch eine andere Möglichkeit. Sie hat wesentlich weniger Code ist dafür aber erheblich langsamer je größer die Quadratzahlen sind.

Von der Quadratzahl wird bei jedem Durchgang eine um 2 erhöhte Zahl subtrahiert (beginnend bei 1) bis das Ergebnis <= 0 ist, also -1, -3, -5, -7, -9, -11, usw. Die Anzahl der Durchgänge werden gezählt und bilden das Ergebnis.

Die Routine wird langsamer als die Funktion aus Version 1 für Quadrahtzahlen ab ca. 27000. Der Einsatz dieser Funktion rentiert sich nur, wenn der Speicherplatz für die Programmcodes knapp wird.

Delphi Routine

Definitionsbereich = 0 ... 2147483647 = MaxInt
Sqrt(2147483647) = 46340,95 => 46340
Function IntSqrt3(w1:dword):integer;
Var k0,i1:dword; k1:word;
Begin
  k0 := 1;
  k1 := 0;
  i1 := w1;
  Repeat
    i1 := i1 - k0;
    If i1 < 0 Then Begin
      Break;
    End Else Begin
      k1 := k1 + 1;
      k0 := k0 + 2;
    End;
  Until i1 <= 0;
  Result := k1;
End;

Assembler Routine für AVR

Definitionsbereich = 0 ... 4294967295
IntSqrt(4294967295) = 65535

Laufzeiten:
SQRT($00000000) = $0000 :         15 Taktzyklen (          0,25µs bei 4MHz)
SQRT($48D0C840) = $8888 : 454401 Taktzyklen (113600,25µs bei 4MHz)
SQRT($FFFFFFFF) = $FFFF : 851988 Taktzyklen (212995,00µs bei 4MHz)

Codelänge: 48 Bytes

;****************************************************************
; Quadratwurzel
;  Input: r16,r17 = Low,High Byte 1
;         r18,r19 = Low,High Byte 2
; Output: r16,r17 = Low,High Byte Quadratwurzel
; Return: -
; Modify: Alle Register > r15
SQRT:      clr   r20          ;Ergbenisspeicher löschen
           clr   r21
           ldi   r22,1        ;Constante 1
           clr   r23
           ldi   r24,2        ;Constante 2
           clr   r25
           clr   r26
           clr   r27
           ldi   r28,1        ;Zwischenspeicher
           clr   r29
           clr   r30
           clr   r31
sqrt01:    sub   r16,r28
           sbc   r17,r29
           sbc   r18,r30
           sbc   r19,r31
           brlo  sqrt02
           add   r20,r22
           adc   r21,r23
           add   r28,r24
           adc   r29,r25
           adc   r30,r26
           adc   r31,r27
           rjmp  sqrt01
sqrt02:    mov   r16,r20
           mov   r17,r21
           ret

5.2.4.3 Version 3

Diese Version arbeitet mit einer konstanten Laufzeit. Sie ist etwas langsamer als die langsamste Laufzeit von Version 1, belegt aber nur halb so viel Programmspeicher. Es wird jedoch zusätzlich eine 32x32 Bit Multiplikation benötigt. Die Codelänge für diese Multiplikation wird jedoch nicht mit eingerechnet, da sie meistens auch von anderen Programmteilen benötigt wird.

Wird sie mit eingerechnet. beträgt die Codelänge 208 Bytes. Sie dadurch länger als die Funktion aus Version 1 und bietet dann nur noch einen Vorteil, wenn der Code optimiert wird.

Delphi Routine

Definitionsbereich = 0 ... 2147483647 = MaxInt
Sqrt(2147483647) = 46340,95
Function IntSqrt(w1:dword):integer;
Var i1:dword; k0,k1:word;
Begin
  k0 := $8000;
  k1 := 0;
  Repeat
    i1 := k1 + k0;
    i1 := i1 * i1;
    If i1 <= w1 Then Begin
      k1 := k1 + k0;
    End;
    k0 := k0 shr 1;
  Until k0 <= 0;
  Result := k1;
End;

Assembler Routine für AVR

Definitionsbereich = 0 ... 4294967295
IntSqrt(4294967295) = 65535

Laufzeiten:
SQRT($00000000) = $0000 : 2338 Taktzyklen (584,50µs bei 4MHz)
SQRT($48D0C840) = $8888 : 2339 Taktzyklen (584,75µs bei 4MHz)
SQRT($FFFFFFFF) = $FFFF : 2338 Taktzyklen (584,50µs bei 4MHz)

Codelänge: 76 Bytes

;****************************************************************
; Quadratwurzel
;  Input: r24,r25 = Low,High Byte 1
;         r26,r27 = Low,High Byte 2
; Output: r24,r25 = Low,High Byte Quadratwurzel
; Return: -
; Modify: Fast alle Register
SQRT:      clr   r12            ;k1 := 0;
           clr   r13
           clr   r14            ;Dummy
           ldi   r28,$80        ;k0 := $8000
           mov   r11,r28
           clr   r10
sqrt01:    mov   r16,r10        ;i1 := k1
           mov   r17,r11
           clr   r18
           clr   r19
           add   r16,r12        ;i1 := i1 + k0
           adc   r17,r13
           adc   r18,r14
           clr   r19
           mov   r20,r16        ;i2 := i1
           mov   r21,r17
           mov   r22,r18
           mov   r23,r19
           rcall umul4          ;i1 := i2 * i1
           cp    r16,r24        ;i1 <= w1 ?
           cpc   r17,r25
           cpc   r18,r26
           cpc   r19,r27
           brlo  sqrt02         ;i1 < w1
           breq  sqrt02         ;i1 = w1
           rjmp  sqrt03
sqrt02:    add   r12,r10        ;k1 := k1 + k0
           adc   r13,r11
sqrt03:    lsr   r11            ;k0 := k0 shr 1
           ror   r10
           tst   r11            ;k0 <= 0 ?
           brne  sqrt01         ;ja
           tst   r10
           brne  sqrt01         ;ja
           mov   r24,r12        ;Result := k1
           mov   r25,r13
           clr   r26
           clr   r27
           ret

Assembler Routine für AVR

Definitionsbereich  Produkt = 0 ... 4294967295
Multiplikant und Multiplikator = 0 ... 65535

Laufzeiten: 119 Taktzyklen (29,75µs bei 4MHz)
Codelänge: 132 Bytes

;****************************************************************
; Multiplikation zweier Vorzeichenlosen 32 Bit Zahlen
;  Input: r16,r17,r18,r19 = Low .. High Multiplikant
;         r20,r21,r22,r23 = Low .. High Multiplikator
; Output: r16,r17,r18,r19 = Low .. High Product Low
;         r20,r21,r22,r23 = Low .. High Product High
; Return: C       = 0 = Das Ergebnis ist 4 Byte
;         C       = 1 = Das Ergebnis ist > 4 Byte
; Modify: -
;                                       r9 r8 r7 r6 r5 r4 r3 r2
; Beispiel: r19r18r17r16*r23r22r21r20 = r23r22r21r20r19r18r17r16
;           10 20 30 40 * 50 60 70 80 =
;           ----------------------------------------------------
; (16)      10          * 50          = 05 00
; (15)      10          *    60       =    06 00 
; (13)      10          *       70    =       07 00   
; (10)      10          *          80 =          08 00
; (14)        20       * 50          =    0A 00
; (12)        20       *    60       =       0C 00
; ( 9)        20       *       70    =          0E 00
; ( 6)         20       *          80 =             10 00
; (11)            30    * 50          =       0F 00
; ( 8)            30    *    60       =          12 00
; ( 5)            30    *       70    =             15 00
; ( 3)            30    *          80 =                18 00
; ( 7)               40 * 50          =          14 00
; ( 4)               40 *    60       =             18 00
; ( 2)               40 *       70    =                1C 00
; ( 1)               40 *          80 =                   20 00
;                                       -----------------------
;                                       05 10 22 3C 3D 34 20 00
umul4:     push   r0
           push   r1
           push   r2
           push   r3
           push   r4
           push   r5
           push   r6
           push   r7
           push   r8
           push   r9

           mul    r16,r20    ;( 1)
           mov    r2,r0
           mov    r3,r1
           mul    r16,r21    ;( 2)
           add    r3,r0
           mov    r4,r1
           mul    r17,r20    ;( 3)
           add    r3,r0
           adc    r4,r1
           mul    r16,r22    ;( 4)
           add    r4,r0
           mov    r5,r1
           mul    r17,r21    ;( 5) 
           add    r4,r0
           adc    r5,r1
           mul    r18,r20    ;( 6)
           add    r4,r0
           adc    r5,r1
           mul    r16,r23    ;( 7)
           add    r5,r0
           mov    r6,r1
           mul    r17,r22    ;( 8)
           add    r5,r0
           adc    r6,r1
           mul    r18,r21    ;( 9)
           add    r5,r0
           adc    r6,r1
           mul    r19,r20    ;(10)
           add    r5,r0
           adc    r6,r1
           mul    r17,r23    ;(11)
           add    r6,r0
           mov    r7,r1
           mul    r18,r22    ;(12)
           add    r6,r0
           adc    r7,r1
           mul    r19,r21    ;(13)
           add    r6,r0
           adc    r7,r1
           mul    r18,r23    ;(14)
           add    r7,r0
           mov    r8,r1
           mul    r19,r22    ;(15)
           add    r7,r0
           adc    r8,r1
           mul    r19,r23    ;(16)
           add    r8,r0
           mov    r9,r1

           mov    r16,r2
           mov    r17,r3
           mov    r18,r4
           mov    r19,r5
           mov    r20,r6
           mov    r21,r7
           mov    r22,r8
           mov    r23,r9

           pop    r9
           pop    r8
           pop    r7
           pop    r6
           pop    r5
           pop    r4
           pop    r3
           pop    r2
           pop    r1
           pop    r0
           ret