A raytracer I worked on had an interesting bug, where a manual optimization caused some unforeseen floating-point issues down the line.
The short story was that we had an expression of the form x = a/b * c/d * e/f in a fairly hot code path. As most developers know, divisions are expensive, so this got rewritten to x = (a * c * e) / (b * d * f).
The problem was that in certain edge cases, the factors were all extremely small. This caused the product b * d to become a denormal number, and when multiplied with f it became zero. This of course caused the division to return infinity.
However, closer inspection of the source of the terms a to f revealed that the pairs a,b etc always had roughly the same magnitude. Thus even though a and b were almost denormal numbers, the division a / b was well-behaved and gave a sensible result, and similar for the others.
Undoing the optimization and using the original form made the code well-behaved again.
Another expression that had to be de-optimized for the same reason was x = a / (b * b), which had to be rewritten as x = (a / b) / b.
So as with subtraction, division of IEEE floating-point numbers can be tricky.
This[1] page contains a nice set of tables with the results of performing IEEE floating-point operations on edge-case inputs.
The short story was that we had an expression of the form x = a/b * c/d * e/f in a fairly hot code path. As most developers know, divisions are expensive, so this got rewritten to x = (a * c * e) / (b * d * f).
The problem was that in certain edge cases, the factors were all extremely small. This caused the product b * d to become a denormal number, and when multiplied with f it became zero. This of course caused the division to return infinity.
However, closer inspection of the source of the terms a to f revealed that the pairs a,b etc always had roughly the same magnitude. Thus even though a and b were almost denormal numbers, the division a / b was well-behaved and gave a sensible result, and similar for the others.
Undoing the optimization and using the original form made the code well-behaved again.
Another expression that had to be de-optimized for the same reason was x = a / (b * b), which had to be rewritten as x = (a / b) / b.
So as with subtraction, division of IEEE floating-point numbers can be tricky.
This[1] page contains a nice set of tables with the results of performing IEEE floating-point operations on edge-case inputs.
[1]: https://www.cs.uaf.edu/2011/fall/cs301/lecture/11_09_weird_f...