If I want to check that positive float A is less than the inverse square of another positive float B (in C99), could something go wrong if B is very small?
I could imagine checking it like
if(A<1/(B*B))
But if B is small enough, would this possibly result in infinity? If that were to happen, would the code still work correctly in all situations?
In a similar vein, I might do
if(1/A>B*B)
... which might be slightly better because B*B might be zero if B is small (is this true?)
Finally, a solution that I can't imagine being wrong is
if(sqrt(1/A)>B)
which I don't think would ever result in zero division, but still might be problematic if A is close to zero.
So basically, my questions are:
EDIT: for those of you who are wondering, I ended up doing
if(B*A*B<1)
I did it in that order as it is visually unambiguous which multiplication occurs first.
If you want to handle the entire range of possible values of A
and B
, then you need to be a little bit careful, but this really isn't too complicated.
The suggestion of using a*b*b < 1.
is a good one; if b
is so tiny that a*b*b
underflows to zero, then a
is necessarily smaller than 1./(b*b)
. Conversely, if b
is so large that a*b*b
overflows to infinity, then the condition will (correctly) not be satisfied. (Potatoswatter correctly points out in a comment on another post that this does not work properly if you write it b*b*a
, because b*b
might overflow to infinity even when the condition should be true, if a
happens to be denormal. However, in C, multiplication associates left-to-right, so that is not an issue if you write it a*b*b
and your platform adheres to a reasonable numerics model.)
Because you know a priori that a
and b
are both positive numbers, there is no way for a*b*b
to generate a NaN, so you needn't worry about that condition. Overflow and underflow are the only possible misbehaviors, and we have accounted for them already. If you needed to support the case where a
or b
might be zero or infinity, then you would need to be somewhat more careful.
To answer your direct questions: (answers assume IEEE-754 arithmetic)
Can 1/X ever be infinity if X is greater than zero (but small)?
Yes! If x is a small positive denormal value, then 1/x
can overflow and produce infinity. For example, in double precision in the default rounding mode, 1 / 0x1.0p-1024
will overflow.
Can X*X ever be zero if X is greater than zero?
Yes! In double precision in the default rounding mode, all values of x smaller than 0x1.0p-538
(thats 2**-578
in the C99 hex format) or so have this property.
Will comparisons with infinity work the way I would expect them to?
Yes! This is one of the best features of IEEE-754.
OK, reposting as an answer.
Try using arithmetically equivalent comparison like if ( A*B*B < 1. )
. You might get in trouble with really big numbers though.
Take a careful look at the IEEE 754 for your corner cases.
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With