7

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

Aziz Shaikh
  • 16,245
  • 11
  • 62
  • 79
Jeremy Salwen
  • 8,061
  • 5
  • 50
  • 73
  • 1
    How small are you talking? You might have some issues with `0.0000000000000000000000003`, so to speak, but it's hard to answer the question without knowing what precision your input floats will be. – Corey Jun 02 '10 at 02:08
  • AFAIK, you could get an arithmetic overflow with small divisors, but not a division by zero error. – Dave Markle Jun 02 '10 at 02:10
  • 2
    How about `if ( A*B*B < 1 )`? – Nikolai Fetissov Jun 02 '10 at 02:11
  • The B input is coming from a user graphical slider or textual input, so I should be able to handle any size gracefully, even if some mischievous person puts in 0.0000000000000000000000003. The A input is from a floating point subtraction calculation of numbers greater than ten or so, so I would also expect it to possibly be very small, but ridiculously small, because it should have lost accuracy. – Jeremy Salwen Jun 02 '10 at 02:19
  • Nikolai, I think you might have the solution. That should never result in any overflow or division by zero of any kind, and is exactly the sort of answer I'm looking for. You should have posted it as an answer! – Jeremy Salwen Jun 02 '10 at 02:23
  • Just to be clear, `0.0000000000000000000000003` isn't even close to being in the ballpark where this would start to be a concern (at least for double precision). In fact, you can't even type a number that *would* cause an issue in a stackoverflow comment field (unless you allow scientific notation). – Stephen Canon Jun 02 '10 at 03:36

4 Answers4

7

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.

Stephen Canon
  • 103,815
  • 19
  • 183
  • 269
5

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.

Martin Thompson
  • 16,395
  • 1
  • 38
  • 56
Nikolai Fetissov
  • 82,306
  • 11
  • 110
  • 171
  • 3
    Actually, as long as he doesn't need to support zero or infinity as values for `a` or `b` (it sounds like he doesn't), then there are no corner case issues; they all fall out correctly. – Stephen Canon Jun 02 '10 at 03:13
3

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.

Ben Robbins
  • 2,889
  • 2
  • 31
  • 32
  • 1
    Nice! Measure twice, cut once. To asker: if this code is not self-explanatory, then add necessary comments to make this clear. – Hamish Grubijan Jun 02 '10 at 02:37
  • Division is somewhat slow, but not toxic. The problem here is that b*b can easily overflow. a*b where a>1 and b<1 is a far better place to start. – Potatoswatter Jun 02 '10 at 03:00
  • @Potatoswatter: Assuming IEEE-754 floating point, the overflow is harmless, and this will give the correct answer anyway. Ditto for underflow. – Stephen Canon Jun 02 '10 at 03:12
  • @stephen: Not quite… if `a` is denormal then `b*b` may overflow to infinity and throw off the result. – Potatoswatter Jun 02 '10 at 03:18
  • @Potatoswatter: ah, didn't catch that this one was `b*b*a` instead of the correct `a*b*b` (which *does* work properly, because multiplication associates left-to-right in C). – Stephen Canon Jun 02 '10 at 03:25
1

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.

Potatoswatter
  • 134,909
  • 25
  • 265
  • 421