9

With float a = ...; and float inva = 1/a; is x / a the same as x * inva?

And what is with this case:

unsigned i = ...;
float v1 = static_cast<float>(i) / 4294967295.0f;
float scl = 1.0f / 4294967295.0f;
float v2 = static_cast<float>(i) * scl;

Is v1 equal to v2 for all unsigned integers?

Danvil
  • 22,240
  • 19
  • 65
  • 88
  • 3
    The first answer here might provide some info for you: http://stackoverflow.com/questions/22621241/what-does-the-constant-0-0039215689-represent – C.B. Mar 26 '14 at 14:37
  • The same question keeps coming up in different forms. See http://floating-point.gui.de/ – devnull Mar 26 '14 at 14:41
  • 3
    @devnull: that site is dead. – Violet Giraffe Mar 26 '14 at 14:42
  • 2
    Not the point here, but `4294967295.0f` is not even representable as a float. – yzt Mar 26 '14 at 14:48
  • @VioletGiraffe ok, this shouldn't be dead: http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html – devnull Mar 26 '14 at 14:50
  • Note: there was a discussion on Reddit recently where the author pointed out that for questions of numerical accuracy in `float`, since they were only 32 bits (4*10^9 combinations) one could easily and exhaustively test all combinations for single operand operations. In your case you have two operands so an exhaustive check might be a bit too much, however you could easily test random floats and check that both methods yield the same result, in a couple of minutes you would know. – Matthieu M. Mar 26 '14 at 14:51
  • I've checked here trying all possible unsigned ints and v1 was always equal to v2. – Jabberwocky Mar 26 '14 at 14:57
  • @yzt That's interesting, because it means `4294967295.0f` “is” a power of two, and therefore `x / 4294967295.0f == x * (1 / 4294967295.0f)` for all `x`. – Pascal Cuoq Mar 26 '14 at 15:09
  • 1
    @MatthieuM. That's for one float. Two floats, as in `x / y == x * (1 / y)`, make just fewer than 2^64 possibilities in general, unless you can argue that most exponents need not be studied, in which case you may end up with a barely practical set of 2^46 possibilities. – Pascal Cuoq Mar 26 '14 at 15:11
  • 1
    @Michael Walz I too am unable to a failing case for a 32-bit `unsigned`. Even tried `unsigned` as a 64-bit - so far no failures. (this is with OP's 2nd case) – chux - Reinstate Monica Mar 26 '14 at 15:13
  • @chux It is normal not to find a counter-example for the division by `4294967295.0f`. – Pascal Cuoq Mar 26 '14 at 15:16
  • @Pascal Cuoq Tis true if `float` is the typical [binary32](http://en.wikipedia.org/wiki/Binary32) but not if `float` is binary64 (do not know a machine where that is true). I wonder about rounding modes? – chux - Reinstate Monica Mar 26 '14 at 15:24
  • @chux In C, the conversion of inexact decimal floating-point literals produces an implementation-defined result that has to be either the floating-point number above or that below. In practice, most compilers round to nearest-even (barring implementation bugs where difficult cases are rounded in the wrong direction, of which examples are on Rick Regan's blog). – Pascal Cuoq Mar 26 '14 at 15:28
  • @Pascal Cuoq If the rounding mode was round toward 0, OP's 2nd case would have many examples where `v1 != v2`. Agree that the usual round-to-nearest would not normally provide a counter example. – chux - Reinstate Monica Mar 26 '14 at 15:38
  • @devnull - Should be called "What Every Computer Scientist Should Know About Floating-Point Arithmetic" but does not. Interestingly, in the answer below it is stated that odd number is a power of two. – SChepurin Mar 26 '14 at 16:26
  • @yzt I did not say anything about 4294967295. I am talking about `4294967295.0f`. On a typical compilation platform with single-precision `float` and `FLT_EVAL_METHOD` 0. Dude. – Pascal Cuoq Mar 27 '14 at 13:45

2 Answers2

7

is v1 equal to v2 for all unsigned integers?

Yes, because 4294967295.0f is a power of two. Division and multiplication by the reciprocal are equivalent when the divisor is a power of two (assuming the computation of the reciprocal does not overflow or underflow to zero).

Division and multiplication by the reciprocal are not equivalent in general, only in the particular case of powers of two. The reason is that for (almost all) powers of two y, the computation of 1 / y is exact, so that x * (1 / y) only rounds once, just like x / y only rounds once.

Pascal Cuoq
  • 79,187
  • 7
  • 161
  • 281
  • 4
    @Danvil The funny thing is that 4294967295 is not a power of two, but `4294967295.0f`, on most compilation platforms, is the single-precision IEEE 754 number nearest 4294967295, that is, 4294967296, and is in effect a power of two. – Pascal Cuoq Mar 26 '14 at 15:30
  • 3
    Note that 4294967295.0f may in fact be considered `==` to 4294967296.0f, since there are loss-of-precision issues trying to represent 2^32-1 as a float... So, while it isn't truly a power of 2 (no odd number can be, by definition, except for 2^0 = 1), it may appear to be... – twalberg Mar 26 '14 at 15:32
5

No, the result will not always be the same. The way you group the operands in floating point multiplication, or division in this case, has an effect on the numerical accuracy of the answer. Thus, the product a*(1/b) might differ from a/b. Check the wikipedia article http://en.wikipedia.org/wiki/Floating_point#Accuracy_problems.

Rontogiannis Aristofanis
  • 8,883
  • 8
  • 41
  • 58