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.
(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.
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.
So somebody tried to calculate, for integer arguments from 0 to 99 inclusive, the CDF of the difference of two Poisson variables with means 4e-6 and 1e-6? I... don’t know if it is at all reasonable to expect an answer to that question. As in, genuinely don’t know—obviously it’s an utterly rotten thing to compute, but at the same time maybe somebody got really interested in that and figured out a way to make it work.
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.
Yeah, unfortunately I have no idea if that was their original goal (which seems unlikely?) or if this is just a minimal example they came up with after tripping over the actual problem in a more realistic setting.
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.
My last bug report I wrote a small C++ program to put all the values between 0x000 .. 0xfff into a tree structure and then iterate over the tree printing out the values.
I’d have loved if the library author replied with “why don’t you just print out the values directly?”
My point wasn’t that the program in the bug report did not do anything meaningful (duh). My point was that most numerical algorithms, once they hit intermediate underflow in the form of subnormals for a set of inputs, do not proceed to compute a result that is in any way related to their stated purpose (with the notable exception of string conversions). Unless a hard limit on the number of iterations is imposed, most iterative algorithms can also be made to hang by giving them unfortunate inputs. In both cases it is the responsibility of the numerical analyst invoking them to not provide such inputs.
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!)
(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.