Rechnung gibt falsches Ergebnis.

[Finnegan]

@Swordfish sagte in Rechnung gibt falsches Ergebnis.:

Mingw-w64 unter MSYS2: https://abload.de/img/untitledvvjye.png

Das hat mich aber jetzt interessiert:

MinGW heruntergeladen und installiert.

$ g++ --version
g++.exe (MinGW.org GCC Build-2) 9.2.0
Copyright (C) 2019 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

$ g++ -dumpmachine
mingw32

$ cat test.cpp
#include<iostream>
#include<math.h>

using namespace std;

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

$ g++ test.cpp -o test

$ file test.exe
test.exe: PE32 executable for MS Windows (console) Intel 80386 32-bit

$ ./test.exe

242 243

LOL ... trägt nicht gerade dazu bei, mein Vertrauen in das originale MinGW zu wiederherzustellen, besonders da diese Werte als FP exakt darstellbar sind (wie alle Integer, die nicht grösser als die Mantisse sind), wie @wob bereits ausgeführt hat

Edit: Ein 32-Bit-Problem scheint es auch nicht zu sein (MinGW-w64/MSYS2):

> i686-w64-mingw32-g++ --version                                                                                                                       Mon Oct 26 16:02:23 2020
i686-w64-mingw32-g++.exe (dev-tools-v2) 9.2.0
Copyright (C) 2019 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

> i686-w64-mingw32-g++ -dumpmachine                                                                                                                    Mon Oct 26 16:03:30 2020
i686-w64-mingw32

> i686-w64-mingw32-g++ test.cpp -o test

> file test.exe
test.exe: PE32 executable (console) Intel 80386, for MS Windows

> ./test.exe
243 243⏎

[Th69]

@Finnegan: Scheint doch ein 32-bit Problem zu sein, denn bei 64-bit kommt doch das richtige Ergebnis 243 heraus?!

Was mich jedoch immer noch wundert, ist das falsche Ergebnis von 242, denn dann müßte ja bei (int)pow(3, 5) als Ergebnis 244 herauskommen, so daß die Subtraktion dann 242 ergibt. Ist da evtl. ein Compilerfehler bei der Subtraktion?

Probiert mal folgenden Code aus (ich habe kein MinGW installiert):

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

[Finnegan]

@Th69 sagte in Rechnung gibt falsches Ergebnis.:

@Finnegan: Scheint doch ein 32-bit Problem zu sein, denn bei 64-bit kommt doch das richtige Ergebnis 243 heraus?!

Für meinen MinGW-w64-Test im Nachtrag habe ich mit i686-w64-mingw32-g++ einen 32-Bit Compiler benutzt. Das ist wie -m32, nur dass ich hier ein selbstgebautes Cross-Compiler Setup verwende (ohne Multilib und sowas, jeder Compiler und die Libs haben nen eigenen Unterverzeichnisbaum).

Was mich jedoch immer noch wundert, ist das falsche Ergebnis von 242, denn dann müßte ja bei (int)pow(3, 5) als Ergebnis 244 herauskommen, so daß die Subtraktion dann 242 ergibt. Ist da evtl. ein Compilerfehler bei der Subtraktion?

Probiert mal folgenden Code aus (ich habe kein MinGW installiert):

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

MinGW.org (32-Bit):

$ g++ --version
g++.exe (MinGW.org GCC Build-2) 9.2.0
Copyright (C) 2019 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

$ cat test2.cpp
#include<iostream>
#include<math.h>

using namespace std;

int main()
{
    int y = 486;
    double d = pow(3, 5);
    int z = d;
    y = y - z;
    cout << y << " " << pow(3,5);
}

$ g++ test2.cpp -o test2

$ ./test2.exe
243 243

MinGW-w64 gibt hier in 64- und 32-Bit dasselbe aus.

Da fragt sich: Warum? Sollte nicht auch y -= pow(3, 5) exakt die selbe Konvertierung durchführen? Das ist schon echt seltsam. Riecht daher eher nach einem Compilerproblem. Auch std::cout << (int) pow(3, 5) gibt 243 aus.

[Th69]

Kannst du den generierten Assemblercode vom Original-Code posten?

[Finnegan]

@Th69 sagte in Rechnung gibt falsches Ergebnis.:

Kannst du den generierten Assemblercode vom Original-Code posten?

Vom fehlerhaften? Ja, muss grad für ne Weile weg, aber später am Abend finde ich sicher Zeit dafür.

[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";

242

vs.

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

243

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 [...]