It's not a compiler bug, it's not a problem in C/C++, or even in IEEE-754 floating point math. The MinGW64 thread library simply forgot to initialize the FPU correctly.
He'd have had the same problem if he wrote his code in assembly language and then ran it in a MinGW64 thread.
TBH I accredit dumb luck more than anything to finding the cause of this. It took many hours.
The equation is the commutative property, not associative. IIRC addition in IEEE-754 is commutative.
The property that fails is (A+B)+C = A+(B+C).
I didn't want to go too far into the background of floating point math in the post as it wasn't ultimately the issue, just something I was mindful of.
> Update: Thanks to gus_massa and wiml @ HackerNews for pointing out I had used associative instead of commutative.... This means that floating point arithmetic is noncommutative [i]n that A + B != B + A.
Floating-point arithmetic is commutative, but not associative. So for floating-point numbers,
A + B = B + A always, but
A + (B + C) != (A + B) + C in many cases.
The problem comes about from significant digits and rounding - A + B might round up to 1, but B + C might round down to zero, etc. The point remains that the order in which you do things actually does matter in floating-point arithmetic, but it's the order you compute the "inner parentheses," not the actual order of numbers in (e.g.) the FPU registers.
The C standard library has functions to do this configuration; search for fesetenv. _fpreset is not a standard function but, at least in Wine's implementation [1], just resets the relevant configuration registers to an initial state.
[1] https://github.com/wine-mirror/wine/blob/e909986e6ea5ecd49b2...
If the thread package was not setting up the "saved" state for a new thread before switching to it then that state could be not what was desired. The best thing to do would probably be to copy the current FP state of the parent at the moment the thread is created.
On a long series of calculations changing the rounding mode could well be enough to cause the difference in results the poster saw.
I would tend to the opinion that redoing the calculation with different rounding modes is in fact a not bad way to figure out how reliable the results are in the first place.
(the rounding mode is not the only thing affected by not saving and restoring the FP state correctly)
Though in general, IEEE-754 math does have a small amount of global state, for stuff like rounding modes and subnormals. Perhaps the spawned thread's fenv was unexpected.
I'm not super convinced by this blog post anyway, since it immediately confuses associativity with commutativity.
And yea you're right on the associative mistake. Will correct this :)
Floating-point rounding modes seem more likely to me. The author should be able to dump the floating point configuration to confirm, I'm sure.
Confirmed: It’s a bug in MinGW64.
When I create a new thread and run floating point operations in that thread, I get slightly different answers."
Interesting!
Annex F
(normative)
IEC 60559 floating-point arithmetic
F.1 Introduction
This annex specifies C language support for the IEC
60559 floating-point standard [...] previously designated
ANSI/IEEE 754-1985.> The Clang compiler does not support IEC 559 math functionality. Clang also does not control and honor the definition of __STDC_IEC_559__ macro. Under specific options such as -Ofast and -ffast-math, the compiler will enable a range of optimizations that provide faster mathematical operations that may not conform to the IEEE-754 specifications. The macro __STDC_IEC_559__ value may be defined but ignored when these faster optimizations are enabled.
If you make use of Clang, you won't find support for Annex F, and in point of fact you can't even macro check to see if it even will support Annex F.
The story with GCC is... More complicated. It guarantees it'll follow Annex F for only some operations [0].
So the upshot is... Most people can't tell when and how Annex F might be followed.
[0] https://gcc.gnu.org/onlinedocs/gcc-7.4.0/gcc/Floating-point-...
Intel C++ seems to be pretty naughty in that it effectively has fast-math on by default.
(That reminds me, I need to harangue the llvm-libc folks to add that pragma to their fpenv code for correctness sake).
This code was specifically C++ because we had the intention of integrating with libraries like CppAD and ADOL-C.
Examples are: Fortran, R, C#, Java, D and many others.
You'd need to refer to a language's specification to find out, though some aren't very well defined. zig comes to mind; it supports the data types and offers strict mode, but doesn't explicitly state its IEE 754:2008 compliance.
Seriously? Anything but floating point if you care about precision and accuracy.