Forum: Mikrocontroller und Digitale Elektronik Integer-Arithmetik - mal wieder


von Walter T. (nicolas)


Angehängte Dateien:

Lesenswert?

Hallo zusammen,

ich will einen 64-Bit-Integer Skalieren, d.h. mit einem ganzzahligen 
Bruch z/n multiplizieren. Der Bruch z/n bewegt sich in der Größenordnung 
von 1, also etwa [0.01 .... 100]. Zähler und Nenner können allerdings 
sehr groß werden.

Weil ich später mit Differenzen weiterarbeite, will ich keinen 
Auflösungsverlust zu den Enden hin. Mit dem sich zwangsläufig ergebenden 
Aliasing/Moiré-Effekt/Jitter oder wie auch immer das in diesem 
Zusammenhang heisst kann ich leben. Also ist Integer-Rechnung anstelle 
double gefragt.

Auf dem ARM Cortex M4 stellt der ARM-GCC keinen 128-Bit-Datentyp zur 
Verfügung. Also läuft es darauf hinaus, die Arithmetik mit 
Mehrwort-Operationen nachzubilden.

Lange Rede kurzer Sinn: Ich habe Integer-Arithmetik implementiert 
(Funktion int64_t scale_s64(int64_t in, int32_t z, int32_t n) ) und bin 
mir nicht sicher, ob ich irgendeinen merkwürdigen Fall (neudeutsch "edge 
case") übersehen habe.

Ganz glücklich bin ich auch mit der Behandlung des Absolutwerts von 
INT_MIN nicht. Die Vergleiche will der Compiler nämlich nicht 
wegoptimieren. Aber mir fällt keine bessere Variante ein.

: Bearbeitet durch User
von Peter D. (peda)


Lesenswert?

Walter T. schrieb:
> Auf dem ARM Cortex M4 stellt der ARM-GCC keinen 128-Bit-Datentyp zur
> Verfügung.

Was sind das denn für super duper genaue Berechnungen, die Du anstellen 
mußt.
Mir reichen in der Regel 24Bit Mantisse (einfaches float) völlig aus.
Mit integer Formaten hatte ich mir früher schon öfter ein Ei gelegt 
(Überläufe, Rundungsfehler).

von Walter T. (nicolas)


Lesenswert?

Peter D. schrieb:
> Was sind das denn für super duper genaue Berechnungen, die Du anstellen
> mußt.

Die sind gar nicht so genau. Sprich: Die Auflösung ist begrenzt. Nur der 
Wertebereich von 32 Bit ist nicht ausreichend und Rundungsfehler dürfen 
sich auch nicht akkumulieren.

Die 52 Bit eines double würden mir wohl reichen.

Mein Plan ist also: Integer-Implementierung und dann 
Geschwindigkeitsvergleich mit double. Aber dafür muss die 
Integer-Variante erst einmal fehlerfrei sein.

: Bearbeitet durch User
von Peter D. (peda)


Lesenswert?

Walter T. schrieb:
> Die 52 Bit eines double würden mir wohl reichen.

Was für extreme Anwendungen sind das denn, die eine so hohe Genauigkeit 
benötigen. Machst Du astronomische Berechnungen?

von Walter T. (nicolas)


Lesenswert?

Peter D. schrieb:
> Machst Du astronomische Berechnungen?

Reine Spielerei. Wir haben Sonntag.

von W.S. (Gast)


Lesenswert?

Peter D. schrieb:
> Mir reichen in der Regel 24Bit Mantisse (einfaches float) völlig aus.

Eben. Mit single hat man etwas mehr als 7 Stellen ist man zumeist gut 
genug bedient. Naja, außer beim Frequenzzähler, wo man das GHz bis auf's 
Hz auflösen will.

Aber der Walter hat mehr als alle anderen hier recht unübliche Probleme, 
die ihm selbst nicht geheuer vorkommen. So wie dieses hier.

Walter T. schrieb:
> ich will einen 64-Bit-Integer Skalieren, d.h. mit einem ganzzahligen
> Bruch z/n multiplizieren. Der Bruch z/n bewegt sich in der Größenordnung
> von 1, also etwa [0.01 .... 100].

Wenn dir double ausreicht (wie du weiter oben schriebest), dann nimm 
auch double. Etwas sinnvolleres kann man dir nicht anraten, es sei denn, 
du willst  dir selber eine 128 Bit Integerarithmetik (in Assembler) 
schreiben.

W.S.

von Walter T. (nicolas)


Angehängte Dateien:

Lesenswert?

Okay...das ist peinlich. Ich habe gerade eine Pause gemacht, und dann 
ist mir eine deutlich einfachere Lösung eingefallen.

Die ganze 
Schul-schriftlich-multiplizieren-addieren-Dividieren-Geschichte lässt 
sich anscheinend auf wenige (wahrscheinlich vier) Zeilen C eindampfen. 
Jetzt muss ich nur noch beweisen, dass der Wertebereich an keiner Stelle 
überschritten wird.
Aber das ist nichts mehr für heute abend. Die Modultests haben zumindest 
schon einmal grünes Licht gegeben.

Vorgehensweise: Erst die Division und dann die Korrektur über das 
Restglied. Dann ist die Behandlung des Vorzeichens noch der einzige 
aufwendige Teil. Da wurmt mich der extra-Vergleich für den INT-MIN Fall 
immer noch.

Edit: Der Nachweis war doch nicht so knifflig. Der ging auch nebenbei 
auf dem WC. Durch den Nachweis sind es jetzt ein paar Zeilen mehr 
geworden, aber kein Vergleich zu vorher. Um so mehr wurmt mich der 
Umgang mit INT_MIN.

(der letzte Cast im Quelltext ist falsch. Ist natürlich int64_t. Kann 
die Datei nur nicht mehr ändern.)

: Bearbeitet durch User
von Walter T. (nicolas)


Angehängte Dateien:

Lesenswert?

Wenn man das Ganze mit vorzeichenbehafteten Integer durchzieht, wird es 
auch nicht kompakter. __builtin_smulll_overflow() ist ganz schön teuer.

Ich hoffe, ich habe hier keinen letalen Fall übersehen.

von Walter T. (nicolas)


Lesenswert?

Da ist noch ein Bug drinne, aber der ist so offensichtlich, dass er 
sofort auffällt, wenn man nicht gerade früh am Morgen noch keinen Kaffee 
hatte.

von Walter T. (nicolas)


Lesenswert?

Ich habe mal ein paar kleine Benchmarks gemacht, indem ich die 
Funktionen mit jeweils 100000 mal mit Zufallsdaten  für die drei 
Eingangsgrößen gefüttert habe.

 - Die Multiwort-Implementierung ist erwartungsgemäß ziemlich langsam. 
1032 Takte pro Durchlauf werden durchschnittlich benötigt.
 - Die unsigned-Variante benötigt durchschnittlich 360 Takte
 - Die signed-Variante benötigt durchschnittlich 530 Takte (*natürlich 
nach Fehlerbehebung, s.o.)
 - Eine double-Variante benötigt durchschnittlich 1480 Takte.
jeweils nach dem Abzug des overheads auf einem STM32F446RE bei 168 MHz.

Der Zufallszahlengenerator rand() scheint das absolute Gegenteil von dem 
zu sein, was man landläufig darunter versteht, und genau für den Zweck 
gemacht. Bei jeweils drei Durchläufen habe ich taktgenau immer das 
Gleiche herausbekommen. (Interrupts waren natürlich aus.)

Das Ergebnis hat mich ein wenig überrascht. Ich hätte nicht gedacht, 
dass double langsamer als die langsame Multiwort-Variante ist. Die 
unsigned-Variante kann trotz der länglichen Fallunterscheidungen am 
Anfang punkten. Die __builtin_mul_overflow() & Co. sind doch recht 
teuer.

von (prx) A. K. (prx)


Lesenswert?

Zu rand() gehört zwingend srand() zur einmaligen Initialisierung beim 
Start des Programms - und zwar mit unterschiedlichen Werten.

Pass auf, dass du die Laufzeit von rand() nicht mitmisst.

: Bearbeitet durch User
von Walter T. (nicolas)


Lesenswert?

(prx) A. K. schrieb:
> Pass auf, dass du die Laufzeit von rand() nicht mitmisst.

Natürlich messe ich sie mit. Aber sie ist ja bei allen Durchgängen 
taktgenau gleich (was mich auch gewundert hat).

Die Zahlen sind Netto mit Overhead (168 Takte) herausgerechnet.

In diesem Fall will ich ja gar keinen unterschiedlichen Seed. Mir kommt 
die deterministische Folge sogar sehr entgegen für einen fairen 
Vergleich.

: Bearbeitet durch User
von (prx) A. K. (prx)


Lesenswert?

Walter T. schrieb:
> Natürlich messe ich sie mit. Aber sie ist ja bei allen Durchgängen
> taktgenau gleich (was mich auch gewundert hat).

Das wundert mich nicht. Algorithmen für Pseudo-Zufallszahlen können 
durchaus eine exakt konstante Laufzeit haben. Beispielsweise wenn sowas 
verwendet wird:
https://de.wikipedia.org/wiki/Linear_rückgekoppeltes_Schieberegister

von Philipp Klaus K. (pkk)


Lesenswert?

(prx) A. K. schrieb:
> Zu rand() gehört zwingend srand() zur einmaligen Initialisierung beim
> Start des Programms - und zwar mit unterschiedlichen Werten.

Laut C-Standard verhält dich rand() so, als gäbe es beim Start des 
Programms ein implizites srand(1)¹.
Dass rand() eine reproduzierbare Sequenz von Zufallszahlen liefert ist 
beim Debuggen hilfreich.

¹ Man kann sich hier allerdings leider nicht ganz darauf verlassen, da 
OpenBSD hier bewusst vom C-Standard abweicht, und unter Umständen echte 
Zufallszahlen liefert, wenn vor dem Aufruf von rand() kein srand() 
explizit aufgerufen wurde.

Beitrag #6811927 wurde vom Autor gelöscht.
von Walter T. (nicolas)


Lesenswert?

Ich schrieb ja oben, die __builtin_XXX_overflow() funktionen seien 
teuer. Andererseits war an der einfachen 96-Bit-Arithmetik ja durchaus 
noch Optimierungspotenzial.

Also habe ich mal
 - Bei der Bildung des Arrays die Endianness berücksichtigt, dass die 
Byte-Reihenfolge im Array zum normalen Layout passt.
 - Bei der Division loop unrolling betrieben. Das spart eine Division 
und eine Restwert-Bildung.
 - Ein bischen Kleinkram.

Jetzt bin ich bei 323 Takten pro Durchlauf. So macht die Sache Spaß.

Ein bischen neugierig bin ich noch, ob sich noch ein ganzes Stück 
einsparen lässt, wenn die Skalierfaktoren jeweils in int16_t passen (tun 
sie in meinem Fall). Die zwei 64-Bit-Divisionen sind jetzt der teuerste 
Teil.

: Bearbeitet durch User
von Walter T. (nicolas)


Lesenswert?

Nachtrag: Wenn man sich auf Skalierungsfaktoren Z/N jeweils 16 Bit 
beschränkt, kommt man auf 138 Takte, allerding dann wirklich mit loop 
unrolling bis zur Unleserlichkeit. Wahrscheinlich könnte ein fähiger 
Assemblist noch ein paar Takte rausholen, aber ich bin jetzt zufrieden.

Im Vergleich zu double mit 1480 Takten pro Konvertierung ist das ganz 
brauchbar und man hat die gleiche Auflösung ueber den gesamten Bereich.

Bitte melde dich an um einen Beitrag zu schreiben. Anmeldung ist kostenlos und dauert nur eine Minute.
Bestehender Account
Schon ein Account bei Google/GoogleMail? Keine Anmeldung erforderlich!
Mit Google-Account einloggen
Noch kein Account? Hier anmelden.