Rechnung gibt falsches Ergebnis.

Finnegan

Weitere Beobachtung: Problem tritt beim MinGW.org-Compiler nur mit -O0 auf. Mit -O2 gibt auch der zweimal die 243 aus.

Habe den Code etwas umgeschrieben damit die relavenaten Stellen im Assembler-Code besser zu finden sind, und der code nicht zu einem write(stdout, "242 243") "optimiert" wird:

test3.cpp:

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

using namespace std;

int test(int y)
{
    y -= pow(3, 5);
    return y;
}

int main()
{
    volatile int do_not_fold_this_constant = 486; 
    int y = do_not_fold_this_constant;
    cout << test(y) << " " << pow(3,5);
}

Bei diesem Code tritt das Problem ebenfalls auf.

Erzeugter Assembler-Code (nur den der relevante test()-Funktion, sind sonst 12.000 Zeilen):

g++ -O0 -masm=intel -fverbose-asm -Wa,-adhln -g test3.cpp (ohne Optimierungen, gibt 242 aus):

   6:test3.cpp     **** int test(int y)
   7:test3.cpp     **** {
  48              		.loc 1 7 1
  49              		.cfi_startproc
  50 0000 55       		push	ebp	 #
  51              		.cfi_def_cfa_offset 8
  52              		.cfi_offset 5, -8
  53 0001 89E5     		mov	ebp, esp	 #,
  54              		.cfi_def_cfa_register 5
  55 0003 83EC28   		sub	esp, 40	 #,
  56              	 # test3.cpp:8:     y -= pow(3, 5);
   8:test3.cpp     ****     y -= pow(3, 5);
  57              		.loc 1 8 13
  58 0006 C7442404 		mov	DWORD PTR [esp+4], 5	 #,
  58      05000000 
  59 000e C7042403 		mov	DWORD PTR [esp], 3	 #,
  59      000000
  60 0015 E8000000 		call	__ZSt3powIiiEN9__gnu_cxx11__promote_2IT_T0_NS0_9__promoteIS2_XsrSt12__is_integerIS2_E7__value
  60      00
  61              	 # test3.cpp:8:     y -= pow(3, 5);
  62              		.loc 1 8 7
  63 001a DB4508   		fild	DWORD PTR [ebp+8]	 # y
  64 001d DEE1     		fsubrp	st(1), st	 #,
  65 001f D97DF6   		fnstcw	WORD PTR [ebp-10]	 #
  66 0022 0FB745F6 		movzx	eax, WORD PTR [ebp-10]	 # tmp88,
  67 0026 80CC0C   		or	ah, 12	 # tmp88,
  68 0029 668945F4 		mov	WORD PTR [ebp-12], ax	 #, tmp88
  69 002d D96DF4   		fldcw	WORD PTR [ebp-12]		 #
  70 0030 DB5D08   		fistp	DWORD PTR [ebp+8]	 # y
  71 0033 D96DF6   		fldcw	WORD PTR [ebp-10]		 #
  72              	 # test3.cpp:9:     return y;
   9:test3.cpp     ****     return y;
  73              		.loc 1 9 12
  74 0036 8B4508   		mov	eax, DWORD PTR [ebp+8]	 # _8, y
  75              	 # test3.cpp:10: }
  76              		.loc 1 10 1
  77 0039 C9       		leave	
  78              		.cfi_restore 5
  79              		.cfi_def_cfa 4, 4
  80 003a C3       		ret

g++ -O2 -masm=intel -fverbose-asm -Wa,-adhln -g test3.cpp (mit Optimierungen, gibt 243 aus):

  6:test3.cpp     **** int test(int y)
   7:test3.cpp     **** {
  81              		.loc 2 7 1 is_stmt 1 view -0
  82              		.cfi_startproc
   8:test3.cpp     ****     y -= pow(3, 5);
  83              		.loc 2 8 5 view LVU3
  84              	 # test3.cpp:7: {
   7:test3.cpp     ****     y -= pow(3, 5);
  85              		.loc 2 7 1 is_stmt 0 view LVU4
  86 0010 83EC0C   		sub	esp, 12	 #,
  87              		.cfi_def_cfa_offset 16
  88              	 # test3.cpp:8:     y -= pow(3, 5);
  89              		.loc 2 8 7 view LVU5
  90 0013 D97C2406 		fnstcw	WORD PTR [esp+6]	 #
  91 0017 DB442410 		fild	DWORD PTR [esp+16]	 # y
  92 001b D8250400 		fsub	DWORD PTR LC0	 #
  92      0000
  93              	LVL2:
   9:test3.cpp     ****     return y;
  94              		.loc 2 9 5 is_stmt 1 view LVU6
  95              	 # test3.cpp:8:     y -= pow(3, 5);
   8:test3.cpp     ****     y -= pow(3, 5);
  96              		.loc 2 8 7 is_stmt 0 view LVU7
  97 0021 0FB74424 		movzx	eax, WORD PTR [esp+6]	 # tmp90,
  97      06
  98 0026 80CC0C   		or	ah, 12	 # tmp90,
  99 0029 66894424 		mov	WORD PTR [esp+4], ax	 #, tmp90
  99      04
 100 002e D96C2404 		fldcw	WORD PTR [esp+4]		 #
 101 0032 DB1C24   		fistp	DWORD PTR [esp]	 # %sfp
 102 0035 D96C2406 		fldcw	WORD PTR [esp+6]		 #
 103              	LVL3:
   8:test3.cpp     ****     y -= pow(3, 5);
 104              		.loc 2 8 7 view LVU8
 105 0039 8B0424   		mov	eax, DWORD PTR [esp]	 # y, %sfp
 106              	 # test3.cpp:10: }
  10:test3.cpp     **** }
 107              		.loc 2 10 1 view LVU9
 108 003c 83C40C   		add	esp, 12	 #,
 109              		.cfi_def_cfa_offset 4
 110              	LVL4:
 111              		.loc 2 10 1 view LVU10
 112 003f C3       		ret

(Leider nicht ganz so hübsch wie mit Godbolt ;-))

Das muss ich mir nochmal genauer ansehen, um zu verstehen, was da jetzt genau passiert. Wahrscheinlich schau ich da morgen nochmal drüber

hustbaer

Wie ist denn pow üblicherweise implementiert? Können CPUs das direkt? Ich würde annehmen dass das über log-mul-exp geht. Und dann würde es mich nicht wundern wenn selbst bei einfachen, schönen Inputs wie zwei kleinen Integern nicht immer exakt das erwartete Ergebnis rauskommt.

hustbaer

Wie auch immer, es sollte helfen explizit zu runden:

y -= round(pow(3, 5));

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