Rechnung gibt falsches Ergebnis.

Finnegan

@hustbaer sagte in Rechnung gibt falsches Ergebnis.:

Wie ist denn pow üblicherweise implementiert? Können CPUs das direkt?

So halb. X86 kennt f2xm1 ( $2^x-1$ ) und fyl2xp1 ( $y \cdot \log_{2}(x + 1)$ ), allerdings mit recht eingeschränkten Werten, die die Operanden haben dürfen. Ich vermute die gängigen pow-Implementierungen nutzen das mit einiges an Code drumherum, was natürlich auch unterschiedliches Verhalten je nach Compiler/Standardlib begünstigt.

@hustbaer sagte in Rechnung gibt falsches Ergebnis.:

Wie auch immer, es sollte helfen explizit zu runden:

Hilft, wenn es denn ein Rundungsproblem wäre, was wir wohl alle zuerst dachten. Letztendlich stellte sich aber heraus, dass der MinGW.org-Compiler bei:

int y = 486;
y -= pow(3, 5);

und

int y = 486;
double d = pow(3, 5);
int z = d;
y = y - z;

zwei unterschiedliche Ergebnisse in y liefert. Das darf eigentlich nicht passieren, egal wie viel Rundungsfehler da drin steckt (allerdings nur mit -O0, mit -O2 siehts gut aus).

hustbaer

@Finnegan sagte in Rechnung gibt falsches Ergebnis.:

Hilft, wenn es denn ein Rundungsproblem wäre, was wir wohl alle zuerst dachten. Letztendlich stellte sich aber heraus, dass der MinGW.org-Compiler bei:
int y = 486;
y -= pow(3, 5);
und
int y = 486;
double d = pow(3, 5);
int z = d;
y = y - z;
zwei unterschiedliche Ergebnisse in y liefert.

Ja. Und?

y -= pow(3, 5)
// ==
y = y - pow(3, 5);
// ==
y = truncate(double(y) - pow(3, 5));

Hier findet die Konvertierung durch Abschneiden der Nachkommastellen nach der Subtraktion statt.
Also bei z.B. y = 486 und pow(3, 5) == 243.1 kommt erstmal 242.9 raus und das wird zu 242 abgeschnitten.

int y = 486;
double d = pow(3, 5);
int z = d;
y = y - z;

Hier wird erstmal 243.1 zu 243 beschnitten, und dann 486 - 243 = 243 gerechnet.

Ist doch alles ganz normal.

ps:
Also nur damit das klar ist: y -= z; entspricht bei eingebauten Typen y = y - z;. D.h. wenn y ein int ist und z ein double, dann gelten die ganz normalen Regeln wie überall sonst. D.h. y wird erstmal nach double konvertiert, und dann die Subtraktion double - double durchgeführt. Danach kommt dann die Zuweisung für die das Ergebnis - durch Abschneiden der Kommastellen - wieder nach int konvertiert wird.

Anderer Code der was anderes macht, anderes Ergebnis. Verstehe die Verwunderung nicht.

Th69

Daran hatte ich überhaupt nicht (mehr) gedacht.

Also muß man explizit einen Cast durchführen, damit eine Integer-Operation (für die Subtraktion) benutzt wird:

y -= static_cast<int>(round(pow(3, 5))); // oder alternativ hier: (int)

Finde ich nicht sehr intuitiv und ich sehe auch keinen Mehrwert in der Verwendung der Fließkommaoperation hier (wenn sowieso das Ergebnis wieder in eine Ganzzahl zurückgeschrieben wird).

PS: Seit Jahren programmiere ich meist in C# und da muß man explizit den Cast durchführen, sonst kommt ein Compilerfehler:

error CS0266: Cannot implicitly convert type 'double' to 'int'. An explicit conversion exists (are you missing a cast?)

Finnegan

@hustbaer sagte in Rechnung gibt falsches Ergebnis.:

Also nur damit das klar ist: y -= z; entspricht bei eingebauten Typen y = y - z;.

Oh, Mist. Ich hab nix gesagt ... von der kompakten -=-Schreibweise verleiten lassen und dann nicht mehr hinterfragt.

Irgendwas seltsames geht allerdings trotzdem vor sich:

int y = 486;
double dy = y;
y = dy - pow(3, 5);
cout << y << "\n";

vs.

int y = 486;
double dy = y;
dy = dy - pow(3, 5);
y = dy;
cout << y << "\n";

Oder hab ich hier auch was übersehen?

Auch interessant:

#include<iostream>
#include<math.h>

using namespace std;

int main()
{
    cout << boolalpha << ((486 - pow(3, 5)) == (486 - pow(3, 5))) << "\n";
}

MinGW.org: false, MinGW-w64: true

Alles mit -O0.

hustbaer

@Th69 sagte in Rechnung gibt falsches Ergebnis.:

Finde ich nicht sehr intuitiv und ich sehe auch keinen Mehrwert in der Verwendung der Fließkommaoperation hier (wenn sowieso das Ergebnis wieder in eine Ganzzahl zurückgeschrieben wird).

Naja. Es macht die Sprachdefinition einfacher wenn man sagt a @= b entspricht bei eingebauten Typen a = a @ b mit @ ist eins aus +, - etc. Weiters glaube ich würde es einiges an Verwunderung stiften wenn sie sich unterschiedlich verhielten.

Th69

Mir geht es hierbei um die automatische Promotion von int a zu double.
Wenn man einen eigenen Datentyp hat (der intern int verwendet!) und dort dann den -=-Operator mittels

type operator -= (const type &a, int b)
{
  a.value -= b;
}

implementiert, so wird ja bei

type a(486);

a -= pow(3,5);

zuerst die Promotion von double nach int ausgeführt und dann der Operator damit aufgerufen (d.h. der Datentyp von a wird nicht verändert).

Swordfish

@Finnegan sagte in Rechnung gibt falsches Ergebnis.:

Auch interessant:

#include<iostream>
#include<math.h>

using namespace std;

int main()
{
    cout << boolalpha << ((486 - pow(3, 5)) == (486 - pow(3, 5))) << "\n";
}

MinGW.org: false, MinGW-w64: true
Alles mit -O0.

DAS macht mir Angst!?

Finnegan

@Swordfish sagte in Rechnung gibt falsches Ergebnis.:

MinGW.org: false, MinGW-w64: true
Alles mit -O0.

DAS macht mir Angst!?

Ja, wenn ich nicht vorher schon was wichtiges übersehen hätte, hätte ich auch laut "Compiler/Standardbibliothek-Bug!" geschrieen

Da gibt es aber nicht viel dran zu rütteln. Ich vermute, dass die pow-Implementierung irgendwie zustandsbehaftet ist - vielleicht eine gobale Variable oder FPU-Flags, die nicht korrekt wiederhergestellt werden.

Ich vermute es ist eher die pow-Implementierung als ein GCC-Fehler.

Edit: Sieht aus als sei das eine MinGW-eigene Implementierung und zwar vermutlich diese hier (Version 5.4.1, die ich hier verwende):

https://osdn.net/projects/mingw/scm/git/mingw-org-wsl/blobs/5.4-trunk/mingwrt/mingwex/math/pow_generic.sx

Das ist eine etwas umfangreichere Assembler-Implementierung und mir defintitiv gerade zu fummelig, die zu debuggen. Mein Bauchgefühl sagt mir aber da ist irgendwo die Ursache zu finden.

*john 0

@Swordfish sagte in Rechnung gibt falsches Ergebnis.:

DAS macht mir Angst!?

Wahrscheinlich wird hier im 32Bit Modus noch das 80Bit Float Format der Intel CPUs (das gab es sonst nur noch bei Motorola 68k und 88k) verwendet, und mit dem 64Bit Modus wird auf SSE/AVX gesetzt was nur IEEE Single und Double kennt.

Finnegan

@john-0 sagte in Rechnung gibt falsches Ergebnis.:

@Swordfish sagte in Rechnung gibt falsches Ergebnis.:

DAS macht mir Angst!?

Wahrscheinlich wird hier im 32Bit Modus noch das 80Bit Float Format der Intel CPUs (das gab es sonst nur noch bei Motorola 68k und 88k) verwendet, und mit dem 64Bit Modus wird auf SSE/AVX gesetzt was nur IEEE Single und Double kennt.

Selbst ein zwölffingriger Kobold im Gehäuse, der mit einem Duodezimalsystem-Rechenschieber arbeitet, muss stets zum selben Ergebnis kommen wenn er exakt die selbe Berechnung durchführt. Was immer hier sonst noch seltsam läuft, hier gibt es mindestens einen Bug.

*john 0

@Finnegan sagte in Rechnung gibt falsches Ergebnis.:

Selbst ein zwölffingriger Kobold im Gehäuse, der mit einem Duodezimalsystem-Rechenschieber arbeitet, muss stets zum selben Ergebnis kommen wenn er exakt die selbe Berechnung durchführt.

Aber das ist nicht der Fall. Es wird im 80Bit Float Format anders gerundet als im IEEE Modus!

Finnegan

@Th69 sagte in Rechnung gibt falsches Ergebnis.:

type operator -= (const type &a, int b)

Ui. Das ist ein gutes Argument. Für einen überladenen Operator ist -= eine Funktion f(int) in diesem Fall. Aus der Perspektive wäre es natürlich auch konsistent, zuerst das Argument nach int zu konvertieren.

Allerdings kann man auch einen int::operator-=(double) implementieren, der die Berechnung so durchführt, dass sie sich exakt wie x = x - y verhält. Das ist hier bei den eigebauten Typen vom Konzept her wohl Fall:

https://en.cppreference.com/w/cpp/language/operator_assignment:

Builtin compound assignment
[...]
The behavior of every builtin compound-assignment expression E1 op= E2 [...] is exactly the same as the behavior of the expression E1 = E1 op E2 [...]

Finnegan

@john-0 sagte in Rechnung gibt falsches Ergebnis.:

Aber das ist nicht der Fall. Es wird im 80Bit Float Format anders gerundet als im IEEE Modus!

Ja, das ist aber eine "andere Berechnung" ich rede hier von "exakt der selben Berechnung" in (486 - pow(3, 5)) == (486 - pow(3, 5)). Die rechte und die linke Seite von == werden bei deinem Argument entweder beide in 80-Bit oder beide in IEEE durchgeführt. Der Vergleich muss also immer true liefern. Tut er aber nicht bei MinGW.org.

*john 0

@Finnegan sagte in Rechnung gibt falsches Ergebnis.:

Ja, das ist aber eine "andere Berechnung" ich rede hier von "exakt der selben Berechnung" in (486 - pow(3, 5)) == (486 - pow(3, 5)). Die rechte und die linke Seite von == werden bei deinem Argument entweder beide in 80-Bit oder beide in IEEE durchgeführt. Der Vergleich muss also immer true liefern. Tut er aber nicht bei MinGW.org.

Nein, so etwas kann man bei Floats leider nicht aussagen. Es gibt etliche Faktoren, die das Ergebnis einer Float Berechnung verändern können, und solange man nicht explizit das IEEE Normverhalten aktiviert hat, kann der Compiler einen ziemlichen Unfug treiben. D.h. exakt die gleichen Rechnungen müssen explizit nicht das Gleiche ergeben – sie können. Lies Dir dazu die Doku von Intel, die IEEE Norm und das Compilerhandbuch durch.

Finnegan

@john-0 sagte in Rechnung gibt falsches Ergebnis.:

[...] D.h. exakt die gleichen Rechnungen müssen explizit nicht das Gleiche ergeben – sie können. Lies Dir dazu die Doku von Intel, die IEEE Norm und das Compilerhandbuch durch.

Auf verschiedenen CPUs, Betriebssystemen oder auch mit anderen Compiler-Flags gebe ich dir da ohne Widerspruch recht. Diese Probleme sind mir bekannt, wodurch man selbst auf AMD- und Intel-x86 mit IEEE-Floats auf standardkonforme Weise leicht unterschiedliche Ergebnisse erhalten kann.

Wir reden hier aber von zwei identischen Berechnungen auf der selben CPU, mit den selben Compiler-Flags auf dem selben OS im selben Programm, ja sogar im selben Ausdruck. Ich kenne die Normen nicht auswendig, aber ich habe starke Zweifel, dass sich speziell dieses Szenario mit den Freiheiten die diese einer Implementation lassen erklären lässt.

Meines wissens sind Float-Berechnungen schon deterministisch, aber eben nicht sonderlich portabel.

SeppJ

Wieso ist das so überraschend für dich? Angenommen das Szenario mit dem 80Bit/64Bit Rechenwerk. Dann ist eine ganz valide Interpretation des Codes folgende Arbeitsanweisung (die auch sehr wahrscheinlich bei O0 genommen wird):

Berechne die linke Seite (80 Bit Präzision)
Speichere das Ergebnis irgendwo (wird zu 64 Bit gerundet), denn wir brauchen Platz im Rechenwerk für die nächste Rechnung
Berechne die rechte Seite (80 Bit Präzision)
Oh, nun soll ich die beiden Werte vergleichen! Der 2. Wert steht ja noch im Rechenwerk (80 Bit), nun muss ich nur noch den 1. Wert wieder in das andere Register vom Rechenwerk laden.
Jetzt wird der 80 Bit-Wert mit dem gerundeten 64-Bit-Wert verglichen.
-> false

Quiche Lorraine

@SeppJ
Ich poste mal diesbezüglich einen passenden Stackoverflow Artikel:

https://stackoverflow.com/questions/7517588/different-floating-point-result-with-optimization-enabled-compiler-bug

Finnegan

@SeppJ sagte in Rechnung gibt falsches Ergebnis.:

Wieso ist das so überraschend für dich?

Danke für die Aufklärung. Das ist in der Tat überraschend für mich, dass hier verglichen wird, ohne vorher sicherzustellen, dass der zweite Wert ebenfalls gerundet wird.

Ich habe mich bisher tatsächlich darauf verlassen, dass ich das exakt selbe Ergebnis erhalte, wenn ich identische Rechenanweisungen in meinem Code durchführe. Egal wie die interne Repräsentation in den CPU-Registern aussieht und wie viele Rundungsfehler meine Berechnungen einbringen.

Ich bin sogar entsetzt, wenn das wirklich alles so Standardkonform sein sollte. Das heisst wenn ich mich darauf verlassen muss, das zwei indentische FP-Berechnugen zum bitgenau selben Ergebnis kommen, verwende ich besser eine wohldefinierte FP-Emulation via Integer - auch wenn ich dasselbe Ergebnis im selben Programm mit exakt den selben Compiler-Flags und ordentlich initialisierter FPU benötige?

Edit: Oder eben dem Compiler verklickern, solche Spässe zu unterlassen: -ffloat-store sorgt hier tatsächlich dafür, dass auch MinGW.org true ausgibt. Danke für den SO-Link!

*john 0

@Finnegan sagte in Rechnung gibt falsches Ergebnis.:

Meines wissens sind Float-Berechnungen schon deterministisch, aber eben nicht sonderlich portabel.

Das ist nur dann garantiert, wenn man explizit das IEEE Normverhalten durch den Compiler aktiviert, und einen Rattenschwanz an Nebenbedingungen einhält. Das wird schon problematisch, wenn ein Autovektorisierer am Werk ist, und wird richtig problematisch, wenn Threading im Spiel ist. Denn die Reihenfolge von Additionen oder Multiplikationen spielt bei der Reproduzierbarkeit der Ergebnisse eine nicht unwesentliche Rolle. Da es so etwas wie out-of-order execution gibt, kann das auch eine Rolle spielen. Was genau der Compiler hier treibt ist unklar, die Masse wird hier true zurückliefern, wie man das erwartet.

Finnegan

@john-0 sagte in Rechnung gibt falsches Ergebnis.:

@Finnegan sagte in Rechnung gibt falsches Ergebnis.:

Meines wissens sind Float-Berechnungen schon deterministisch, aber eben nicht sonderlich portabel.

Das ist nur dann garantiert, wenn man explizit das IEEE Normverhalten durch den Compiler aktiviert, [...]

Jo danke, wieder was dazugelernt. Das mit der 80bit-Repräsentation und IEEE waren mir zwar bekannt, aber ich dachte bisher tatsächlich, dass sich das nur auswirkt, wenn man zwei mit unterschiedlichen Compiler-Flags gebaute Programme hat, oder sogar andere CPUs/Architekturen. Separate Threads hätte ich auch noch verstanden, wenn sie auf verschiedenen Kernen laufen. Aber im selben Ausdruck finde ich schon erschreckend

Nicht leicht, wenn man z.B. bei einer komplexeren physikalischen Simulation mit identischen Eingabewerten exakt den selben Endzustand reproduzieren will.