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
– Steve Summit, Commented Jul 20, 2021 at 16:33
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
– Steve Summit, Commented Jul 20, 2021 at 16:41
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". — Steve Summit
– Steve Summit, Commented Jul 21, 2021 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 github.com/bminor/glibc/blob/master/sysdeps/ieee754/dbl-64/… so as to not cause overflows and maintain accuracy — Sam Mason
– Sam Mason, Commented Jul 21, 2021 at 11:33

chux · Accepted Answer · 2021-07-20 17:13:51Z

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 · Accepted Answer · 2021-07-20 18:26:08Z

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.

Collectives™ on Stack Overflow

Explanation for floating point rounding effect

2 Answers 2

Comments

Comments

Your Answer

Linked

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

Comments

Comments

Your Answer

Sign up or log in

Post as a guest

Linked

Related