Someone’s Been Messing with My Subnormals
moyix.blogspot.com
moyix.blogspot.com
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.
It seems they changed this behavior in Direct3D 10:
https://microsoft.public.win32.programmer.directx.graphics.n...
I was making a game using D3D, Lua and Chipmunk physics, and some of the behaviour of the game was being odd.
So I started to try printing random stuff with Lua, eventually I just tried: print(5+5), and to my surprise my console outputted "11".
I went into Lua's irc channel to talk about this, and everyone said I was nuts, that the number was too small to trigger precision issues, that I was a troll and so on.
After a lot of searching I found out about this D3D bug, so I switched the game to use OpenGL instead there it was, 5+5 = 10 again!
Now why fiddling with the FPU could make 5+5 become 11, I have no idea.
Many of these same DLLs would also hijack SetUnhandledExceptionFilter() for their custom exception support, which would also result in hard fastfail crashes when they failed to unhook properly. Ended up having to hotpatch SetUnhandledExceptionFilter() Detours-style to prevent my crash reporting filter from being overridden. Years later, Microsoft revealed that Office had done the same thing for the same reasons.
The new version of this problem is DLLs that use AVX instructions and then don't execute a VZEROALL/VZEROUPPER instruction before returning. This is more sinister as it doesn't cause a failure, it just causes SSE2 code to run up to four times slower in the thread.
It is nice that diverse vendor-specific calling conventions and ABIs are less common these days.
This is basically the reason compiler autovectorization doesn't do MMX.
I hope we learned our lessons on this specific question in the design of Wasm. There are subnormals in Wasm and you can't turn them off for performance.
Of course that's why floating point math is mostly a big no-no in driver code. Unless the driver preserves and restores FPU state on its own.
https://github.com/JuliaPackaging/BinaryBuilderBase.jl/blob/...
In user codes you can do `@fastmath`, but it's at the semantic level so it will change `sin` to `sin_fast` but not recurse down into other people's functions, because at that point you're just asking for trouble. There's also calls to rename it `@unsafemath` in Julia, just to make it explicit. In summary, "Fastmath" 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.
Sure, automatic FMA can change the result, but to my knowledge it always gives a more accurate result, not a less accurate one, and the way in which the results may differ is bounded.
For reference, to handle this well we use MuladdMacro.jl which is a semantic transformation that turns x*y+z into muladd expressions, and it does not recurse into functions so it does not change the definitions of the callers inside of the macro scope.
https://github.com/SciML/MuladdMacro.jl
This is something that will always increase performance and accuracy (performance because muladd in Julia is an FMA that is only applied if hardware FMA exists, effectively never resorting to a software FMA emulation) because it's targeted to do only a transformation that has that property.
As a side note, this kind of thing is why I think a good title for a fast-math would be "Fast math, or how I learned to start worrying and hate floating point."
Keep in mind that twiddling these flags is going to require saving the MXCSR register to memory, or'ing or and'ing bits in memory, and then reading that memory back into MXCSR. And both saving and reading the MXCSR requires stalls, because floating point operations both read and write that register. So you require, minimum, 4 L1 cache hits and 2 partial pipeline flushes to twiddle a MXCSR bit.
(As far as I'm aware, modern microarchitectures generally don't register-rename the floating-point status register.)
Note that you wouldn't necessarily need to do a read-modify-write -- it'd suffice in most cases to just to save the old value and then reset the whole MXCSR for the scope requiring special treatment.
And implicit locales should just die. It's sad that even newer functionality like std::format relies on them. Sadder that if it didn't then you'd probably have compiler developers pulling shit like GCC/libstdc++ does for std::from_chars [0] which defeats the entire point of that function but hey, at least they can mark that as implemented. Not like there are suitably licensed implementations of the functionality [1] available that they could use instead if they don't want to implement float parsing themselves.
[0] https://github.com/gcc-mirror/gcc/blob/master/libstdc++-v3/s...
Somewhere on this page there's a discussion of dirty upper half states in Intel vector instructions, and something similar might help there.
For example, if you got "x = (a / b) * (c / d)" one might think that rewriting it as "x = (a * c) / (b * d)" will save you a division and gain you speed. It will and it might, respectively.
However it will also potentially break an otherwise safe operation. If the numbers are very small, but still normal, then the product (b * d) might result in a denormalized number, and dividing by it will result in +/- infinity.
However, the code might guarantee that the ratios (a / b) and (c / d) are not too small or too large, so that multiplying them is guaranteed to lead to a useful result.
Anyone have any references on how the current state of affairs on modern AMD/Intels?
[1]: ARM Cortex-M4 for example can have a hardware FPU, but where division and sqrt are optional, see https://developer.arm.com/documentation/102832/latest/
Looks like one FP divider on modern intel. Though you can pack multiple divisions into an instruction.
For AMD I can find throughput numbers but not how many there are, in a brief search. I'd guess two??
I might be misremembering, but I think fastmath was one of the flags explicitly warned against in the Gentoo manual.
The CPU flags was less interesting to me compared to being able to disable features like X.
Still, there are always (admittedly diminishing) returns if you can target the exact CPU you have - pretty sure no binary distro ever had variants for AMD's TBM instructions [0] but my Gentoo install made use of them (which was made absolutely clear when trying to run any of that on a Zen 2 machine lol).
[0] https://en.wikipedia.org/wiki/X86_Bit_manipulation_instructi...
I've long since stopped worrying about it because on the systems I run, which are not top-of-the-line but aren't RPis either, it's not worth worrying about anymore for most programs. At most maybe you should target the one particular program you use that could use a boost.
It is, here: https://wiki.gentoo.org/wiki/GCC_optimization#But_I_get_bett...
Though compiling for the exact CPU was definitely not insignificant before the amd64 switch which gave a new baseline - many distros still targeted i586.
[1] https://clang.llvm.org/docs/UsersManual.html#cmdoption-ffast...
Actually it wasn't so bad :) All the current debs add up to a mere 113GB, and the .so files are just 46GB once extracted. Quick work. Here's the list:
https://moyix.net/~moyix/debian_unstable_amd64_ffast_results...
And here's a visualization of the reverse dependency graph:
What about -fno-unsafe-math-optimizations?
$ gcc -Ofast -fno-unsafe-math-optimizations -fpic -shared foo.c -o foo.so
$ objdump -j .text --disassemble=set_fast_math foo.so
foo.so: file format elf64-x86-64
Disassembly of section .text:
0000000000001040 <set_fast_math>:
1040: f3 0f 1e fa endbr64
1044: 0f ae 5c 24 fc stmxcsr -0x4(%rsp)
1049: 81 4c 24 fc 40 80 00 orl $0x8040,-0x4(%rsp)
1050: 00
1051: 0f ae 54 24 fc ldmxcsr -0x4(%rsp)
1056: c3 retqhttps://stackoverflow.com/questions/11714827/how-can-i-turn-...
-Weverything is everything but that is only really useful to discover new warning flags (in combination with lots of -Wno...) for which you are better off reading GCC/Clang's release notes and/or documentation.
Do people ever actually test their fixes?
Edit: right, via btrees:
https://github.com/zopefoundation/BTrees/pull/179/files
https://github.com/zopefoundation/meta/commit/2a33089508792e...
Though they use the term denormalized, which AFAICT is the a synonym for subnormal.
Edit: Thanks for this great blogpost :)
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.
Don't you love your fun safe math?
(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.)
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.
https://github.com/gevent/gevent/pull/1820
I haven't examined the code of scipy.stats.skellam.sf so I can't say for sure that it's not converging, but it's clearly some kind of pathological behavior.
Anyhow, my spelunking was cut off by sleep, so the best I can tell that would end up in the CDFLIB[1] routine CUMCHN with X = 8e-6, PNONC = 2e-6, DF from 0 to 99. The insides don’t really look like the kind of magic that is held up by Sterbenz’s lemma and strategically arranged to take advantage of gradual underflow, so at first glance I wouldn’t trust anything subnormal-dependent that it would compute, but maybe it still is? Sleep.
[1] https://people.sc.fsu.edu/~jburkardt/f_src/cdflib/cdflib.f90
I think it suffices to show that the behavior of FTZ/DAZ caused an actual problem for someone, though. I agree that the vast majority of numerical code won't care about FTZ/DAZ, but when it's enabled thread-wide you have no idea what kind of code you'll end up affecting.
I’d have loved if the library author replied with “why don’t you just print out the values directly?”
So it seems likely (if not certain) that the bug report in question says, essentially, that a function inside scipy, given certain arguments, returns one sequence of meaningless bits quickly when gevent is not loaded and a different sequence of meaningless bits very slowly when it is. While that is certainly not nice, and I’d accept it as a bug in gevent, it still does not count as a piece of code succeeding at a reasonable task with subnormals but failing without them.
(I do not have any ideological opposition to such code—it would certainly be interesting to see some!)
It has global effects like those in TFA, and even locally you no longer know if a line or two of arithmetic will become more precise (e.g., by using higher precision intermediate results), less precise, or become complete gibberish (e.g., because it thinks it can prove you're now dividing by zero and thus can just return whatever it wants).
-fgoodenough-math?
-fbroken-but-fast-math
-Ofast
Disregard strict standards compliance. ...
There's strict standards compliance and then there's the crazy grab bag of code changes that is `-ffast-math`. Further, I'd say gevent can defensibly say that -ffast-math is okay for them given what the manual says: -ffast-math
... it can result in incorrect output for programs that depend on an exact
implementation of IEEE or ISO rules/specifications for math functions.
It may, however, yield faster code for programs that do not require the
guarantees of these specifications.
This is 100% on the compiler people. For the option name, the documentation, and the behavior.https://gcc.gnu.org/onlinedocs/gcc-12.2.0/gcc/Optimize-Optio...
That said, I don't see why the -Ofast option even needs to exist, except backwards compatibility, as -ffast-math and the others can (and should IMO) be specified explicitly.
Untrue. The doc entry for -ffast-math says "can result in incorrect output for programs that depend on an exact implementation of IEEE or ISO rules/specifications for math functions". Emphasis mine.
So they clearly say that the entire program can turn invalid when -ffast-math is used.
You and some other people here act like the docs say "translation unit" or something like that, instead of "program", but this is simply not the case.
Furthermore, the entry for -ffast-math points to entries for suboptions that -ffast-math turns on (located right below in the man page), e.g. -funsafe-math-optimization. These also make clear how dangerous they can be even when turned on one at a time.
They: Sir, that is dash F unsafe!
> -cl-unsafe-math-optimizations
> Allow optimizations for floating-point arithmetic that (a) assume that arguments and results are valid, (b) may violate IEEE 754 standard and (c) may violate the OpenCL numerical compliance requirements as defined in section 7.4 for single-precision floating-point, section 9.3.9 for double-precision floating-point, and edge case behavior in section 7.5. This option includes the -cl-no-signed-zeros and -cl-mad-enable options.
While it stops short of saying "this will likely break your code" (maybe because it doesn't have the nonlocal effects of -ffast-math), it makes it much more clear that this flag is generally unsafe and fragile, except under rather specific circumstances. Also, it is reasonably exact about what those circumstances are. I'm not sure -ffast-math is documented with enough precision for a programmer to even know whether it will break their code. Best you can do is try and see if the program still works.
-ffast-math:
> This option is not turned on by any -O option besides -Ofast since it can result in incorrect output for programs that depend on an exact implementation of IEEE or ISO rules/specifications for math functions.
It also point to the -funsafe-math-optimizations sub-option, where it is said that:
> Allow optimizations for floating-point arithmetic that (a) assume that arguments and results are valid and (b) may violate IEEE or ANSI standards. When used at link time, it may include libraries or startup files that change the default FPU control word or other similar optimizations. [...]
It's hard to come up with a similar name that isn't long.
The suggestion given elsewhere in these comments to call it "unsafe math" instead of "fast math" sounds good. It's nearly as short, and properly conveys the "you must know what you're doing" aspect of these flags. It's even better if you're used to Rust.
> Please don't post shallow dismissals, especially of other people's work. A good critical comment teaches us something.
Linking in code with undefined (in this case, redefined) behavior doesn’t automatically invalidate the entire program. But thats the language used because once the undefined behavior is hit at runtime, the spec no longer defines what the behavior is and what the program will do afterwards.
IIRC Swift changed -Ofast to -Ounchecked after I complained about it.
-fcorrupt-quietly
Your own compiled library functions will use the optimizations, but won't force the weird FPU register modes on the rest of the process.
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.
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.
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.
It's -ffast-math isn't it?
...
Yep. That option is a candy coated foot gun.
I believe it’s a necessary price in this case, but it does highlight how suboptimal it is to pay the price in other cases.
(I don't say this because I want to excuse dynamic linking, which I also generally dislike! Only that I think the problem is somewhere else in this particular case.)
Dynamic linking solves a real problem, especially in this space. It comes with new problems of its own but so does the alternative.
But for something like a Python extension it’s what we’ve got.
Which has the ancillary benefit of surfacing stuff like this.