Can bad stuff happen when dividing 1/a very small float?

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 p开发者_如何学编程ossibly 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:

• Can 1/X ever be infinity if X is greater than zero (but small)?
• Can X*X ever be zero if X is greater than zero?
• Will comparisons with infinity work the way I would expect them to?

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.

You want to avoid divisions so the trick is to modify the equation. You can multiply both sides of your first equation by (b*b) to get:

b*b*a < 1.0

This won't have any divisions so should be ok.

Division per se isn't so bad. However, standard IEEE 754 FP types allow for a greater negative negative range of exponents than positive, due to denormalized numbers. For example, float ranges from 1.4×10^-45 to 3.4×10^-38, so you cannot take the inverse of 2×10^-44.

Therefore, as Jeremy suggests, start by multiplying A by B, where one has a positive exponent and the other has a negative exponent, to avoid overflow.

This is why A*B*B<1 is the proper answer.