Dividieren durch Multiplizieren

[Mampf F. (mampf)]

Guten Morgen,

da das Thema ja ein bisserl Mathe ist, hab ich es hier gepostet :)

Ich hab einen Dividierer (im FPGA) und der ist langsam (ca 42 
clock-zyklen @ 100MHz).

Da immer nur durch 27 dividiert wird, könnte man auf die Idee kommen, 
mittels Festpunktarithmetik mit dem Kehrwert zu multiplizieren.

Mich würde nun interessieren, ob jemand weiß, wieviele 
"Nachkommastellen" man benötigt, damit das Ergebnis vor'm Komma nach der 
Multiplikation immer korrekt ist.

Beispiel 27/27:

27*2**16 * trunc((1/27 * 2**16)) / (2**32) = 0,9998931884765625

Sieht soweit nicht schlecht aus, aber runden müsste man wohl noch.

Was denkt ihr? Reichen 16Bit?

Viele Grüße,
Mampf

[A. S. (Gast)]

Wie ist dein Eingangswertebereich? (16 Bit? 0... 65535?)

Wie die Ausgabe? Also auch wieviel Bits und was bedeutet 1lsb.

Dann kannst du ausrechnen, wie du ohne Fehler hinkommst.

[A. S. (Gast)]

Uns nicht trunc, sondern aufrunden, zumindest bei Insignien. Sonst 
machst Du zweimal trunc

[Mampf F. (mampf)]

A. S. schrieb:
> Wie ist dein Eingangswertebereich? (16 Bit? 0... 65535?)

Das sind 32Bit. Also quasi 32+16Bit Nachkomma wären dann 48Bit.

> Wie die Ausgabe? Also auch wieviel Bits und was bedeutet 1lsb.

Die ist auch 32Bit. Den Nachkomma-Teil würde ich wegschneiden (nach dem 
Runden) und würde dann durch Subtraktion den Rest berechnen (brauch 
leider beides). Die Frage nach dem LSB verstehe ich nicht ganz.

Im Prinzip soll ein 32Bit int reingehen und ein 32Bit int rauskommen mit 
5bit Rest.

A. S. schrieb:
> Uns nicht trunc, sondern aufrunden, zumindest bei Insignien. Sonst
> machst Du zweimal trunc

Du meinst das trunc aus dem Beispiel? Das hab ich nur verwendet, damit 
es mir die Nachkomma-Stellen der Nachkomma-Stellen abschneidet. Sonst 
wär 1 rausgekommen.

Vermutlich geht bei Fixpunkt runden genauso wie bei 32bit ints? Also 0,5 
addieren und dann trunc?

[Mampf F. (mampf)]

Ah wieder was gelernt ... Der Nachkomma-Teil muss mindestens 32Bit 
haben, damit das klappt.

leider scheint __uint128_t mit dem GCC nicht zu funktionieren, sodass 
ich das nicht mal testen kann :/

[Mampf F. (mampf)]

So jetzt hab ich es ...

Also ... die Anzahl der Bits für den Divisor müssen 32 + 5 + 1 sein 
(32Bit Input, 6Bit Divisor (5 für 0-27 + 1Bit für das +0,5) mit einer 
Nachkommastelle)

Man addiert nach der Multiplikation mit 2**30 (eben diese 0,5) und 
truncated.

Dann kommt das richtige raus :-)

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <math.h>
int main() {
    __uint128_t shift2 = 36;
    __uint128_t divisor = (__uint128_t) (((__uint128_t)1<<shift2) / 27.0);
//    printf("%032llx\n", (long long unsigned int) divisor);
    
    for (int i=0;i<100000;i++) {
        __uint128_t dividend = rand();
        
        __uint128_t result = dividend * divisor;
        result += ((__uint128_t)1<<30);
        result >>= (__uint128_t) (shift2);

        uint32_t cmpres = (uint32_t) dividend / 27;
    
        if ((uint32_t) result != cmpres) {
            printf("error\n");
            printf("%d %d\n", (int) result, cmpres);
        }
    }    
    return 0;
}

[Yalu X. (yalu) (Moderator)]

Du kannst auch den GCC fragen, wie er das Problem lösen würde:

  y = x / 27;

wird übersetzt zu:

  movl  x(%rip), %ecx
  movl  $795364315, %edx
  movl  %ecx, %eax
  mull  %edx
  subl  %edx, %ecx
  shrl  %ecx
  addl  %ecx, %edx
  shrl  $4, %edx
  movl  %edx, y(%rip)

Die magische Konstante 795364315 ist dabei gleich ceil(5 * 2**32 / 27).
Rückübersetzt nach C entspricht dieser Code dem folgenden:

  uint32_t tmp = (uint64_t)x * 795364315 >> 32;
  y = (tmp + ((x - tmp) >> 1)) >> 4;

Im Vergleich zu deiner Lösung ist hier die magische Konstante nur 30
statt 32 Bit groß, dafür werden aber insgesamt mehr Rechenoperationen
benötigt. Auf dem PC ist deine Lösung vermutlich die effizientere (die
128 Bit sind unnötig, 64 Bit für das Multiplikationsergebnis und 32 Bit
für den Rest sollten genügen). Welche Lösung auf dem FPGA (wo die Shifts
nichts kosten und die Addierer exakt auf die jeweils benötigte Bitzahl
optimiert werden können) bzgl. der Anzahl der Logikblöcke die bessere
ist, müsste noch untersucht werden.

Wie kommst du eigentlich in

        result += ((__uint128_t)1<<30);

auf die 1<<30?

Damit kompensierst du wohl die Tatsache, dass divisor auf Grund der
endlichen Bitzahl etwas zu klein ist, denn ohne diese Korrektur wären
die Ergebnisse teilweise falsch.

Dur erwähnst in deinen Beiträgen auch immer wieder das Thema Runden.
Aber dein Verfahren liefert ja – genauso wie auch der C-Ausdruck
dividend/27 – immer abgerundete Ergebnisse. Damit kann der Wert 1<<30
also nichts zu tun haben, zumal für klassisches Runden (d.h. Aufrunden
bei einem Nachkommaanteil von ≥0.5) dieser Wert nicht 1<<30, sondern
1<<35 sein müsste.

[Mampf F. (mampf)]

Yalu X. schrieb:
> Dur erwähnst in deinen Beiträgen auch immer wieder das Thema Runden.
> Aber dein Verfahren liefert ja – genauso wie auch der C-Ausdruck
> dividend/27 – immer abgerundete Ergebnisse. Damit kann der Wert 1<<30
> also nichts zu tun haben, zumal für klassisches Runden (d.h. Aufrunden
> bei einem Nachkommaanteil von ≥0.5) dieser Wert nicht 1<<30, sondern
> 1<<35 sein müsste.

Glaube das ist nicht ganz richtig ... Da der Divisor ja schon 
ge-truncated wird, muss man runden, da 27/27 sonst zB 0.999999 geben 
könnte (und auch tut).

Das 1<<30 kommt wohl daher:

Divisor muss 36bit haben ... 32Bit Größe, weil der Dividend so groß ist, 
5Bit Größe, weil der Divisor so groß ist und noch 1Bit Nachkomma des 
Divisors zum Runden.

Dieses letzte Bit ist dann genau auf der Stelle 1<<30 und ohne das sind 
oft die Ergebnisse falsch.

Mathematisch kann ich das nicht beweisen, aber es gibt ein konsistentes 
Muster.

Aber geniale Idee, den GCC zu befragen, wie er das machen würde.

In der Realität hat bei mir der Dividend 40Bit, der Fixpunkt-Divisor 
damit 46Bit und die +1 kommt auf das Bit 40.

Da reichen dann leider die 64Bit nicht mehr aus - aber dem FPGA ist das 
eh fast egal :)

[A. S. (Gast)]

Machen wir es einfacher:

Schreib Du hin, welche Division du machen möchtest, und ich schreib das 
als Multiplikation.

[Mampf F. (mampf)]

A. S. schrieb:
> Machen wir es einfacher:
>
> Schreib Du hin, welche Division du machen möchtest, und ich schreib das
> als Multiplikation.

Können wir probieren :)

trunc(dividend / divisor). Der Dividend hat 40Bit, der Divisor 5Bit 
(fester Wert 27).

Im FPGA läuft es schon. Allerdings schafft das FPGA 40*46Bit 
Multiplikation nicht ohne 1-stufige Pipeline @ 100MHz. Der eine 
Clock-Cycle ist nicht so wichtig, aber wenn es ohne geht, wäre natürlich 
besser :)

[A. S. (Gast)]

Mampf F. schrieb:
> 40Bit

Also Werte von 1 bis 1 Billionen?

Und das Ergebnis dann Werte von 1 bis 30E6?

[Yalu X. (yalu) (Moderator)]

Mampf F. schrieb:
> trunc(dividend / divisor). Der Dividend hat 40Bit, der Divisor 5Bit
> (fester Wert 27).

So gehts mit nur einer Multiplikation ohne zusätzliche Addition:

$$y=\left\lfloor\frac {kx}{2^{43}}\right\rfloor$$

wobei

$$k=\left\lceil\frac{2^{43}}{27}\right\rceil=325781223045$$

Der Faktor k ist 39 Bit lang.

Folgendes Testprogramm sucht das kleinste x, für das der obige Ausdruck
ein falsches Ergebnis liefert:

#include <stdio.h>
#include <stdint.h>

#define N 43
#define K 325781223045  // = ⌈2⁴³/27⌉

int main(void) {
  for(uint64_t x=0; ; x++) {
    uint64_t y1 = x / 27;
    uint64_t y2 = ((__uint128_t)x * K) >> N;
    if(y1 != y2) {
      printf("Fehler bei x=%lu y1=%lu y2=%lu\n", x, y1, y2);
      break;
    }
  }
}

Ausgage (nach ca. 20 Minuten):

Fehler bei x=1256584717466 y1=46540174720 y2=46540174721

Somit rechnet das Verfahren bis einschließlich x=1256584717465 korrekt,
was oberhalb der gewünschten 2⁴³-1=1099511627775 liegt.

> Im FPGA läuft es schon. Allerdings schafft das FPGA 40*46Bit
> Multiplikation nicht ohne 1-stufige Pipeline @ 100MHz. Der eine
> Clock-Cycle ist nicht so wichtig, aber wenn es ohne geht, wäre natürlich
> besser :)

Die obige Multiplikation lässt sich in 16 Additionen/Subtraktionen
zerlegen, die in 4 sequentiell aufeinanderfolgenden Ebenen ausgeführt
werden können. Die gesamte Multiplikation benötigt also ohne Pipelining
das Vierfache der Zeit einer einzelnen Addition/Subtraktion.

[Nils P. (torus)]

Im Buch Hackers Delight ist der Algorithmus zum errechnen der Inversen 
für konstante Divisionen dokumentiert.

Auf der Homepage kann man sich auch die Konstanten direkt errechnen 
lassen.

 https://www.hackersdelight.org/
 https://www.hackersdelight.org/magic.htm

[Mampf F. (mampf)]

Yalu X. schrieb:
> wobei
> k=⌈24327⌉=325781223045
> k=\left\lceil\frac{2^{43}}{27}\right\rceil=325781223045
>
> Der Faktor k ist 39 Bit lang.

Ahja, das Ceiling ist interessant!

Wenn ich die Konstante um 1 erhöhe, dann muss ich später dieses 1<<40 
addieren nicht mehr machen.

Ansonsten sieht es schon ähnlich aus zu dem, was ich gemacht habe ... 
Größtenteils anders geschrieben :)

Vielen Dank!

[Michael W. (Gast)]

Mampf F. schrieb:
> Da immer nur durch 27 dividiert wird, könnte man auf die Idee kommen,
> mittels Festpunktarithmetik mit dem Kehrwert zu multiplizieren.

Das ist aber jetzt nicht dein Ernst? Das  ist das erste, was man tut, um 
sich Divisionsrechenleistung zu sparen.

> Wenn ich die Konstante um 1 erhöhe, dann muss ich später dieses 1<<40
> addieren nicht mehr machen.
Das tut indirekt der Compiler für dich und ob das die jeweils optimale 
Lösung fürs Runden in FPGAs ist, sei dahin gestellt.

> Ansonsten sieht es schon ähnlich aus zu dem, was ich gemacht habe ...
> Größtenteils anders geschrieben :)
Nun ja.