Effiziente Division durch Multiplikation

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.

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.

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.

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 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

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

Das hat sich gerade gedoppelt:

2*4294967296 = 8589934592

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

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.

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.

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.