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

Question

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

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, because it is visually unambiguous which multiplication occurs first.

Answer 1

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.

Answer 2

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 745 for your corner cases.

Answer 3

Granted, this does not exactly answer the question directly, but anybody doing floating point math should read this:
What Every Computer Scientist Should Know About Floating-Point Arithmetic, by David Goldberg, published in the March, 1991 issue of Computing Surveys. Copyright 1991, Association for Computing Machinery, Inc.

Answer 4

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.