Forum: Mikrocontroller und Digitale Elektronik Effiziente Division durch Multiplikation


Announcement: there is an English version of this forum on EmbDev.net. Posts you create there will be displayed on Mikrocontroller.net and EmbDev.net.
von Walter T. (nicolas)


Bewertung
0 lesenswert
nicht lesenswert
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:
1
int32_t divideByB(int32_t a)
2
{
3
    // Vorbereitung
4
    const uint32_t b = 12345L;
5
    assert( b!=0 );
6
    uint32_t oneDividedByB_q32;
7
    if( b == 1 )
8
        oneDividedByB_q32 = UINT32_MAX; // Naeherung
9
    else
10
        oneDividedByB_q32 = (1ULL<<32)/b; // Exakt
11
12
    // Berechnung
13
    int64_t result = (int64_t) a * oneDividedByB_q32; // Kann nicht ueberlaufen, da 31 Bit * 32 Bit
14
    return result/(1LL<<32); // Nur vorderes Register zurueckgeben
15
}

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.

: Bearbeitet durch User
von Joe F. (easylife)


Bewertung
0 lesenswert
nicht lesenswert
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.

: Bearbeitet durch User
von Dr. Sommer (Gast)


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

von holger (Gast)


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

von Jim M. (turboj)


Bewertung
0 lesenswert
nicht lesenswert
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.

Beitrag #5298648 wurde von einem Moderator gelöscht.
von LIO (Gast)


Bewertung
0 lesenswert
nicht lesenswert
Warnt der Compiler nicht davor?

von avr (Gast)


Bewertung
0 lesenswert
nicht lesenswert
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.

von Walter T. (nicolas)


Bewertung
0 lesenswert
nicht lesenswert
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.

von Walter T. (nicolas)


Bewertung
0 lesenswert
nicht lesenswert
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?

: Bearbeitet durch User
von Dr. Sommer (Gast)


Bewertung
0 lesenswert
nicht lesenswert
Walter T. schrieb:
> Bei dem erwarteten Wertebereich habe ich so
> keine Genauigkeitsverluste zu verzeichnen.
1
#include <stdio.h>
2
#include <assert.h>
3
#include <stdint.h>
4
#include <inttypes.h>
5
6
int32_t divideByB(int32_t a, uint32_t b)
7
{
8
    // Vorbereitung
9
    assert( b!=0 );
10
    uint32_t oneDividedByB_q32;
11
    if( b == 1 )
12
        oneDividedByB_q32 = UINT32_MAX; // Naeherung
13
    else
14
        oneDividedByB_q32 = (1ULL<<32)/b; // Exakt
15
16
    // Berechnung
17
    int64_t result = (int64_t) a * oneDividedByB_q32; // Kann nicht ueberlaufen, da 31 Bit * 32 Bit
18
    return result/(1LL<<32); // Nur vorderes Register zurueckgeben
19
}
20
21
22
int main (void) {
23
  int32_t a = 10;
24
  uint32_t b = 5;
25
  printf ("a/b=%" PRId32 "\ndivideByB(%" PRId32 ", %" PRIu32") = %" PRId32 "\n", a/b, a, b, divideByB (a, b));
26
}

ergibt:
1
a/b=2
2
divideByB(10, 5) = 1
Hmm...

von Walter T. (nicolas)


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

Und das Schlimme: Ich verstehe nicht warum.
1
int32_t divideByB(int32_t a)
2
{
3
    // Vorbereitung
4
    const uint32_t b = 5L;
5
    assert( b!=0 );
6
    uint32_t oneDividedByB_q32;
7
    if( b == 1 )
8
        oneDividedByB_q32 = UINT32_MAX; // Naeherung, Liefert immer a - 1 zurueck
9
    else
10
        oneDividedByB_q32 = (1ULL<<32)/b; // Exakt
11
12
    printf("oneDividedByB = %u ", oneDividedByB_q32);
13
14
    // Berechnung
15
    int64_t result = (int64_t) a * (uint64_t) oneDividedByB_q32; // Kann nicht ueberlaufen, da 31 Bit * 32 Bit
16
    printf("zwischen = %lli ", result);
17
    printf("1<<32=%lli ", 1LL<<32);
18
    return (result/(1LL<<32)); // Nur vorderes Register zurueckgeben
19
}
20
21
int main(void)
22
{
23
    printf("result=%li", divideByB(10));
24
    return 1;
25
}
1
oneDividedByB = 858993459 zwischen = 8589934590 1<<32=4294967296 result=1

: Bearbeitet durch User
von Dr. Sommer (Gast)


Bewertung
0 lesenswert
nicht lesenswert
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.

von avr (Gast)


Bewertung
0 lesenswert
nicht lesenswert
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.

von Martin S. (docmartin)


Bewertung
0 lesenswert
nicht lesenswert
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

von Walter T. (nicolas)


Bewertung
0 lesenswert
nicht lesenswert
Das hat sich gerade gedoppelt:

2*4294967296 = 8589934592

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

von Walter T. (nicolas)


Bewertung
0 lesenswert
nicht lesenswert
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
1
int32_t r, a, z, n
2
3
[...]
4
5
r = (int64_t) a*z/n;
(siehe Beitrag "Einflußfaktoren zur Rechendauer bei Divisionen"), aber bislang 
konnte ich es jetzt immer irgendwie umgehen.

von Johann L. (gjlayde) Benutzerseite


Bewertung
0 lesenswert
nicht lesenswert
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:
1
div_7:
2
  ldr  r3, .L3
3
  umull  r2, r1, r0, r3
4
  sub  r0, r0, r1
5
  add  r0, r1, r0, lsr #1
6
  lsr  r0, r0, #2
7
  bx  lr
8
  .align  2
9
.L3:
10
  .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.

von avr (Gast)


Bewertung
0 lesenswert
nicht lesenswert
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.

von Walter T. (nicolas)


Bewertung
0 lesenswert
nicht lesenswert
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.

: Bearbeitet durch User

Antwort schreiben

Die Angabe einer E-Mail-Adresse ist freiwillig. Wenn Sie automatisch per E-Mail über Antworten auf Ihren Beitrag informiert werden möchten, melden Sie sich bitte an.

Wichtige Regeln - erst lesen, dann posten!

  • Groß- und Kleinschreibung verwenden
  • Längeren Sourcecode nicht im Text einfügen, sondern als Dateianhang

Formatierung (mehr Informationen...)

  • [c]C-Code[/c]
  • [code]Code in anderen Sprachen, ASCII-Zeichnungen[/code]
  • [math]Formel in LaTeX-Syntax[/math]
  • [[Titel]] - Link zu Artikel
  • Verweis auf anderen Beitrag einfügen: Rechtsklick auf Beitragstitel,
    "Adresse kopieren", und in den Text einfügen




Bild automatisch verkleinern, falls nötig
Bitte das JPG-Format nur für Fotos und Scans verwenden!
Zeichnungen und Screenshots im PNG- oder
GIF-Format hochladen. Siehe Bildformate.
Hinweis: der ursprüngliche Beitrag ist mehr als 6 Monate alt.
Bitte hier nur auf die ursprüngliche Frage antworten,
für neue Fragen einen neuen Beitrag erstellen.

Mit dem Abschicken bestätigst du, die Nutzungsbedingungen anzuerkennen.