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.