[ASM] (schnelle) Integer Wurzel 32bit

[Läubi .. (laeubi)]

Habe leztens mal die Wurzelfunktion aus der Codesammlung ausprobiert und 
war etwas "geschockt", bei großen Werten braucht die doch seeeeeeehr 
lange (über 850.000 clocks! ~ 100ms je nach Operandengröße).
Habe mich dann etwas nach Alternativen umgesehen und eine 
Implementierung gefunden die mit einer konstanten Zeit auskommt und habe 
das ganze mal in ASM umgesezt.

*Ergebnis*: unabhängig vom Operanden etwa 1300 Takte

Im Anhang befindet sich die ASM Datei + eine 
geschwindigkeitsoptimierte Variante welche nur 600 Takte braucht für 
beliebig große 32bit Operanden (verbraucht dann allerdings etwa 
1300bytes an code gegenüber 150byte der nicht optimierten Variante)

Im Source kann ausgewählt werden ob Platz/Speed optimierte Variante 
Verwendung finden soll, und ob das Device den movw Befehl unterstüzt.

Hier mal der Code der platzoptimierten Version, im ASM-File ist das 
vieleicht für den nicht so geübten ASM/AVR Studio Programmierer etwas 
unübersichtlich da ich viel mit den Assemblerdirektiven rumgespielt habe 
:)

wurzel32fast:
/* #######################################################################
Based on an Algorithm by Jim Ulery
  http://www.azillionmonkeys.com/qed/ulerysqroot.pdf

  Ported to ASM by Läubi, 2008

  Input r5:r4:r3:r2
  Output r1:r0

  ! Achtung! Die Eingabe wird hierbei zerstört!
  ! Wenn diese noch benötigt wird vorher sichern!

  C-Version:
  static unsigned julery_isqrt(unsigned long val) {
      unsigned long temp, g=0, b = 0x8000, bshft = 15;
      do {
          if ( val >= ( temp = (((g << 1) + b)<<bshft--) ) ) {
             g += b;
             val -= temp;
          }
      } while (b >>= 1);
      return g;
  }
####################################################################### */
//Register sichern
  push XL
  push XH
  push r17
  push r16
  push r6
  push r7
  push r8
  push r9
  push r10
  //X = b
  ldi XL,  LOW(0x8000)
  ldi XH, HIGH(0x8000)
  //r1:r0 = g = 0
  clr r0
  clr r1
  //temp = bshft
  ldi r17, 15
  clr r16
  wurzel32_2_wloop:
    //temp = g
    clr r9
    clr r8
    mov r7, r1
    mov r6, r0
    //g << 1
    lsl r6
    rol r7
    rol r8
    //leztes rol kann ausgelassen werden da g maximal 16 bit breit
    //und r9 = 0
    //+b
    add r6, XL
    adc r7, XH
    adc r8, r16
    adc r9, r16
    // << bshft
    mov r10, r17
    wurzel32_2_wl:
      lsl r6
      rol r7
      rol r8
      rol r9
      dec r17
    brne wurzel32_2_wl
    //bshft--
    mov r17, r10
    dec r17
    //val >= temp
    cp r2, r6
    cpc r3, r7
    cpc r4, r8
    cpc r5, r9
    brlo wurzel32_2_wd
    //g + = b
    add r0, XL
    adc r1, XH
    //val -=temp
    sub r2, r6
    sbc r3, r7
    sbc r4, r8
    sbc r5, r9
  //b >> 1 != 0?
  wurzel32_2_wd:
    lsr XH
    ror XL
    cpi XL, 1
    cpc XH, r16
    brne wurzel32_2_wloop
  //Register wieder herstellen
  pop r10
  pop r9
  pop r8
  pop r7
  pop r6
  pop r16
  pop r17
  pop XH
  pop XL
  ret

Hinweise sind gerne willkommen, ich möchte aber gleich darauf hinweisen 
das ich NICHT bis auf den lezten Takt optimiert habe (es ist also noch 
Potential drin wenn jemand spaß daran hat das noch weiter zu führen)

[Läubi .. (laeubi)]

Und hier nochmal die "C-Variante" in einer compilierbaren Version welche 
ich zum testen verwendet habe.

[sum (Gast)]

Hi Läubi,

sehr schick. Hast du mal verglichen, was der C-Compiler daraus macht?

Grüesse,
 sum

[Benedikt K. (benedikt)]

Die hier ist schneller und kleiner:
http://elm-chan.org/docs/avrlib/sqrt32.txt
http://elm-chan.org/docs/avrlib/sqrt16.txt

[sum (Gast)]

Habe gerade nochmal nachgeschaut, bei mir macht der Compiler 140 Byte 
draus. Wieviele Takte das braucht, weiss ich natürlich nicht.

sum

[Läubi .. (laeubi)]

sum wrote:
> Hi Läubi,
>
> sehr schick. Hast du mal verglichen, was der C-Compiler daraus macht?
>
> Grüesse,
>  sum
nein, ich hab kein gcc oder ähnliches bei mir installiert. Vieleicht 
kanns ja mal jemand probieren :)

Benedikt K. wrote:
> Die hier ist schneller und kleiner
Kannte ich nicht und hab ich jezt auch nicht getestet, kleiner ja auf 
jedenfall, aber es werden keine Register gerettet wenn du das bei mir 
noch abziehst kommt man auf den selben "speed".
Und ich fands auch mal interesannt sowas selbst zu machen und die Lösung 
aus der Codesammlung scheint mir nicht sehr brauchbar für 32bit (es sei 
den man hat vieeeel Zeit)

[Benedikt K. (benedikt)]

Läubi [e-mail] wrote:

> Benedikt K. wrote:
>> Die hier ist schneller und kleiner
> Kannte ich nicht und hab ich jezt auch nicht getestet, kleiner ja auf
> jedenfall, aber es werden keine Register gerettet wenn du das bei mir
> noch abziehst kommt man auf den selben "speed".

Retten braucht man nichts.
Ich hab die Wurzelroutine von ELM Chan etwas angepasst, so dass man ihn 
direkt in den Compiler einbinden kann. Man muss nur 1 Register sichern. 
Macht 55 Words...

Nutzen lässt sich das ganze dann so:
unsigned int sqrt32(unsigned long);

[Läubi .. (laeubi)]

> var[12:02] = <broken>
hab ich jezt so interpretiert das die Arbeitsregister 12 bis 2 zerstört 
werden.

Gibts eigentlich ne Erklärung WARUM die Routine funktioniert? Stehen ja 
leider keien Kommentare drann oder so :(

[Benedikt K. (benedikt)]

Läubi [e-mail] wrote:
>> var[12:02] = <broken>
> hab ich jezt so interpretiert das die Arbeitsregister 12 bis 2 zerstört
> werden.

Ja, dürfte so sein. Aber das interessiert im Normalfall ja nicht weiter, 
da der Compiler davon ausgeht, dass diese Register überschrieben werden.

> Gibts eigentlich ne Erklärung WARUM die Routine funktioniert? Stehen ja
> leider keien Kommentare drann oder so :(

Ich habe es auch schon aufgegeben die Programme von dem zu verstehen. 
Mir reicht es, wenn ich diese nutze kann.
Ich habe noch die 16bit Wurzelfunktion angehängt:
unsigned char sqrt16(unsigned int);

Die Multiplikatikations und Divisionsroutinen von dem sind auch nett:
Hier die komplette Seite von dem:
http://elm-chan.org/cc_e.html

[Läubi .. (laeubi)]

So, ich hab jezt doch nochmal etwas dran rumgebastelt.

Als erstes hab ich etwas mit der Mathematik drauf eingeschlagen :P

unsigned int julery_isqrt2(unsigned int val) {
  unsigned int temp, g, b, bshft;
  g=0;
  b = 0x8000;
  bshft = 15;
  int i;
  for (i = 0; i < 15; i++) {
    temp = g;
    temp = temp + (b >> 1);
    temp = temp <<(bshft+1);
    if ( val >= temp) {
      g += b;
      val -= temp;
    }
    bshft--;
    b = b >> 1;
  }
  temp = (g<<1);
  temp = temp|1;
  if ( val >= temp) {
    g += b;
    val -= temp;
  }
  return g;
}

Man spart sich hier etwas das rumschieben von g, b und bshft können 
konstant angenommen werden wenn man die schleife ausrollt.
Der C-Compiler sollte es auch etwas einfacher haben das zu optimieren. 
Wenn man eine Architektur mit barrelshifter hat ist der Algorithmus auch 
noch etwas effizienter zu implementieren.

Ergebnis: Schnelle Variante 300-400 Takte, 810bytes code.

Sehr viel Code im vergleich zu gerademal etwa 100 Takten ersparnis im 
Vergleich zu Elm-Chans Version aber fals jemand es besonders eilig hat 
und etwas Codespace übrig ist es vieleicht ne Alternative, und ansonsten 
hat das rumbasteln auch Spaß gemacht :)
Der Code enthält jezt alle Optiemierungen für das schieben die mir 
eingefallen sind (0, 8, 16) der rest wird dann per links/rechts 
nachgeholt. da geht auch die meiste Zeit für drauf.

[Daniel C. (cagara)]

Was macht ihr da eigentlich :) schmunzel

Wurzel in 800 Takten?? uiuiui

Warum nicht so? (x ist anfang der Wert, und am schluss seine wurzel)

int a,b;
b = x;
a = x = 0x3f;
x = b/x;
a = x = (x+a)>>1;
x = b/x;
a = x = (x+a)>>1;
x = b/x;
x = (x+a)>>1;

diese annäherung ist doch ne gute obere schranke für die echte wurzel, 
und sogut wie identisch!

Beschrieben in diesem Paper: 
http://supp.iar.com/FilesPublic/SUPPORT/000419/AN-G-002.pdf

[Andreas S. (andreas) (Admin)]

Und wie lange dauert eine 32 Bit-Division auf dem AVR?

[Löwe (Gast)]

Zu faul zum Lesen - Daniel? Elm Chan braucht max. 185 Takte (16-Bit) und 
das Ergebnis ist genau.

[Gerry M. (gerry501)]

Läubi .. schrieb:
>> var[12:02] = <broken>
> hab ich jezt so interpretiert das die Arbeitsregister 12 bis 2 zerstört
> werden.
>
> Gibts eigentlich ne Erklärung WARUM die Routine funktioniert? Stehen ja
> leider keien Kommentare drann oder so :(

Es wird ein Verfahren zum schriftlichen Wurzelziehen verwendet (ähnlich 
wie schriftlich dividieren):
http://de.wikipedia.org/wiki/Schriftliches_Wurzelziehen

hier ist das ganze für Binär-Zahlen noch verständlicher dargestellt:
http://www.onlinemathe.de/forum/Binaersystem-Zahlentheorie

Elm-Chan hat das Verfahren sehr trickreich umgesetzt. Ein Tip für 
Interessierte (16bit squareroot): Das Bit2 in seinem Register var4 
entspricht bei den 8 Durchläufen dem jeweiligen 'Lösungs-Bit'.

[Johann L. (gjlayde)]

Läubi .. schrieb:

> Im Anhang befindet sich die ASM Datei + eine
> geschwindigkeitsoptimierte Variante welche nur 600 Takte braucht für
> beliebig große 32bit Operanden (verbraucht dann allerdings etwa
> 1300bytes an code gegenüber 150byte der nicht optimierten Variante)
>
> Hinweise sind gerne willkommen, ich möchte aber gleich darauf hinweisen
> das ich NICHT bis auf den lezten Takt optimiert habe (es ist also noch
> Potential drin wenn jemand spaß daran hat das noch weiter zu führen)

Hi, im Wiki gibt's jetzt eine Version die nur noch 315 Ticks braucht und 
zusammen mit avr-gcc läuft:
   http://www.mikrocontroller.net/articles/AVR_Arithmetik#avr-gcc_Implementierung_.2832_Bit.29

[Michael H. (overthere)]

Sieht schon verdammt brauchbar aus!
Könnte noch jemand die 16 Bit variante implementieren und ins Wiki 
schreiben?

[Johann L. (gjlayde)]

Michael H. schrieb:
> Sieht schon verdammt brauchbar aus!
> Könnte noch jemand die 16 Bit variante implementieren und ins Wiki
> schreiben?

Das hat eigentlich alles schon  Ruud v Gessel gemacht. Den Link zu den 
Implementierungen hab ich im Wiki zugefügt. Müsste also nur jemand 
wikitizieren und am besten umstellen auf avr-gcc ABI-konforme 
Registerbelegung. Und noch nen C-Header analog zur 32-Bit Wurzel.

Johann