Effiziente Division durch Multiplikation

[Walter T. (nicolas)]

Guten Abend zusammen,

ich will auf dem ARM Cortex M3/M4 eine Division durch eine effizientere 
Multiplikation ersetzen. Bei dem erwarteten Wertebereich habe ich so 
keine Genauigkeitsverluste zu verzeichnen. Momentan bin ich mir jedoch 
unsicher, ob mir die Integer-Promotion bei negativem Zähler keinen 
Strick drehen kann.

Gegeben ist ein Bruch a/b. b ändert sich selten. b ist positiv. Der 
Übersicht halber habe ich alles in eine Funktion geschrieben:

int32_t divideByB(int32_t a)
{
    // Vorbereitung
    const uint32_t b = 12345L;
    assert( b!=0 );
    uint32_t oneDividedByB_q32;
    if( b == 1 )
        oneDividedByB_q32 = UINT32_MAX; // Naeherung
    else
        oneDividedByB_q32 = (1ULL<<32)/b; // Exakt

    // Berechnung
    int64_t result = (int64_t) a * oneDividedByB_q32; // Kann nicht ueberlaufen, da 31 Bit * 32 Bit
    return result/(1LL<<32); // Nur vorderes Register zurueckgeben
}

Die Hilfsvariable wurde q32 gewählt, damit das fertige Ergebnis direkt 
in einem Register steht und der Bitshift erst gar nicht ausgeführt 
werden muß. Bei der Berechnung der Hilfsvariablen wird der Nenner zu 
UINT64 befördert und die Division ausgeführt (was lange dauert - hier 
aber nicht weiter stört).

a ist Int64. oneDividedByB_q32 passt in ein Int64, wird also zu einem 
Int64 befördert und die Multiplikation ausgeführt. Wenn b!=1, ist das 
Ergebnis immer identisch mit der direkten Division.

Habe ich das richtig zusammengepuzzlet, oder habe ich irgendeinen 
möglichen Überlauf übersehen?

Viele Grüße
W.T.

[Joe F. (easylife)]

Macht natürlich nur Sinn, wenn dein Divisor immer gleich ist, und 
"Vorbereitung" nicht jedesmal berechnet werden muss.

Statt

return result/(1LL<<32);

würde ich eher

return (result >> 32);

hinschreiben sonst baut dir der Compiler u.U. hier noch eine Division 
ein.

[Dr. Sommer (Gast)]

Und das ist wirklich schneller als die 12 Takte, die Cortex-M3/4 für 
eine Integer-Division braucht?

[holger (Gast)]

Schiebeoperationen würde ich bei signed Typen
grundsätzlich nicht machen.

[Jim M. (turboj)]

Joe F. schrieb:
> return (result >> 32);
>
> hinschreiben sonst baut dir der Compiler u.U. hier noch eine Division
> ein.

Blöd nur das Bitshifts auf negativen Integers undefiniertes Verhalten in 
C sind. Ich verstehe den OP so, dass er das nicht komplett ausschließen 
kann.

[LIO (Gast)]

Warnt der Compiler nicht davor?

[avr (Gast)]

Dr. Sommer schrieb:
> Und das ist wirklich schneller als die 12 Takte, die Cortex-M3/4 für
> eine Integer-Division braucht?

Wenn der Compiler das entsprechend umsetzt, braucht das vier Takte.

[Walter T. (nicolas)]

Mein

Joe F. schrieb:
> würde ich eher
>
> return (result >> 32);
>
> hinschreiben

Hmneee, das Ergebnis wäre implementation defined. Ich will das 
Vorzeichenbit durchaus mit Sicherheit behalten.

Dr. Sommer schrieb:
> Und das ist wirklich schneller als die 12 Takte, die Cortex-M3/4 für
> eine Integer-Division braucht?

Unmittelbar: Kaum. Mittelbar kann ich im Faktor oneDividedByX_q32 
allerdings das Ergebnis einer längeren Rechnung zwischenspeichern, und 
dann will ich den Teil mit der Multiplikation auch "noch eben" 
mitnehmen, weil ich mir dann in der "schnellen" Funktion die 
Fallunterscheidung mit O sparen kann.

[Walter T. (nicolas)]

Bitte noch einmal zurück zu meiner Frage:

1. Ich bin mir bei der Multiplikation Uint32_t/Int32 unsicher, ob ich 
die Integer-Promotion-Regeln richtig ausgelegt habe - kann ich hier 
immer sicher sein, daß nichts überlaufen kann?

2. Läßt sich das Vorgehen auf einen Divisor mit Vorzeichen erweitern?

[Dr. Sommer (Gast)]

Walter T. schrieb:
> Bei dem erwarteten Wertebereich habe ich so
> keine Genauigkeitsverluste zu verzeichnen.

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

int32_t divideByB(int32_t a, uint32_t b)
{
    // Vorbereitung
    assert( b!=0 );
    uint32_t oneDividedByB_q32;
    if( b == 1 )
        oneDividedByB_q32 = UINT32_MAX; // Naeherung
    else
        oneDividedByB_q32 = (1ULL<<32)/b; // Exakt

    // Berechnung
    int64_t result = (int64_t) a * oneDividedByB_q32; // Kann nicht ueberlaufen, da 31 Bit * 32 Bit
    return result/(1LL<<32); // Nur vorderes Register zurueckgeben
}


int main (void) {
  int32_t a = 10;
  uint32_t b = 5;
  printf ("a/b=%" PRId32 "\ndivideByB(%" PRId32 ", %" PRIu32") = %" PRId32 "\n", a/b, a, b, divideByB (a, b));
}

ergibt:

a/b=2
divideByB(10, 5) = 1

Hmm...

[Walter T. (nicolas)]

Dr. Sommer schrieb:
> a/b=2
> divideByB(10, 5) = 1Hmm...

Und das Schlimme: Ich verstehe nicht warum.

int32_t divideByB(int32_t a)
{
    // Vorbereitung
    const uint32_t b = 5L;
    assert( b!=0 );
    uint32_t oneDividedByB_q32;
    if( b == 1 )
        oneDividedByB_q32 = UINT32_MAX; // Naeherung, Liefert immer a - 1 zurueck
    else
        oneDividedByB_q32 = (1ULL<<32)/b; // Exakt

    printf("oneDividedByB = %u ", oneDividedByB_q32);

    // Berechnung
    int64_t result = (int64_t) a * (uint64_t) oneDividedByB_q32; // Kann nicht ueberlaufen, da 31 Bit * 32 Bit
    printf("zwischen = %lli ", result);
    printf("1<<32=%lli ", 1LL<<32);
    return (result/(1LL<<32)); // Nur vorderes Register zurueckgeben
}

int main(void)
{
    printf("result=%li", divideByB(10));
    return 1;
}

oneDividedByB = 858993459 zwischen = 8589934590 1<<32=4294967296 result=1

[Dr. Sommer (Gast)]

Walter T. schrieb:
> Und das Schlimme: Ich verstehe nicht warum.
Ist doch klar, bei der ersten Division wird abgerundet, weil (1<<32) 
nicht glatt durch 5 teilbar ist. Somit ist das Ergebnis der 
Multiplikation etwas zu klein und wird wieder abgerundet. Man könnte die 
beiden Divisionen korrekt rundend ausführen aber auch dann müssten sich 
Werte finden lassen bei denen es nicht passt.
Wenn du weniger als 22 signifikante Binär-Stellen hast, könntest du 
float verwenden, da dauert die Multiplikation 1 Takt.

[avr (Gast)]

Du musst vor dem Shift, oder vor der Multiplikation runden durch 
Addition eines Offsets. Der shift rundet nämlich ab. Hab mir jetzt aber 
keine Gedanken über den negativen Bereich gemacht. Ich hab die Methode 
früher für vorzeichenlose Divisionen auf avrs benutzt.

[Martin S. (docmartin)]

Tja, ein typischer Rundungsfehler.
4294967296 * 2 = 8589934592 ...

d.h., dein Ergebnis ist 1,9999999995343387126922607421875 und
für Integer wird einfach alles nach dem Komma entfernt --> 1

Ahoi, Martin

[Walter T. (nicolas)]

Das hat sich gerade gedoppelt:

2*4294967296 = 8589934592

d.h. hier wird korrekt abgerundet. Habe es gerade selbst herausgefunden.

[Walter T. (nicolas)]

avr schrieb:
> oder vor der Multiplikation runden durch Addition eines Offsets.

Mit dem Offset riskiere ich natürlich wieder einen Überlauf, wenn ich a 
nicht begrenze. Für den einfachen Fall mit x/b mit int32_t x, b ist die 
Multiplikation mit dem Kehrwert dann einfach zu kompliziert, weil 
Fallunterscheidungen nötig sind.

Schade - aber nicht schlimm. Wie ja schon oben richtig geschrieben 
wurde, kostet mich die nackte Division zweier 32-Bit-Zahlen gerade mal 
maximal 12 Takte.


Generell suche ich auch noch immer nach einer sinnvollen allgemeinen 
Lösung zum allgemeineren Problem

int32_t r, a, z, n

[...]

r = (int64_t) a*z/n;

(siehe Beitrag "Einflußfaktoren zur Rechendauer bei Divisionen"), aber bislang 
konnte ich es jetzt immer irgendwie umgehen.

[Johann L. (gjlayde)]

Walter T. schrieb:
> Generell suche ich auch noch immer nach einer sinnvollen allgemeinen
> Lösung zum allgemeineren Problem

Ist nicht so ganz trivial, zumindest wenn man sich den Code im GCC dazu 
anschaut.

Und einfach mit dem "Kehrwert" zu multiplizieren ist auch nicht die 
Lösung, zumindest wenn das Ergebnis das gleiche wie bei einer Division 
sein soll.  Eine unsigned Division durch 7 übersetzt gcc zum Beispiel 
so:

div_7:
  ldr  r3, .L3
  umull  r2, r1, r0, r3
  sub  r0, r0, r1
  add  r0, r1, r0, lsr #1
  lsr  r0, r0, #2
  bx  lr
  .align  2
.L3:
  .word  613566757

Es gibt also noch zusätzliche Arithmetik, um das Ergebnis anzupassen.

Ohne diese Anpassung wird dein Ergebnis i.d.R. nicht mit dem Quotienten 
einer entsprechenden Division übereinstimmen, und deine Hausaufgabe ist 
dann, den (maximalen) Fehler abzuschätzen und zu entscheiden ob dieser 
in der Anwendung tolerierbar ist.

Was du momentan versuchst ist eine Art JIT Ansatz, jedoch ohne über den 
eigentlichen Algorithmus zu verfügen.  Und selbst mit dessen Kenntnis 
müsste immer noch eine Verzweigung anhand des Divisors geschehen, der 
zusätzliche Laufzeit kostet.  Oder der Algorithmus wied so allgemein 
gehelten, dass er für alle Divisoren gültig ist, was u.U. einen 
deitlich größeren Overhead bedeutet.

Überhaupt scheint das nur sinnvoll, wenn alle Divisoren bekannt sind, 
und dann ein switch-case eine schnelle Verzweigung zum Code hinbekommt. 
Oder wenn alle Divisoren eine gemeinsame Struktur haben, welche sich 
dann in einem einheitlichen Algorithmus niederschlägt.

[avr (Gast)]

Walter T. schrieb:
> Mit dem Offset riskiere ich natürlich wieder einen Überlauf, wenn ich a
> nicht begrenze.

Der M4 hat doch auch saturierende Arithmetik. Auf dem AVR hab ich damals 
übrigens meine Algorithmen für jede Eingabe auf dem PC getestet, um 
sicher zu gehen, dass sie funktionieren.

[Walter T. (nicolas)]

Mit 8 Stunden Abstand fällt mir auf: Eigentlich ist das Ergebnis gar 
nicht so schlecht. Über den kompletten Wertebereich liegt der Fehler im 
Vergleich zur Rechnung mit double-Genauigkeit zwischen -1.5 und +0.5. Im 
16-Bit-Bereich sogar zwischen -1 und 0.

Bei 10 und 5 sieht es nur drastisch falsch aus. Und exakt ist es 
natürlich wirklich nicht.