Division von unsigned char

Man könnte 4 Divisionen a 4 Bit über Tabellen abwickeln. Wird aber auch 
eng.

>Man kann zwar zur Vereinfachung von  Multiplikationen beide Operanden
>in kleinere Einheiten splitten, bei Divisionen geht das aber nur mit
>dem Dividenden.

Man kann dieses Splitting, eg. Karatsuba/Toom-Cook Trick bei 
Multipikationen und auch Divisionen anwenden. Also auch bei der Division 
kann man nach dem "Teile&Herrsche" Prinzip (Divide&Conquer) arbeiten. 
Siehe:

File: TechRep.ps,
"Fast Recursive Division"
Christoph Burnikel & Joachim Ziegler
MPI-I-98-1-022     October 1998
Forschungsbericht     Research Report
Max Planck Institut für Informatik

Allerdings macht dies auf einem AVR keinen Sinn. Selbst Karatsuba für 
Multiplikationen macht auf dem AVR wenig Sinn da die Zahlenbereiche zu 
klein sind und das Taktverhältnis von 1:2 -> ADD:MUL zu gering ist. Hat 
man also eine schnelle HW-MUL dann bestimmt das Verhältnis an Takten die 
man für Additionen/Subtraktionen/Shifts im Vergleich zur 
HW-Multiplikation benötigt über die Sinnhaftigkeit solcher Algortihmen. 
Sie substituieren immer eine math. Operation durch andere meistens 
schnellere Oparationen. Karatsuba spart durch Divide&Conquer Techniken 
eine bestimmte Anzahl an Multiplikationen ein und ersetzt diese durch 
einen erhöhten Aufwand an Additionen/Subtrationen/Shifts und 
Speicherkopierungen. Wenn die Multiplikation aber sehr effizient ist 
dann machen solche Algos. keinen Sinn mehr.

Gruß Hagen

PS: http://cr.yp.to/bib/1998/burnikel.ps

Ich werd mir nachher nochmal 'n Kopf drüber machen, aber im Grunde 
sollte es wie schriftliche Division funktionieren, nur dass jede Ziffer 
einem Nibble entspricht. Man bräuchte dann 2 Tabellen a 256 Einträge.

Aber wie schon oben gesagt, es wird auch damit eng.

Ich reiße mal kurz an, wie ich mir das gedacht hatte. Variablen stehen 
für Ziffern.

Für eine Division AB/CD berechnet man erstmal AB/C0. Der Quotient 
entspricht dem von A0/C0=:Q, der Rest ergibt sich über den Rest von 
A0/C0 + B=:R.

Das Ergebnis hat jetzt erstmal wenig mit AB/CD zu tun. Die Differenz zum 
Quotienten und Rest von AB/CD müsste man aber über 2 LUTs abhängig von Q 
und D ermitteln können (eine für Quotienten, eine für Rest).

Ok, das wären dann 4 Tabellen, nicht 2. Ich bin mir auch noch nicht 
sicher ob das wirklich immer funktioniert, ein paar Testrechnungen haben 
jedenfalls geklappt.
Aber mein Gefühl sagt mir, dass das irgendwie mit einer LUT gehen muss. 
Bewiesen ist natürlich noch garnix.

Nach der Rechnung könnten auch negative Reste auftauchen, und Reste 
größer dem Quotienten. Das muss man dann noch umrechnen.

Vielleicht kann das ja mal jemand von euch nachprüfen. Beispielrechnung 
(dezimal):

56:12 = 4 Rest 8

56:10 = 5 Rest 6 (über 5:1, Rest+6)

Nun braucht man den LUT-Eintrag für Q=5, D=2.

Die LUT kann man füllen indem man einfache irgendeine Testrechnungen 
macht, die auf den Eintrag passt (natürlich 1mal am anfang und nicht 
während der Rechnung).
Im einfachsten Fall bleibt A "ähnlich", also zB. 58:12. 58:10 gibt 5 
Rest 8, passt also. 58:12 macht 4 Rest 10.
Für den LUT-Eintrag ergibt sich damit 0 für Q und -10 für R (R sollte 
immer negativ sein).

Wir waren bei 56:10 = 5 Rest 6. 5 bleibt 5, 6-10 macht -4.
Also ergibt sich für 56:12 = 5 Rest (6-10) = 5 Rest -4 = 4 Rest 8 - 
stimmt offensichtlich.


Komplizierteres Beispiel mit LUT-Eintrag Q=3, D=2:

30:12 = 2 Rest 6, 30:10=3 Rest 0.

LUT-Eintrag mal über was ganz anderes ermittelt: 60:22. Hat mit 60:20=3 
Rest 0 den selben Eintrag, 60:22 = 2 Rest 16. Ergibt sich 0 für Q und -6 
für R.

30:10=3 Rest 0 -> 30:12=3 Rest (0-6) = 3 Rest -6 = 2 Rest 6 - passt 
auch.


Ob man die LUT für Q wirklich braucht weis ich grad nicht, glaub eher 
nicht. Dann wären's noch 3 Tabellen, 2 für 4bit Div und Rest, und eine 
für die Verschiebung vom Rest. Dann sollten die Werte für die 
Restverschiebung auch positiv sein können.

Für Division durch einstellige Werte funktioniert der Algorithmus so 
nicht.

Ok, ganz so einfach geht's doch nicht, bei großen D und kleinen C gibt's 
Probleme. Bei etwa der hälfte aller möglichen Divisionen sollte es 
funktionieren wie oben beschrieben... mal schauen ob mir noch was für 
die andere Hälfte einfällt.
Vielleicht funktioniert's wenn man Q rundet und bei zB. 50:19 nicht 
50:10 sondern 50:20 rechnet, und einen LUT-Eintrag Q=2, D=-1 nachschaut. 
Die Tabelle wird dadurch nicht größer, da D nach wie vor nur 10 Werte 
(in Dez.) annehmen kann. Die Restverschiebungen werden dadurch aber 
kleiner.

Martin wrote:
> Hallo,
>
> kann mit jemand einen schnellen Algorithmus für eine Division zweier
> unsigned char Variablen nennen? Als Controller habe ich den ATMEGA8. Es
> stehen nur 16 Takte zur Verfügung.
>
> Geht das überhaupt? Meine erste Idee war eine Tabelle zu verwenden.
> Diese wird aber viel zu groß.
>
> Martin

Die Sterne stehen schlecht...

Mit folgendem Verfahren (16-Bit Reziprok-Tabelle) braucht's trotz der 
Tabelle über 20 Ticks. Dafür stimmen die Ergebnisse exakt mit a/b 
überein (unsigned).

1

static unsigned short quo[256] =

2

{

3

    0, 0, 384, 22101, 320, 13363, 43818, 37668, 288, 29212,

4

    39449, 17943, 22037, 45587, 18962, 4625, 272, 4111, 14606, 31245,

5

    52492, 12556, 41739, 8715, 43786, 15882, 55561, 31753, 9481,

6

    54280, 35080, 17160, 264, 49671, 34823, 20743, 7431, 60422, 48390,

7

    37126, 26374, 16134, 6406, 62725, 53765, 45317, 37125, 29445,

8

    22021, 14853, 7941, 1541, 60676, 54532, 48644, 43012, 37636,

9

    32260, 27140, 22276, 17668, 13060, 8708, 4356, 260, 61699, 57603,

10

    54019, 50179, 46595, 43267, 39939, 36611, 33283, 30211, 27139,

11

    24323, 21507, 18691, 15875, 13315, 10755, 8195, 5635, 3331, 1027,

12

    64258, 61954, 59650, 57602, 55554, 53506, 51458, 49410, 47618,

13

    45570, 43778, 41986, 40194, 38402, 36866, 35074, 33538, 32002,

14

    30466, 28930, 27394, 25858, 24322, 23042, 21506, 20226, 18946,

15

    17410, 16130, 14850, 13570, 12546, 11266, 9986, 8962, 7682, 6658,

16

    5378, 4354, 3330, 2306, 1282, 258, 64769, 63745, 62721, 61697,

17

    60673, 59905, 58881, 57857, 57089, 56065, 55297, 54529, 53505,

18

    52737, 51969, 51201, 50177, 49409, 48641, 47873, 47105, 46337,

19

    45825, 45057, 44289, 43521, 42753, 42241, 41473, 40705, 40193,

20

    39425, 38913, 38145, 37633, 36865, 36353, 35585, 35073, 34561,

21

    33793, 33281, 32769, 32257, 31489, 30977, 30465, 29953, 29441,

22

    28929, 28417, 27905, 27393, 26881, 26369, 25857, 25345, 24833,

23

    24321, 23809, 23297, 22785, 22529, 22017, 21505, 20993, 20737,

24

    20225, 19713, 19201, 18945, 18433, 18177, 17665, 17153, 16897,

25

    16385, 16129, 15617, 15361, 14849, 14593, 14081, 13825, 13313,

26

    13057, 12545, 12289, 12033, 11521, 11265, 10753, 10497, 10241,

27

    9729, 9473, 9217, 8705, 8449, 8193, 7937, 7425, 7169, 6913, 6657,

28

    6401, 5889, 5633, 5377, 5121, 4865, 4609, 4097, 3841, 3585, 3329,

29

    3073, 2817, 2561, 2305, 2049, 1793, 1537, 1281, 1025, 769, 513

30

};

31

32

#define DIV8_TAB(__a,__b)                       \

33

({                                              \

34

    unsigned char __a2 = __a;                   \

35

    unsigned char __b2 = __b;                   \

36

    unsigned short * __p1b = & quo[__b2];       \

37

    unsigned char __ab;                         \

38

    asm ("lpm  r0, z+"             "\n\t"       \

39

         "mul  %[a], r0"           "\n\t"       \

40

         "mov  %[c], r1"           "\n\t"       \

41

         "mov  r18, r0"            "\n\t"       \

42

         "lpm  r0, z"              "\n\t"       \

43

         "mul  %[a], r0"           "\n\t"       \

44

         "add  r18, r1"            "\n\t"       \

45

         "clr  __zero_reg__"       "\n\t"       \

46

         "adc  %[c], __zero_reg__"              \

47

         : [c] "=&r" (__ab), "+z" (__p1b)       \

48

         : [a] "r"   (__a2)                     \

49

         : "18");                               \

50

    __ab;                                       \

51

})

52

53

unsigned char div8_tab (unsigned char a, unsigned char b)

54

{

55

    if (b != 1)

56

        a = DIV8_TAB (a,b);

57

58

    return a;

59

}

Getestet hab ich's aber nur für i386 in C und für AVR mit Brain 0.9b 
simuliert.

Der erzeugte Assembler sieht so aus:

1

.global  div8_tab

2

  .type  div8_tab, @function

3

div8_tab:

4

/* prologue: frame size=0 */

5

/* prologue end (size=0) */

6

  cpi r22,lo8(1)   ;  b,   ;  12  cmpqi/2  [length = 1]

7

  breq .L2   ; ,   ;  13  branch  [length = 1]

8

  mov r30,r22   ;  __p1b, b   ;  21  zero_extendqihi2/2  [length = 2]

9

  clr r31   ;  __p1b

10

  add r30,r30   ;  __p1b, __p1b   ;  23  *addhi3/1  [length = 2]

11

  adc r31,r31   ;  __p1b, __p1b

12

  subi r30,lo8(-(quo))   ;  __p1b,   ;  24  *addhi3/4  [length = 2]

13

  sbci r31,hi8(-(quo))   ;  __p1b,

14

/* #APP */

15

  lpm  r0, z+

16

  mul  r24, r0   ;  a

17

  mov  r25, r1   ; 

18

  mov  r18, 0

19

  lpm  r0, z

20

  mul  r24, r0   ;  a

21

  add  r18, r0

22

  clr  __zero_reg__

23

  adc  r25, __zero_reg__   ; 

24

/* #NOAPP */

25

  mov r24,r25   ;  a,   ;  52  *movqi/1  [length = 1]

26

.L2:

27

  clr r25   ;  <result>   ;  40  zero_extendqihi2/1  [length = 1]

28

/* epilogue: frame size=0 */

29

  ret

viel unnötig ist da nicht...

Vielleicht kannst Du die Echtzeit-Anforderung Deiner Anwendung 
entschärfen. Ausserdem steht bei Frequenzmessung im Zähler womöglich ne 
Konstante?

Ich hab meine Variante mal implementiert und alle möglichen 
uchar-divisionen durchgetestet. So wie's oben dasteht stimmen 58737 von 
65536 Rechnungen (div durch 0 als stimmt immer angenommen, div durch 
einstelliges immer als falsch (betrifft 2304 Rechnungen)), Quotient 
und Rest.

Schaut also schonmal ganz gut aus, das sind 89.6%. Irgendwie wird man 
den Rest auch noch an's Laufen bekommen.

Während der Rechnung lässt sich übrigens absehen, ob das Resultat gültig 
sein wird. Da könnte man abbrechen und eine normale Division 
durchführen. Bringt dann halt nur in 90% aller Fälle deutlich mehr 
Geschwindigkeit. Aber in 16 Takte wird es, wie gesagt, nicht passen.

Die einstelligen Divisionen betreffen 3840, nicht 2304 Rechnungen.

Da ich kein AVR Assembler kann, hier mal der Code in x86 asm. Ihr könnt 
ja dann schauen ob das in 16 Takte passt:

Gerechnet wird ax/bx, edi zeigt auf 4bit div lut, esi auf 4bit mod lut 
(jeweils a/b -> index=a*16+b), und ebp auf die mit der Restverschiebung.
Auf'm AVR gehen einige der Operationen sicher mit weniger Befehlen als 
auf x86. Vielleicht ließen sich die Tabellen auch besser anordnen.

1

MOV  CL, AL    // AB/CD 

2

AND  CL, 0xF0

3

4

MOV  CH, BL

5

SHR  CH, 4

6

AND  CL, CH    // CL=AC

7

8

MOV  DL, [EDI+CL]  // div - &0xF0 is always 0

9

MOV  DH, [ESI+CL]  // mod

10

      // DL=A0/C0

11

      // DH=A0%C0=:Q1

12

AND  AL, 0xF    

13

ADD  DH, AL    // DH=AB%C0

14

15

16

MOV  CL, DL    // EBP tabelle umsortieren

17

SHL  CL, 4    // spart instruktionen

18

MOV  CH, BL

19

AND  CH, 0x0F

20

AND  CL, CH    // CL=Q1D

21

22

ADD  DH, [EBP+CL]  // restverschiebung

23

24

JNBZ  noadjust

25

26

DEC  DL

27

28

noadjust:

So in etwa müsste das aussehen, für die Korrektheit übernehm ich keine 
Garantie, hab das nicht ausprobiert (es wird nur der Quotient berechnet, 
damit der Rest stimmt müssen noch 1 oder 2 Anweisungen dazu). Das 
Displacement mit CL funktioniert glaub ich so nicht (hab auf x86 schon 
ewig nimmer mit 8 und 16bit Registern rumgemacht), ist aber egal, beim 
AVR hat man das Problem nicht.

Es sind genau 16 Befehle, es fehlt aber die Erkennung ob die Operation 
überhaupt möglich ist. Es gibt 2 Punkte wo's scheitern kann. Wenn 
BL&0xF0==0 (könnte man in der DIV oder MOD Tabelle mit einem extra Wert, 
zB. 255 markieren), und wenn die Restverschiebung zu groß ist (in der 
Restverschiebungstabelle markieren, zB. auch 255).

Vielleicht könnte man das auch zusammenfassen, indem man zB. in die DIV 
Tabelle einen Wert schreibt, der dann in der Restverschiebungstabelle 
einen ungültigen Wert gibt. Dann reicht eine Abfrage.

Wenn man den Rest mitberechnet kann man die Erfolgsquote mit einem 2. 
Adjust Durchlauf (wenn Rest <0 Rest+Divisor, Quotient--) auf etwa 95% 
steigern.

Wenn man weiß in welchen Fällen der Algorithmus fehlschlägt, ist doch 
alles in Ordnung.

Simon K. wrote:
> Wenn man weiß in welchen Fällen der Algorithmus fehlschlägt, ist doch
> alles in Ordnung.

Eben, und dann kann man die langwierige Division durchführen. Wobei sich 
für einstellige Divisoren sicher auch schnellere Varianten finden 
lassen.

Wenn's um einen Algorithmus geht wo viele Divisionen mit 
gleichverteilten Parametern vorkommen, und der Algorithmus von oben 1/4 
der Zeit einer normalen Division braucht, reduziert sich die benötigte 
Zeit auf 30%.

Bei CPUs mit Divisionsbefehl dauert diese übrigens auch nicht immer 
gleich lang.

Mit einer zusätzlichen 256 Byte Tabelle müssten auch einstellige 
Divisionen funktionieren. So ungefähr, am Beispiel 51/4:

51/4 = 50/4+1/4

1/4 kann man über die vorhandene Tabelle berechnen, 0 Rest 1

50/4 geht so: Über die vorhandene Tabelle 5/4, 1 Rest 1. Der Rest 
braucht nie mehr als 4 Bit (bzw. eine Stelle Dez.), da Divisior immer 
einstellig.

Der Rest und der Quotient von 1/4 beschreiben nun die letzte Stelle vom 
Quotienten von 51/4, der Quotient aus 5/4 die erste Stelle.

Man müsste die letzte Stelle über eine Tabelle, die den Divisor und den 
Rest von 5/4 auf die letzte Stelle mappt, berechnen können - dazu dann 
noch der Quotient von 1/4.

In dem Fall hier wäre das dann Q=4, R=1 -> 2.

Andere Rechnung mit dem Wert: 59/4. 5/4 ist wieder 1R1, 9/4=2R1. 59/4 
ist dann also (1*10+2)+2=14 - stimmt. Der Rest stimmt nicht, um den zu 
berechnen bräuchte man noch eine Tabelle.


Auf die Weise stimmen schonmal 63453 von 65536 möglichen Rechnungen (mit 
2fachem Test statt 1fachem wie im Code oben), also 96.8%. Insg. braucht 
man dafür 1kB Tabellen.

Wenn man so ca. 4kB verheizt bekommt man übrigens noch schnellere 
Algorithmen. Damit könnte man dann vll sogar unter 16 Takte kommen. Aber 
4kB ist schon einiges.

Ich schrieb "nicht immer", nicht "stets". Auf den MIPS Kollegen dauert 
das zB. unterschiedlich lang. Aber ich nehm an du bist der Typ aus dem 
Fusionsthread - erledigt.


@crazy horse

Ist ja nicht ausgeschlossen, dass sich noch was schnelles findet. Wenn 
sich hier eine gute Lösung findet, kann man die ja auch in anderen 
Problemstellungen verwenden. Ich find das Thema ansich auch nicht 
uninteressant.

Man könnte es auch so probieren:

Herkömmlich wäre zB. so (a/b): b shift left bis oberes Bit gesetzt, 
Schiebeweite merken (=:l). Ergebnis e:=0.

b>=a? wenn ja: a-=b, e|=0x80
b>>=1;
b==0 -> springe zu ende

b>=a? wenn ja: a-=b, e|=0x40
b>>=1;
b==0 -> springe zu ende

usw, die Maske verschieb sich jeweils um 1 (man könnte auch immer das 
gleiche Bit in e setzen, und dann schieben - ist aber langsamer).

ende:

e>>=l;


Im dümmsten Fall braucht man 7 Schritte (a/1), im einfachsten Fall nur 
einen (b>=128) und man ist sehr schnell fertig.
Für a<16 braucht man 4..7 Schritte, die Divisionen durch einstellige 
Divisoren braucht also am längsten.

Man nimmt sich nun eine Tabelle her, die für alle 256 Divisoren die 
Schrittzahl enthält (256 Byte) - damit kann man auch auf l schließen und 
muss das nicht berechnen.

Bei <= 4 Schritten rechnet man wie oben durch und hat für nachgucken in 
der Tabelle, springen, und 4 Schritte 16 Takte. Könnte reichen, bleibt 
aber eng.

Bei > 4 Schritten benutzt man einen anderen Algorithmus. Man könnte eine 
fertige Tabelle für 8bit/4bit benutzen, die braucht 4kB. Damit sollte es 
für jede Rechnung in <=16 Takten gehen.


Wenn man's schafft die Werte einzuschränken ist allerdings wirklich sehr 
viel gewonnen. Je mehr man optimiert, desto weiter muss man sich von 
allgemeingültigen Algorithmen entfernen.

crazy horse wrote:
> ich halte nix von solchen Ausquetschereien. Was bringt ein Algorithmus,
> der in 90% der Fälle die Arbeit in der geforderten Zeit erledig, in den
> restlichen 10% aber versagt?

Nein, das muss man völlig anders sehen. Möglicherweise lassen sich die 
Spezialfälle, in denen der Algorithmus nicht funktioniert ganz einfach 
per IF-Abfrage abfangen und sofort ein festes Ergebnis zurückgeben. Dann 
wäre die Laufzeit des Algorithmus unter Umständen sogar kürzer!

Außerdem kann man diesen "halb funktionierenden" Algorithmus ja auch 
einstellen, damit andere mit der Idee weiterarbeiten können und sich der 
Algorithmus mglw. noch verbessert.

Ich erinnere da nur an die Algorithmen, die es alle schon gab um 
Primzahlen hervorzusagen, welche allerdings nur bis zu einer bestimmten 
Zahl ging und ab da sehr häufig Fehler aufwiesen.

Siggi wrote:
> @Johann
>
> Dein Algorithmus funktioniert nicht. 200 geteilt durch 7 ergibt 4!?
>
> Getestet mit GCC und AVR-Studio 4.


hmmm...

An Stelle 7 der Tabelle steht 37668, also LOW=36, HIGH=141.

Zunächst wird das LOW-Byte mit 200 multipliziert, gibt 200*36 = 28*256 + 
32.

Nach dem ersten MUL müsste R1=%[c] also 28=0x1c sein und R0=R18=32=0x20.

Dann wird gerechnet 200*141=110*256 + ...

110 + R18 gibt kein 8-Bit-Überlauf, daher bleibt das Ergebnis auf 28, 
dem korrekten Ergebnis.

Von den Werten her stimmt es, es müsste ein blöder Fehler im asm sein, 
den ich nicht sehe (hab leider keine Testplattform ausser Brain 1.0)

Im asm-Template fehlen ein paar Unterstriche, die lässt die 
Foren-Software gerne großzügig unter den Tisch fallen. Hast Du die 
nachgebessert? Ich weiß nicht, wie ich's hier hinschreiben soll, weil 
die Unterstriche ja flöten gehen... So wie ~~zero~reg~~ muss es heissen.

Hier ne bessere Version, die nicht explizit R18 belegt:

#define DIV8_TAB(__a,__b)                       \
({                                              \
    unsigned char __a2 = __a;                   \
    unsigned char __b2 = __b;                   \
    unsigned short * __p1b = & quo[__b2];       \
    unsigned char __ab;                         \
    unsigned char __d;                          \
    asm ("lpm  r0, z+"             "\n\t"       \
         "mul  %[a], r0"           "\n\t"       \
         "mov  %[c], r1"           "\n\t"       \
         "mov  %[d], r0"           "\n\t"       \
         "lpm  r0, z"              "\n\t"       \
         "mul  %[a], r0"           "\n\t"       \
         "add  %[d], r1"           "\n\t"       \
         "clr  __zero_reg__"       "\n\t"       \
         "adc  %[c], __zero_reg__"              \
         : [c] "=&r" (__ab), "+z" (__p1b), [d] "=&r" (__d) \
         : [a] "r"   (__a2));                   \
    __ab;                                       \
})