Explanation for floating point rounding effect

Question

#include <stdio.h>
int main() {
    printf("%.14f\n", 0.0001f * 10000000.0f);  // 1
    printf("%.14f\n", 0.001f * 1000000.0f);  // 2
    printf("%.14f\n", 0.01f * 100000.0f);  // 3
    return 0;
}

Godbolt

The output of this code is:

1000.00000000000000
1000.00006103515625
1000.00000000000000

I know, decimal fractions cannot be represented exactly with floats. But why are lines 1 and 3 calculated correctly and 2 is not? Do you have a clear explanation of what is going on here in detail?

Eric will be along with a better answer shortly, so I'll just say: sometimes you get round-up in one place and round-down in the other, and they cancel. Also, in general you wouldn't necessarily say that lines 1 and 3 are calculated *correctly*, just that the remaining error is smaller than `%.14f` will show. (But, yes, for a conventional `double`, `%.14f` should show you just about all there is.) — Steve Summit, Jul 20 '21 at 16:33
It's a little like asking, "I'm terrible at darts. I'm usually in the outer circle, if I don't hit the wall outside the dartboard. But yesterday I got two bullseyes out of three! How could that happen?" — Steve Summit, Jul 20 '21 at 16:38
Expanding on my first comment, the nearest float to 0.01 is `0.0099999997`, the nearest float to 0.001 is `0.001000000047`, and the nearest float to 0.0001 is `0.000099999997`. So there's a bit of a pattern. — Steve Summit, Jul 20 '21 at 16:41
After the multiplications, two of them would come out to 999.99997, except that's not representable as a `float`; the nearest value is 1000 exactly. But the other one comes out to 1000.000047, which isn't representable as a `float` either, but the nearest value is 1000.000061. That is, 1000.000061 is closer to 1000.000047 than 1000 is. — Steve Summit, Jul 20 '21 at 16:48
@SteveSummit If an outcoming result is not representable as a float, how can an FPU handle this result and how is the rounding done? — strange-corner, Jul 21 '21 at 07:03
@strangecorner In most floating-point schemes, the rule is that when a result cannot be represented exactly, the nearest representable number should be used instead. That's why 1000.000047 turns into 1000.000061. — Steve Summit, Jul 21 '21 at 11:22
The people who design floating-point algorithms spend a lot of time trying to meet this requirement, because it's not always easy. Generally you have to work with a few more (sometimes a lot more) bits internally, to hold your intermediate results, so that the final result can be rounded correctly. The final result should be correct to within half a "[unit in the last place](https://en.wikipedia.org/wiki/Unit_in_the_last_place)". — Steve Summit, Jul 21 '21 at 11:27
as an example, `hypot(x,y)` is defined to just do `sqrt(x*x + y*y)`, but is implemented in glibc as the much more complicated https://github.com/bminor/glibc/blob/master/sysdeps/ieee754/dbl-64/e_hypot.c so as to not cause overflows and maintain accuracy — Sam Mason, Jul 21 '21 at 11:33
I think you know this, but: the end result for you or me does not look like "rounding". When 1000.000047 goes to 1000.000061, what kind of rounding is that? It looks weird because, internally, computer floating-point numbers aren't base 10, they're base 2. The `float` numbers in the vicinity of 1000 are, in decimal: 999.999877, 999.999938, 1000.000000, 1000.000061, 1000.000122. The fractions look ridiculous, but they're nice and "even" in binary. In hex they're `0x3e7.fff8`, `0x3e7.fffc`, `0x3e8.0000`, `0x3e8.0004`, `0x3e8.0008`. — Steve Summit, Jul 21 '21 at 11:44

chux - Reinstate Monica · Accepted Answer · 2021-07-20T17:13:51.480

Sometimes the cumulative roundings (3 steps each in OP samples) result in the same as the mathematical/decimal one, sometimes not. @Steve Summit, @Steve Summit

a clear explanation of what is going on here in detail?

There are 3 steps of potential rounding for each line of code:

Source code to float. Recall common float are of the form: some_limited_integer * 2^some_power.
float multiplication with a rounding
Printing of a float rounded to 14 decimal places^*1.

For printf("%.14f\n", 0.0001f * 10000000.0f); // 1

Code 0.0001f to a float with a value of 0.0000999999974737875163555145263671875
0.0000999999974737875163555145263671875 * 10000000.0 --> 999.999974737875163555145263671875 --> rounded to nearest float --> 1000.0
1000.0 --> "1000.00000000000000".

For printf("%.14f\n", 0.001f * 1000000.0f); // 2

Code 0.001f to a float with a value of 0.001000000047497451305389404296875
0.001000000047497451305389404296875 * 1000000.0 --> 1000.000047497451305389404296875 --> rounded to nearest float --> 1000.00006103515625
1000.00006103515625 --> "1000.00006103515625".

In #1 the rounding was down, then up - tending to cancel.
In #2, the rounding was up and up - resulting in noticeable double rounding effect.

Roughly, each step may inject up to a 1/2 ULP error.

Other considerations: 1) alternative rounding mode. Above uses round to nearest. 2) Weak libraries. Above assumes a quality printf().

^*1 In OP's samples, there was no rounding error. In general, printing float with "%f" can round.

Steve Summit · Answer 2 · 2021-07-20T18:26:08.163

The other way of answering this is to say that you did not get one "wrong" answer and two "right" answers. You actually got three "right" answers, where "right" means, "as good as could have been expected".

Type float only gives you about 7 decimal digits of precision. So for numbers in the range of 1000, that's three places past the decimal. So change the program like this:

printf("%.3f\n", 0.0001f * 10000000.0f);
printf("%.3f\n", 0.001f * 1000000.0f);
printf("%.3f\n", 0.01f * 100000.0f);

The output is:

1000.000
1000.000
1000.000

No discrepancy, all answers apparently correct.

Or, do it in exponential notation, with one digit before the decimal and 6 after.

printf("%.6e\n", 0.0001f * 10000000.0f);
printf("%.6e\n", 0.001f * 1000000.0f);
printf("%.6e\n", 0.01f * 100000.0f);

gives

1.000000e+03
1.000000e+03
1.000000e+03

Again, all answers the same.

This might seem like "cheating": we happen to know there's some "interesting" things going on in the digits off to the right of what we can see, so isn't it wrong to suppress them in this way, and make it look like all the answers were the same? I'm not going to answer that question, other than to point out that when you're doing floating-point work, there's almost always something off there to the right, that you're rounding off and suppressing -- it's just a matter of degree.

Explanation for floating point rounding effect

2 Answers2

Linked