Someone’s Been Messing with My Subnormals

At a previous company I worked at, we had an issue with our software (Windows-based, written in a proprietary language) randomly crashing. After some debugging, we found that this happened whenever the user made some specific actions, but only if, in that session, the user had previously printed something or opened a file picker. The culprit was either a printer driver or a shell extension which, when loaded, changed the floating point control word to trap. That happened whenever the culprit DLL had been compiled by a specific compiler, which had the offending code in the startup routine it linked into every DLL it produced.<p>Our solution was the inverse of the one presented in this article: instead of wrapping our routines to temporarily set the floating point control word to sane values, we wrapped the calls to either printing or the file picker, and reset the floating point control word to its previous (and sane) value after these calls.

The Julia package ecosystem has a lot of safeguards against silent incorrect behavior like this. For example, if you try to add a package binary build which would use fast math flags, it will throw an error and tell you to repent:<p><a href="https:&#x2F;&#x2F;github.com&#x2F;JuliaPackaging&#x2F;BinaryBuilderBase.jl&#x2F;blob&#x2F;186b1350fa925ff9722d320b4c8549ec6a5735db&#x2F;src&#x2F;Runner.jl#L289" rel="nofollow">https:&#x2F;&#x2F;github.com&#x2F;JuliaPackaging&#x2F;BinaryBuilderBase.jl&#x2F;blob&#x2F;...</a><p>In user codes you can do `@fastmath`, but it&#x27;s at the semantic level so it will change `sin` to `sin_fast` but not recurse down into other people&#x27;s functions, because at that point you&#x27;re just asking for trouble. There&#x27;s also calls to rename it `@unsafemath` in Julia, just to make it explicit. In summary, &quot;Fastmath&quot; is overused and many times people actually want other optimizations (automatic FMA), and people really need to stop throwing global changes around willy-nilly, and programming languages need to force people to avoid such global issues both semantically and within its package ecosystems norms.

The problem here is that enabling FTZ&#x2F;DAZ flags involves modifying global (technically thread-local) state that is relatively expensive to do. Ideally, you&#x27;d want to twiddle these flags only for code that wants to work in this mode, but given the relative expense of this operation, it&#x27;s not entirely practicable to auto-add twiddling to every function call, and doing it manually is somewhat challenging because compilers tend to support accessing the floating-point status rather poorly. Also, FTZ&#x2F;DAZ aren&#x27;t IEEE 754, so there&#x27;s no portable function for twiddling these bits as there is for other rounding mode or exception controls. I will note that icc&#x27;s -fp-model=fast and MSVC&#x27;s &#x2F;fp:fast correctly do not link code with crtfastmath.<p>As a side note, this kind of thing is why I think a good title for a fast-math would be &quot;Fast math, or how I learned to start worrying and hate floating point.&quot;

Global state is the root of so many evils! FPU rounding mode, FPU flush-to-zero mode, C locale, errno, and probably some other things should all be eliminated. The functionality should still exist but not as global flags.

Denormalized numbers is one reason why you really want to think carefully if you try to optimize code by rewriting expressions involving multiplication and division.<p>For example, if you got &quot;x = (a &#x2F; b) * (c &#x2F; d)&quot; one might think that rewriting it as &quot;x = (a * c) &#x2F; (b * d)&quot; will save you a division and gain you speed. It will and it might, respectively.<p>However it will also potentially break an otherwise safe operation. If the numbers are <i>very</i> small, but still normal, then the product (b * d) might result in a denormalized number, and dividing by it will result in +&#x2F;- infinity.<p>However, the code might guarantee that the ratios (a &#x2F; b) and (c &#x2F; d) are not too small or too large, so that multiplying them is guaranteed to lead to a useful result.

I ran Gentoo back in the good old days. The biggest draw was that after about a week of compiling my system ran a lot faster because of all the compiler optimisations one could enable because it only had to work on your CPU.<p>I might be misremembering, but I think fastmath was one of the flags explicitly warned against in the Gentoo manual.

-funsafe-math is neither fun nor safe.

See also <a href="https:&#x2F;&#x2F;simonbyrne.github.io&#x2F;notes&#x2F;fastmath&#x2F;" rel="nofollow">https:&#x2F;&#x2F;simonbyrne.github.io&#x2F;notes&#x2F;fastmath&#x2F;</a> for a similar story in julia, where ffast-math is now banned for C&#x2F;C++&#x2F;Fortran dependencies

Wow, I am surprised that -ffast-math triggers a mode switch in the FPU in part due to the author&#x27;s library problem, but also because the documentation for clang at least[1] does not say it impacts behaviour of denormals and in fact has a separate mode switch for that, which is not explicitly called out as being implied by -ffast-math.<p>[1] <a href="https:&#x2F;&#x2F;clang.llvm.org&#x2F;docs&#x2F;UsersManual.html#cmdoption-ffast-math" rel="nofollow">https:&#x2F;&#x2F;clang.llvm.org&#x2F;docs&#x2F;UsersManual.html#cmdoption-ffast...</a>

Does this only affect pypi, or should I now worry about shared libraries shipped with my distro as well? Debian is not crazy enough to ship shared libs compiled with -ffast-math, right? RIGHT?

&gt; it turns out that when you use -Ofast, -fno-fast-math does not, in fact, disable fast math. lol. lmao.<p>What about -fno-unsafe-math-optimizations?

For anyone (like me) not knowing what subnormal floats are, this StackOfverflow question and answer explain it quite nicely: <a href="https:&#x2F;&#x2F;stackoverflow.com&#x2F;questions&#x2F;15140847&#x2F;denormalized-numbers-ieee-754-floating-point#" rel="nofollow">https:&#x2F;&#x2F;stackoverflow.com&#x2F;questions&#x2F;15140847&#x2F;denormalized-nu...</a><p>Though they use the term denormalized, which AFAICT is the a synonym for subnormal.<p>Edit: Thanks for this great blogpost :)

&gt; Finally, I am legally obligated to inform you that I used GNU Parallel to scan all the wheels<p>FWIW, this is tangential to this awesome article, but if the author is here or anyone else who cares: that statement isn’t true, you are not legally obligated to mention GNU Parallel. It is nice to do though! The link even says this explicitly, and also separately mentions the citation notice is only asking for scientific paper citations, not blog citations. I love parallel, but boy has that citation notice thing caused a lot of confusion over the years.

Oh man, great job digging through all that. This is exactly the kind of content I want to see.<p>Don&#x27;t you love your fun safe math?

This is a rockstar quality post. It is astonishing how much detective work was involved.

That’s…terrifying. This is a fantastic find: big, big respect to @moyix, this is going to save people’s ass.

10&#x2F;10 yak shave. Would certainly read again

Following the article’s links, I fail to find an actual example of anything failing to converge in flush-subnormals mode. I mean, I’m sure one could be squeezed out, but the justification given amounts to “Sterbenz’s lemma [the one that rephrases “catastrophic cancellation” as “exact differences”] fails, maybe something somewhere also will”. And my (shallow but not nonexistent) experience with numerical analysis is that proofs lump subnormals with underflow, and most of them don’t survive even intermediate underflows.<p>(AFAIU the original Intel justification for pushing subnormals into 754 was gradual underflow, i.e. to give people at least something to look at for debugging when they’ve ran out of precision.)<p>So, yes, it’s not exactly polite to fiddle with floating-point flag bits that are not yours, and it’s better that this not happen for reproducibility if nothing else, but I doubt it actually breaks any interesting numerics.

-Ofast isn&#x27;t a good name for the option, but in GCC&#x27;s defense the manual is pretty clear about all this, and there&#x27;s no excuse for blindly turning on compiler options - they literally change the semantics of your code.

If you are willing to accept the various caveats that come with -ffast-math for your own library, then it looks like it&#x27;s okay to use -ffast-math during <i>compilation</i>, but not during <i>linking</i> (because it links in crtfastmath.so, as the author points out).<p>Your own compiled library functions will use the optimizations, but won&#x27;t force the weird FPU register modes on the rest of the process.

That thread-local MXCSR register is particularly entertaining in a thread pool environment, such as OpenMP. OSes carefully preserve that piece of thread state across context switches.<p>I tend to avoid touching that value, even when it means extra instructions like roundpd for specific rounding mode, or shuffles to avoid division by 0 in the unused lanes.

In D you can opt into specific float algorithms locally rather than a compiler flag.<p>Use of fast math can really really really bite you sometimes, so just being able to opt into using fma and nothing else is awesome.

A decorator is a nice idea for this.<p>I was going to suggest another package that just resets the MXCSR when imported, but I guess... hypothetically... some function might actually want the FTZ behavior.

&gt; built with an appealing-sounding but dangerous compiler option<p>It&#x27;s -ffast-math isn&#x27;t it?<p>...<p>Yep. That option is a candy coated foot gun.

I thought the purpose of Python was to make development simple and predictable. Needing to track down the compilation and linker flags of every single shared library reveals the fallacy of this abstraction.

Dynamic linking is the root of all kinds of evil, enough said.