8 Bit Berechnung schnell, schneller am schnellsten

Das lässt sich doch bestimmt verallgemeinern.

Vielleicht solltest Du verraten, um was für einen 8 Bit-Computer es 
geht?
Ich vermute, es gibt bessere Orte, um nachzufragen...

Und allenfalls sollte man sich ueberlegen, ob double wirklich noetig 
ist. Worum geht es denn ?

Idee: Addition von Logarithmen statt Multiplikation.

Es gilt:

und somit

Wenn man das Logarithmieren und die e Funktion durch Tabellen ersetzt 
bleibt eine Addition. Evtl. muss man aus Genauigkeitsgruenden mit 16 bit 
Zwischenwerten rechnen.

Andere Idee:

Vermutlich sind ja die wx, wz und wy fuer viele Punkte konstant (wie 
immer fehlen die Details hierzu). Man koennte fuer die konstanten 
Multiplikationen Code erzeugen, der die Multiplikation mittels 
Additionen macht.
Kann man Code erzeugen, oder muss aller Code ins FLASH, da es nicht 
genug RAM gibt ? Wir wissen es nicht, denn wir wissen nicht fuer welchen 
Prozessor !

ZigZeg

Weitere Ideen:
Naiv ist hat eine Multiplikationstabelle 256*256=64k Eintraege, also zu 
viel. Meist handelt es sich um vorzeichenbehaftete Werte. Wenn man die 
Vorzeichen separat behandelt, so braucht man nur noch 128*128=16k 
Eintraege. Wenn noetig kann man noch z.B. bei einem Multiplikanden ein 
Bit Genauigkeit einsparen, und man hat nur noch 8k Eintraege.

Auch:
Wird der Z-Wert wirklich benoetigt ? Wenn nicht spart man immerhin 1/3 
des Berechnungsaufwandes.

Und:
Wenn es Symmetrien in der Punktwolke gibt, so reicht es oft z.B. ein 
Viertel der Punkte zu transformieren. Erst wenn man die Punkte 
transformiert hat "spielt man ein bischen an den Vorzeichen", und kommt 
(nach addition eines Mittelpunkts) auf die volle Punktwolke.

Man benötigt:

6 add/sub
6 lookups
6 mul

Lookup braucht min. 4 takte, die adda einen. Bleiben noch knapp 30 je 
mul.

Sportlich, aber nicht unmöglich. Setz mich heut Abend mal dran.

E: Es geht um einen tiny, richtig?

Eine 16x16 Multiplikation benoetigt etwa 16 shifts und 16 additionen bei 
16bit registern. Bei 8 bit registern ist es etwas mehr.

Um nur Punkte auf dem Kreis zu haben, wuerde ich mir N Punkte als 
Koordinaten auf einem 45 Grad Stueck merken, der Rest des Kreises ist 
triviale Spiegelung. Zwischen diesen Punkten linear interpolieren. 
Welche Genauigkeit der Punkte ist denn noetig ?

Hier reichen aber 8*8 bit mit 8 bit Ergebnis.

Also 8 add 8 shift und 8 vergleiche plus evtl etwas  nacharbeit.

Noch schneller geht es mit folgender algebraischen Identität:

(a+b)² - (a-b)² = a^2 + 2ab + b² - a² + 2ab - b² = 4 ab

Man benötigt nur recht kleine Tabellen für die Quadratzahlen.

Auf einem C64 kommt man aber auch damit nicht sehr weit...

Falls du die Multiplikation wirklich per Tabelle abwickeln willst, dann 
so:

Um a*b zu berechnen schreibt man a=u+v und b=u-v, somit ist ab = u²-v²; 
man braucht also nur Quadrate abzuspreichern, was Platz O(n) kostet 
anstatt O(n²) wie oben vorgeschlagen.  Die Abbildung von a, b auf u, v 
ist eine simple lineare Transformation.

c-hater schrieb:
> So schön der Ansatz nämlich theoretisch auch aussieht, kostet die
> praktische Umsetzung doch mehr Takte als eine echte
> Multiplikationsroutine. Und die braucht überhaupt keine Tabellen...

Soso... Zahlen, Daten, Fakten bitte. Zielprozessor scheint übrigens ein 
6510 zu sein.

Peter, was ist es denn jetzt für ein Controller und was für ein 
Datenformat brauchst du bei der AddPixel Funktion? 8Bit oder 16Bit?

8,16-Bit, integer, festkomma?

c-hater schrieb:
> Tim    schrieb:
>
>> Soso... Zahlen, Daten, Fakten bitte.
>
> AVR 8x8 ohne MUL-Instruktion: 40 Takte. Das darfst du gern mit deinem
> Ansatz unterbieten, ich wäre keinesfalls sauer, wenn du es könntest.
> Multiplikation auf'm Tiny ist ein ewiger Krampf,da ist jede effizientere
> Umsetzung hochwillkommen...

1

uint16_t Quadrate[256];  //Enthält Quadratzahlen geteilt durch 4

2

3

uint8_t a,b;

4

uint16_t result;

5

6

result=Quadrate[a+b]-Quadrate[a-b];

Kannst Du Dir jetzt selber in Assembler umsetzen. Aber selbst mit einem 
C-Compiler werden da wohl kaum 40 Takte herauskommen.

Peter: Du denkst noch nicht in 6510 Assembler... Aber Du kannst Dir das 
natürlich auch mit Variablen in der Zeropage realisieren. Sollte 
funktionieren, wird aber nicht schnell.

Tim    schrieb:
> c-hater schrieb:
>> Tim    schrieb:
>>
>>> Soso... Zahlen, Daten, Fakten bitte.
>>
>> AVR 8x8 ohne MUL-Instruktion: 40 Takte. Das darfst du gern mit deinem
>> Ansatz unterbieten, ich wäre keinesfalls sauer, wenn du es könntest.
>> Multiplikation auf'm Tiny ist ein ewiger Krampf,da ist jede effizientere
>> Umsetzung hochwillkommen...
> uint16_t Quadrate[256];  //Enthält Quadratzahlen geteilt durch 4
>
> uint8_t a,b;
> uint16_t result;
>
> result=Quadrate[a+b]-Quadrate[a-b];
>
> Kannst Du Dir jetzt selber in Assembler umsetzen. Aber selbst mit einem
> C-Compiler werden da wohl kaum 40 Takte herauskommen.

Ich sollte ergänzen, dass das in der Form erst einmal eine 7:7->15 
Multiplikation ergibt. Außerdem kann es Ungenauigkeiten bei den LSB 
geben.

Dafür muss man sich, je nach Anwendung, eine Sonderbehandlung einbauen. 
Und/oder zusätzlich shiften.

c-hater schrieb:
> Tim    schrieb:
>
>> Soso... Zahlen, Daten, Fakten bitte.
>
> AVR 8x8 ohne MUL-Instruktion: 40 Takte. Das darfst du gern mit deinem
> Ansatz unterbieten, ich wäre keinesfalls sauer, wenn du es könntest.

Bitteschön:

AVR hier zwar OT, aber hab's trotzdem mal ausgetextet:

1

#define A0  24

2

#define A1  A0+1

3

4

#define B0  22

5

#define B1  B0+1

6

7

#define Z0  30

8

#define Z1  Z0+1

9

10

#define AB0  24

11

#define AB1  AB0+1

12

.text

13

14

;; unsigned mul_quad (unsigned a, unsigned b)

15

;; return a * b

16

17

mul_quad:

18

  movw  Z0, A0

19

  add  Z0, B0

20

  adc  Z1, B1

21

  bst  Z0,0

22

  andi  Z0, ~1

23

  subi  Z0, lo8 (-(lupo))

24

  sbci  Z1, hi8 (-(lupo))

25

  lpm  AB0, Z+

26

  lpm  AB1, Z

27

28

  sbc  Z0, B0  ;; sic!

29

  sbc  Z1, B1

30

  sub  Z0, B0

31

  sbc  Z1, B1

32

  lpm  0, Z+

33

  sub  AB0, 0

34

  lpm  0, Z

35

  sbc  AB1, 0

36

37

  brtc  0f

38

  add  AB0, B0

39

  adc  AB1, B1

40

41

0:  ret

42

43

.global  mul_quad

44

.type  mul_quad, @function

45

.size  mul_quad, .-mul_quad

46

47

.macro  quads from=0, to

48

  .word \from * \from

49

  .if   \to-\from

50

    quads "(\from+1)",\to

51

  .endif

52

.endm

53

54

.section .progmem.data.lupo, "a", @progbits

55

56

lupo:

57

  quads 0, 99

58

  quads 100, 199

59

  quads 200, 255

60

61

.type  lupo, @object

62

.size  lupo, .-lupo

> Multiplikation auf'm Tiny ist ein ewiger Krampf,da ist jede effizientere
> Umsetzung hochwillkommen...

Obige Routine braucht weniger als 30 Ticks (ohne RET und CALL). 
Übrigens ist es keine 7*7 Routine sondern sie kann 8*8 = 16, allerdings 
nur unsigned.

Zugegeben, das hat 2 Schönheitsfehler:

1) Die riesige Tabelle (512 Bytes).

2) Es muss a >= b sein, d.h. 1. Argument >= 2. Argument

Mit Behebung von 2) dürfte das allerdings auch in den Bereich von 40 
Ticks kommen.

c-hater schrieb:
>> uint16_t Quadrate[256];  //Enthält Quadratzahlen geteilt durch 4
^^^^ hier verbirgt sich ein wichtiger hinweis  ^^^^^^^^

>> C-Compiler werden da wohl kaum 40 Takte herauskommen.
>
> a=1
> b=2

((1+2)²)/4-((1-2)^2)/4= 9/4 - 1/4

Mit Nachkommastellen: 9/4-1/4=2
Nur Integer: 9/4-1/4=2-0=2

> a=255
> b=2

Hier muss man noch über Vorzeichen diskutieren.

Johann L. schrieb:
> 2) Es muss a >= b sein, d.h. 1. Argument >= 2. Argument
>
> Mit Behebung von 2) dürfte das allerdings auch in den Bereich von 40
> Ticks kommen.

Ach wat. Das ist ein Vergleich plus Branch, also 3 Takte mehr. Es müssen 
ja nur beide Faktoren vertauscht werden.

Wofür eigentlich die 16 Bit Eingangswerte?

Du kannst außerdem noch Code sparen, in dem Du die lookup-Table auf eine 
Pagegrenze legst.

Tim    schrieb:
> Johann L. schrieb:
>> 2) Es muss a >= b sein, d.h. 1. Argument >= 2. Argument
>>
>> Mit Behebung von 2) dürfte das allerdings auch in den Bereich von 40
>> Ticks kommen.
>
> Ach wat. Das ist ein Vergleich plus Branch, also 3 Takte mehr. Es müssen
> ja nur beide Faktoren vertauscht werden.

Ja, stimmt.  Aus den 28 Ticks werden dann 32.

> Wofür eigentlich die 16 Bit Eingangswerte?

Weil man um diese Uhrzeit keinen Code mehr schreiben sollte ;-)

Allerdings kostet es in der Variante nicht mehr, d.h. anstatt einer 
8*8=16 bekommt man eine 16*16=16 Multiplikation für lau.

Mit der Version, die keine Sortierung der Argumente voraussetzt, kostet 
die 16*16=16 dann allerdings 2 Ticks mehr als die 8*8=16.

So geht es auf dem 6502:

1

; Description: Unsigned 8-bit multiplication with unsigned 16-bit result.

2

;                                                                        

3

; Input: 8-bit unsigned value in T1                                      

4

;        8-bit unsigned value in T2                                      

5

;        Carry=0: Re-use T1 from previous multiplication (faster)        

6

;        Carry=1: Set T1 (slower)                                        

7

;                                                                        

8

; Output: 16-bit unsigned value in PRODUCT                               

9

;                                                                        

10

; Clobbered: PRODUCT, X, A, C                                            

11

;                                                                        

12

; Allocation setup: T1,T2 and PRODUCT preferably on Zero-page.           

13

;                   square1_lo, square1_hi, square2_lo, square2_hi must be

14

;                   page aligned. Each table are 512 bytes. Total 2kb.    

15

;                                                                         

16

; Table generation: I:0..511                                              

17

;                   square1_lo = <((I*I)/4)                               

18

;                   square1_hi = >((I*I)/4)                               

19

;                   square2_lo = <(((I-255)*(I-255))/4)                   

20

;                   square2_hi = >(((I-255)*(I-255))/4)                   

21

.proc multiply_8bit_unsigned                                              

22

                bcc :+                                                    

23

                    lda T1                                                

24

                    sta sm1+1                                             

25

                    sta sm3+1                                             

26

                    eor #$ff                                              

27

                    sta sm2+1                                             

28

                    sta sm4+1                                             

29

                :                                                         

30

31

                ldx T2

32

                sec   

33

sm1:            lda square1_lo,x

34

sm2:            sbc square2_lo,x

35

                sta PRODUCT+0   

36

sm3:            lda square1_hi,x

37

sm4:            sbc square2_hi,x

38

                sta PRODUCT+1   

39

40

                rts

41

.endproc

Von hier: http://codebase64.org/doku.php?id=base:6502_6510_maths

(Google! Hätte ich auch vorher dran denken können...)

Der Trick an dem Code oben ist, dass ein Faktor per 
selbstmodifizierendem Code geschrieben wird.

Bei einer Vektor/Matrix Rotation kann das z.B. die Rotationsmatrix sein, 
die sich nur selten verändert - und plötzlich landet man bei 4-6 
Befehlen pro Multiplikation. Die Addition in der Matrix gibt es dann 
noch kostenlos dazu.

Hier ist eine komplette Implementation einer 3D Rotation mit 
perspektivischer Projektion in 6502 Assembler:

http://codebase64.org/doku.php?id=base:3d_rotation

AVR-Assembler ist für Anfänger, kann ich dazu nur sagen...

Johann L. schrieb:
> Tim    schrieb:
>> Johann L. schrieb:
>>> 2) Es muss a >= b sein, d.h. 1. Argument >= 2. Argument
>>>
>>> Mit Behebung von 2) dürfte das allerdings auch in den Bereich von 40
>>> Ticks kommen.
>>
>> Ach wat. Das ist ein Vergleich plus Branch, also 3 Takte mehr. Es müssen
>> ja nur beide Faktoren vertauscht werden.
>
> Ja, stimmt.  Aus den 28 Ticks werden dann 32.

Hier noch die AVR-Version mit den 32 Ticks und ohne Restriktionen an die 
Parameter, wie oben gemäß avr-gcc ABI aber etwas mehr kommentiert.

Die Berechnung verwendet

Wie gesagt brauch das 32 Ticks; da scheint irgendwie die Schallgrenze zu 
sein.  Tabelle auf 2^8 Bytes alignen ist m.E. Overkill, aber wer vor 
einer Tabelle zur Multiplikation nicht zurückschreckt, den stört das 
"bisschen" Mehrverbrauch um einen einzigen Tick zu sparen wohl auch 
nicht.

Die Routine lässt sich auch auf eine 16*16=16 Multiplikation erweitern, 
und falls der Tiny ein MOVW hat dann geht das ebenfalls mit 32 Ticks. 
Allerdings wächst die Tabelle exponentiell mit der gewünschten Bitbreite 
der Eingabe, und das wird dann rasch einfach nur noch Aua...

Man mag einwenden, dass die verwendete Formel zu kompliziert sei und nur 
einem Theoretiker einfallen könne, aber immerhin läuft das 20% schneller 
als die oben gesetzten 40 Ticks — und zwar in der Praxis :-)

Einige der o.g. Verfahren verwenden (a+-b)^2 was dann eine Tabelle mit 
den Quadraten (bzw. 1/4 der Quadrate) von 0..511 erfordert.  Der 
folgende Ansatz braucht nur die Quadrate von 0..255, und mit einer 
größeren Tabelle wird's auch nicht schneller.  Zumindest sehr ich da auf 
Anhieb nix mehr — zumindest nicht für AVR.

1

zero = 1

2

tmp  = 0

3

4

.text

5

6

#define A0  24

7

#define B0  22

8

9

#define Z0  30

10

#define Z1  Z0+1

11

12

#define AB0  24

13

#define AB1  AB0+1

14

15

;; uint16_t mul_quad (uint8_t a, uint8_t b)

16

;; Return: a * b

17

;; Passend zur  avr-gcc ABI

18

19

mul_quad:

20

    ;; Z = zero-extend (A) + zero-extend (B)

21

    mov     Z0, A0

22

    add     Z0, B0

23

    clr     Z1

24

    adc     Z1, zero

25

    ;; B = min (A, B)

26

    cp      A0, B0

27

    brcc    1f

28

    mov     B0, A0

29

1:

30

    ;; T = (A + B) mod 2

31

    bst     Z0, 0

32

    ;; Z -= T  -->  Z mod 2 = 0

33

    andi    Z0, ~1

34

    ;; lupo[n] = n^2; n = 0 ... 255

35

    ;; Z = lupo + (A + B) div 2

36

    subi    Z0, lo8 (-(lupo))

37

    sbci    Z1, hi8 (-(lupo))

38

    ;; AB = lupo[(A + B) div 2] = (A + B - T)^2 / 4

39

    lpm     AB0, Z+

40

    lpm     AB1, Z

41

42

    ;; x + y - 2 min (x, y) = |x - y|  mit  x, y >= 0

43

    ;; --> Z = lupo + |A - B - T| / 2

44

    sbc     Z0, B0  ;; sic!  wegen C = 1 wird das Z+ von oben aufgehoben

45

    sbci    Z1, 0

46

    sub     Z0, B0

47

    sbci    Z1, 0

48

49

    ;; AB -= lupo[|A - B - T| / 2]

50

    ;; AB = ((A + B - T)^2 - |A - B - T|^2) / 4

51

    ;;    = A * B - T * B

52

    lpm     tmp, Z+

53

    sub     AB0, tmp

54

    lpm     tmp, Z

55

    sbc     AB1, tmp

56

57

    ;; AB += T * B  als  AB += T ? B : 0  -->  AB = A * B

58

    brtc    2f

59

    add     AB0, B0

60

    adc     AB1, zero

61

2:

62

    ret

63

64

.macro  quads from=0, to

65

    .word \from * \from

66

    .if   \to-\from

67

      quads "(\from+1)",\to

68

    .endif

69

.endm

70

71

.section .progmem.data.lupo, "a", @progbits

72

73

lupo:

74

    quads 0, 99

75

    quads 100, 199

76

    quads 200, 255

77

78

.global  mul_quad

79

.type  mul_quad, @function

80

.size  mul_quad, .-mul_quad

81

82

.type  lupo, @object

83

.size  lupo, .-lupo

- Mit Alignment auf eine Pagegrenze würde man 2 Zyklen sparen.

- Bei der min(a,b) Sonderbehandlung kann man noch einmal einen Takt 
sparen, wenn man den kompletten Code mit dem anderen Register 
dupliziert.

- Wenn man die Sonderbehandlung für das LSB weglässt, und mit etwas 
Ungenauigkeit leben kann, dann spart man noch einmal 5 Zyklen (worst 
case).

Macht 28, bzw 23 Taktzyklen.

Insgesamt leidet der Code aber darunter, dass der AVR keine indirekt 
indizierten Speicherzugriffe kann (wie LDA $mmmm,y oder lda ($zz),y beim 
6502). Ist halt RISC.

Hier merkt man wieder, das RISC und 8 Bit nicht so recht miteinander 
wollen. Die Addressberechnungen benötigen zu viele Befehle. Streng 
genommen wird auch der Orthogonalitätsgedanke hinter RISC nicht 
eingehalten, da Address- und Datenbusbreite unterschiedlich sind.

Johann L. schrieb:
> Man mag einwenden, dass die verwendete Formel zu kompliziert sei und nur
> einem Theoretiker einfallen könne, aber immerhin läuft das 20% schneller
> als die oben gesetzten 40 Ticks — und zwar in der Praxis :-)

Wie man an den C64 Beispielen sieht, ist es alles andere als Theorie. Da 
hat sich eher der Urheber des Spruches oben disqualifiziert.

dbdb schrieb:
> Kranker scheiss, Respekt!

+1

dbdb schrieb:
> Kranker scheiss, Respekt!
>
> Wie kommt man auf die Formel? Die quadratische Identität kann ich ja
> noch nachvollziehen, aber das ist ja echt hammer

Die Formel fasst nur in der üblichen, international verständlichen und 
Hardware- und Software-unabhängigen Notation zusammen, was die konkrete 
Implementierung macht.  Wenn man will, kann man damit nachrechnen, dass 
die Implementierung das macht, was sie machen soll.

Bei der Grundidee ab = u² - v² hat man u = (a+b)/2 und v = (a-b)/2 wobei 
es unterschiedliche Strategien gibt, sich der lästigen halbzahligen 
Werte zu entledigen.

Da aus einer Tabelle gelesen wird, in der Worte (16-Bit Werte) stehen 
und Adressen Byte-Adressen sind, muss man für den Zugriff auf die 
Tabelle mit 2 multiplizieren.  Diese Multiplikation hebt sich mit den 
Divisionen durch 2 auf -- zumindest fast, denn von ungeraden (relativ 
zum Tabellenanfang) Byte-Adressen zu lesen liefert natürlich Unsinn.

Also setzt man das untere Bit des Offsets einfach auf 0 (per ANDI ~1). 
Falls a+b (und damit auch a-b) gerade ist, hat sich dadurch nix 
geändert.  Ist es jedoch ungerade, so liefert unbeirrtes Weitermachen 
das Ergebnis für u = (a+b-1)/2 und v = (a-b-1)/2.  Nachrechnen ergibt, 
dass der Wert dann um b zu klein ist, d.h. für diesen Fall addiert man 
noch b zum Ergebnis.  Das a * b erhält man also mit u² - v² + tb mit u = 
(a+b-t)/2, v = (a-b-t)/2 und t (per BST gespeichert im T-Flag) ist eben 
a+b mod 2.

Weil man die Adresse zum Wert von (a+b-t)/2 bereits hat (Tabellenanfang 
+ a+b-t) und man danach den Wert für (a-b-t)/2 braucht (also Tab + 
a-b-t) subtrahiert man von der alten Adresse 2b und erhält so die 
Adresse für v.

Was verbleibt ist das Problem mit dem Vorzeichen von a-b.  a+b-t ist 
davon unberührt, und a+b-t - 2min(a,b) ist a-b-t wenn a >= b und b-a-t 
wenn b >= a.  In Software ist das so gelöst dass b nach Berechnung von 
a+b mit a überschrieben wird falls a kleiner ist als b.  Dieses 
Überschreiben drückt sich in der Formel dadurch aus, dass an allen 
Stellen, wo b nach den Überschrieb verwendet wird, ein min(a,b) 
auftaucht.


Possetitjel schrieb:

> Andererseits ist (a-b)² gleich a² - 2ab + b².
>
> Der Ausdruck a² - (a-b)² + b² liefert somit, wie man durch
> Nachrechnen leicht verifiziert, direkt 2ab. a² und b² lassen
> sich, da sowohl a als auch b durch 8 Bit darstellbar sind,
> ebenfalls ohne tieferes Nachdenken in der Tabelle nachschlagen.
> Wenn man von links nach rechts auswertet, tritt auch kein
> Unterlauf auf. Der mögliche Überlauf kann durch den ohnehin
> notwendigen Rechts-Shift behandelt werden.
>
> Der einzige Nachteil besteht darin, dass drei Tabellenzugriffe
> notwendig sind (und natürlich muss durch eine Verzweigung ganz
> am Anfang die Bedingung a>b sichergestellt werden).
>
> Mich würde ein Geschwindigkeitsvergleich zu Deiner Variante
> interessieren, aber da ich AVR-Assembler fast nicht lesen
> und gar nicht programmieren kann, bleibt dies vermutlich ein
> frommer Wunsch...

Rein intuitiv würde ich sagen, dass das für AVR nix bringt, eben wegen 
des bereits von Tim angemerkten Rumgehudels mit den Adressberechnungen. 
Zudem kostet es bereits 6 Ticks, einen 16-Bit Wert aus dem Flash zu 
lesen.  Und da ist noch kein einziger Befehl Adressberechnung dabei.

Auf 65xx Assembler bin ich zu lange raus; da ist der Nerv-Faktor eher 
dass der nur 3 Register hat.  Aber was Tim gezeigt hat wird sich nur 
sehr schwer -- falls überhaupt -- unterbieten lassen.

Wer dran interessiert ist kann ja weiter dran rumspielen, ich jedenfalls 
brauche so eine Routine nicht, weder für AVR noch für 65xx noch für 
sonst einen µC.

Tim    schrieb:
> Insgesamt leidet der Code aber darunter, dass der AVR keine indirekt
> indizierten Speicherzugriffe kann (wie LDA $mmmm,y oder lda ($zz),y beim
> 6502). Ist halt RISC.

Viele RISCs besitzen die Fähigkeit, die Adresse von Speicherzugriffen 
als Summe zweier Register zu bilden.