7

Is division by zero possible in the following case due to the floating point error in the subtraction?

float x, y, z;
...
if (y != 1.0)
    z = x / (y - 1.0);

In other words, is the following any safer?

float divisor = y - 1.0;
if (divisor != 0.0)
    z = x / divisor;
R. Martinho Fernandes
  • 228,013
  • 71
  • 433
  • 510
Max
  • 1,203
  • 1
  • 14
  • 20
  • For the second code you posted - do you mean: `if (y != 1.0)`? – jrd1 Oct 10 '12 at 18:01
  • It both results in the same (I assume you meant divisor != 0.0f, different might be with something like y = 13.0f-12.0f; (y-1.0f) could work – SinisterMJ Oct 10 '12 at 18:01
  • I believe these are the same. And yes, *never* check equality with floating-point numbers: `y != 1.0` is off-limits. – ev-br Oct 10 '12 at 18:02
  • Oops, I meant 'if (divisor != 0.0)'. I'll update. Thanks! – Max Oct 10 '12 at 18:03
  • If your using the lazy initialization the first one is best, else there is no difference. – andre Oct 10 '12 at 18:04
  • 3
    Why are you checking for exactly `0.0` (or `1.0`)? If `y - 1.0` is small enough or `x` is big enough you will get `inf` even for a `divisor != 0`, so it doesn't help prevent that and it would probably be easier to check for `inf` after the division then working out when the division will give you infinity. So even if both codepaths are identical, it doesn't make any difference for the safety. – Grizzly Oct 10 '12 at 18:07
  • 2
    Note that `1.0` isnt a `float` -- it's a `double` – John Dibling Oct 10 '12 at 18:15
  • 1
    @Grizzly: That's *if* you can be sure that division by zero will yield infinity (or NaN for 0.0/0.0). The language doesn't guarantee that behavior. – Keith Thompson Oct 10 '12 at 18:23

3 Answers3

6

Assuming IEEE-754 floating-point, they are equivalent.

It is a basic theorem of FP arithmetic that for finite x and y, x - y == 0 if and only if x == y, assuming gradual underflow.

If subnormal results are flushed to zero (instead of gradual underflow), this theorem holds only if the result x - y is normal. Because 1.0 is well scaled, y - 1.0 is never subnormal, and so y - 1.0 is zero if and only if y is exactly 1.0, regardless of how underflow are handled.

C++ doesn't guarantee IEEE-754, of course, but the theorem is true for most "reasonable" floating-point systems.

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

This will prevent you from dividing by exactly zero, however that does not mean still won't end up with +/-inf as a result. The denominator could still be small enough so that the answer is not representable with a double and you will end up with an inf. For example:

#include <iostream>
#include <limits>

int main(int argc, char const *argv[])
{
    double small = std::numeric_limits<double>::epsilon();
    double large = std::numeric_limits<double>::max() / small;
    std::cout << "small: " << small << std::endl;
    std::cout << "large: " << large << std::endl;
    return 0;
}

In this program small is non-zero, but it is so small that large exceeds the range of double and is inf.

David Brown
  • 13,336
  • 4
  • 38
  • 55
  • 2
    `y - 1` can't ever be `denorm_min()`, the smallest it can be is `-epsilon()/2`, and it won't be infinity. – kennytm Oct 10 '12 at 18:27
  • @KennyTM That is a good point. The answer could still be `inf` if the numerator is large enough though. I'll update my answer. – David Brown Oct 10 '12 at 18:31
2

There is no difference between the two code snippets () - in fact, the optimizer could even optimize both fragments to the same binary code, assuming that there are no further uses of the divisor variable.

Note, however, that division by a floating point zero 0.0 does not result in a run-time error, but produces an inf or -inf instead.

Sergey Kalinichenko
  • 714,442
  • 84
  • 1,110
  • 1,523
  • 2
    I wouldn't bet on compilers doing any apparently "obvious" transformations on floating point calculations, at least not without being asked to with special command line flags. Floating point arithmetic is very tricky. For example, it's [not associative and GCC respects that](http://stackoverflow.com/q/6430448/395760). –  Oct 10 '12 at 18:18