> IIRC dividing a subnormal by a subnormal returns infinity.
This should not be true in a conformant implementation of IEEE 754-2008. You can get infinity by dividing by a subnormal (due to overflow), but that should only be possible with sufficiently large numerators.
> The problem was that in certain circumstances the terms were all small, but of comparable magnitude. Thus the latter lead to subnormal divided by subnormal. The fix was to undo the optimization and do the divisions first, which lead to the multiplication of three numbers of order 1.
My guess is that this problem was actually the numerator being subnormal and the denominator underflowing all the way to zero and therefore the result going to infinity.
It's been over a decade, so I struggle to recall the details.
> My guess is that this problem was actually the numerator being subnormal and the denominator underflowing all the way to zero and therefore the result going to infinity.
This should not be true in a conformant implementation of IEEE 754-2008. You can get infinity by dividing by a subnormal (due to overflow), but that should only be possible with sufficiently large numerators.
> The problem was that in certain circumstances the terms were all small, but of comparable magnitude. Thus the latter lead to subnormal divided by subnormal. The fix was to undo the optimization and do the divisions first, which lead to the multiplication of three numbers of order 1.
My guess is that this problem was actually the numerator being subnormal and the denominator underflowing all the way to zero and therefore the result going to infinity.