Your program won't print anything, assuming round-ties-to-even rounding mode. The essence of the argument is as follows:
We're assuming that both f
and divisor
are between 1.0
and 2.0
. So f = a / 2^23
and divisor = b / 2^23
for some integers a
and b
in the range [2^23, 2^24)
. The case divisor = 1.0
isn't interesting, so we can further assume that b > 2^23
.
The only way that (float)(f * (1.0 / divisor))
could give the wrong result would be for the exact value f / divisor
to be so close to a halfway case (i.e., a number exactly halfway between two single-precision floats) that the accumulated errors in the expression f * (1.0 / divisor)
push us to the other side of that halfway case from the true value.
But that can't happen. For simplicity, let's first assume that f >= divisor
, so that the exact quotient is in [1.0, 2.0)
. Now any halfway case for single precision in the interval [1.0, 2.0)
has the form c / 2^24
for some odd integer c
with 2^24 < c < 2^25
. The exact value of f / divisor
is a / b
, so the absolute value of the difference f / divisor - c / 2^24
is bounded below by 1 / (2^24 b)
, so is at least 1 / 2^48
(since b < 2^24
). So we're more than 16
double-precision ulps away from any halfway case, and it should be easy to show that the error in the double precision computation can never exceed 16 ulps. (I haven't done the arithmetic, but I'd guess it's easy to show an upper bound of 3 ulps on the error.)
So f / divisor
can't be close enough to a halfway case to create problems. Note that f / divisor
can't be an exact halfway case, either: since c
is odd, c
and 2^24
are relatively prime, so the only way we could have c / 2^24 = a / b
is if b
is a multiple of 2^24
. But b
is in the range (2^23, 2^24)
, so that's not possible.
The case where f < divisor
is similar: the halfway cases then have the form c / 2^25
and the analogous argument shows that abs(f / divisor - c / 2^25)
is greater than 1 / 2^49
, which again gives us a margin of 16
double-precision ulps to play with.