Yes this is technically true, but it manifests itself through use of the MinGW compiler. When I first noticed the problem, there was very little indication as to the cause.

TBH I accredit dumb luck more than anything to finding the cause of this. It took many hours.

Op Here. Yes you are correct. Been very sleep derived when I wrote this (newborns eh). I'll update this.

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.

The processor itself supports configuring some aspects of floating point behavior; on x86 this is the MXCSR register. The most important one is rounding mode: you can choose round-to-nearest, round-down, round-up, or round-towards-zero. Less common options include whether to raise an exception on various out-of-range conditions, and whether to treat denormal numbers as zero (faster but less accurate and not IEEE 754 compliant).

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

At a minimum, when there is a switch between threads the floating point state of the old thread needs to be saved and the floating point state of the incoming thread loaded. This includes such things as rounding mode.

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)

Things like rounding mode and denornal behavior could be controlled by CPU flags. Imy not familiar with x86, but this is defiy the case in PPC.

Some kind of internal state specific to Microsoft's libm? It doesn't look like _fpreset() is present anywhere else.

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.

TBH. I had no idea this was a thing either. I was incredibly confused when I identified a single piece of code that was producing an variation in result.

Probably the author is accidentally getting 80-bit "extended double" arithmetic.

> Support for Annex F (IEEE-754 / IEC 559) of C99/C11

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

C "supports" whatever the hardware and compiler feel like supporting. It absolutely does not mandate anything. In particular, there are a number of particularly simple compiler optimizations that are not forbidden, though they are not technically correct according to IEEE 754, such as algebraic reassociation and simple commutativity. Moreover, C allows subexpressions to be computed in higher precision (e.g. 80 bit "long double"), which is observable. That last one is primarily due to the x87 FPU coprocessor design that has given us a good 35 years of headaches. Good riddance to that!

The clang docs explicitly say that they don't support Annex F at all. Is this story finally going to change? So that you have at least partial support?

No, they consistently perform optimizations that violate IEEE 754 rules, such as reassociation and commutation of expressions. These are difficult to observe (unless you are in the habit of observing -0 vs 0), but yes, you can observe them. They are very precisely specified by IEEE 754 and C/C++ compilers perform illegal optimizations.

The correct - but useless - answer is: any language that's IEC 60559:xxxx conformant.

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.

COBOL

Seriously? Anything but floating point if you care about precision and accuracy.