4

I am familiar with the fact that decimal fractions often don't convert neatly into binary, and so an approximation is stored and when converted back into decimal is a bit off from the original number. I believe 0.21 would be an example of such a fraction.

How and why is this program compiled in gcc showing 7.21 accurately?

And (the important second part to this question), why given that it's showing 7.21 accurately, is it showing 839.21 inaccurately?

I can understand why it'd show 0.21 inaccurately, or any number point 21 inaccurately, since it takes a lot of bits to put it in binary accurately, if one even at all. But i'd expect it to consistently show inaccurately show n.21 regardless of what integer n is

printf("%.100f\n", 7.21);  
printf("%.100f\n", 839.21); 
produces 7.210000000000000000000... 
839.2100000000000400000....

Any idea why it does one accurately and one inaccurately? (on my gcc compiler anyway!)

if relevant, gcc --version shows gcc.exe (rubenvb-4.5.4) 4.5.4 Copyright (C) 2010 Free Software Foundation, Inc.

I spoke to one person that said he got both as inaccurate, which is more what one would expect.. So i've taken a screenshow showing the behaviour

enter image description here

Note- I notice that my cygwin gcc implementation is not the latest cygwin gcc implementation.. I just did where gcc it was running it from C:\Perl64\site\bin\gcc.exe so it wasn't running a cygwin one. It probably was running an old ming one(Which is what R suggested). chux's is cygwin n GNU C11 (GCC) version 6.4.0 (i686-pc-cygwin)`.

barlop
  • 12,887
  • 8
  • 80
  • 109
  • 1
    Try `volatile double x = 7.21; printf("%.100f\n", x);` to see if optimizations are in play. – chux - Reinstate Monica Nov 28 '17 at 22:44
  • @chux volatile double also prints 7.21 exactly – barlop Nov 28 '17 at 22:59
  • 1
    Thaks for checking that. I am confident that the issue is a weak `printf()` with FP arguments. If possible with your compiler, use `printf("%a\n", x);` to see a hexadecimal significand. – chux - Reinstate Monica Nov 28 '17 at 23:07
  • @chux `0x1.cd70a4p+2` what does a weak printf() mean? – barlop Nov 28 '17 at 23:21
  • barlop, Is `printf("%.100f\n", 7.21);` the true code or are you printing a `float`? It is possible that even your `printf("%a\n"...` is weak too. `0x1.cd70a4p+2` is expected with `volatile float x = 7.21; printf("%a\n", x);` – chux - Reinstate Monica Nov 28 '17 at 23:24
  • A "weak" `printf()` mean that the code is a bit sloppy/inaccurate with the lesser significant digits. – chux - Reinstate Monica Nov 28 '17 at 23:27
  • @chux just literals and doubles. And the code I am running is what I have said I have run . So if you asked me to run it with a double variable I ran it with a double variable and if my screenshot shows a line where it runs with a literal then that is what ran, and if I said I ran it with a literal then I did. And as I said, I get the same result whether the variable is volatile or not, and whether it's a literal or a variable. – barlop Nov 28 '17 at 23:31
  • Then concerning output `0x1.cd70a4p+2` aka (7.2100000381...), code is internally converting to `float` and then printing. `0x1.cd70a3d70a3d7p+2` is expected for an accurate `double` print. – chux - Reinstate Monica Nov 28 '17 at 23:35
  • Let us [continue this discussion in chat](http://chat.stackoverflow.com/rooms/160036/discussion-between-barlop-and-chux). – barlop Nov 28 '17 at 23:42
  • Sounds like you're using mingw not cygwin. This is a known "feature" of Microsoft's printf implementation... – R.. GitHub STOP HELPING ICE Nov 29 '17 at 05:01

2 Answers2

4

How and why is this program compiled in gcc showing 7.21 accurately?

The older printf() handling on the least significant decimal digits is limited.

printf("%.100f\n", x); need not print x accurately passed a certain digit. There really is no C spec on this yet decimal output should be at least correct for DBL_DIG digits (typically 15) - this would be "weak" and problematic.

A better, and often acceptable, printf() will be correct for at least DBL_DECIMAL_DIG digits (typically 17). Without going to unlimited precision, getting the last digit correct can be troublesome. See Table-maker's dilemma. Padding the "right" with zeros rather than the correct digits is not uncommon. This is what OP's code did. It went to the correct rounded 17 digits and then zero padded.

A high quality printf() will print double x = 839.21; and x = 7.21; correctly for all digits. E.g.: 

839.2100000000000363797880709171295166015625...
839.210000000000000019984014443252817727625370025634765625.... (as long double)
vs OP's
839.2100000000000400000....

123.4567890123456789 // digit place

7.20999999999999996447286321199499070644378662109375...
7.210000000000000000034694469519536141888238489627838134765625000... (as long double)
7.210000000000000000000

OP's printf() is only good up to 16 or so digits.

The output of 7.210000000000000000000.... only looks great because the printf() went out to a point and then padded with zeros. See @Eric Postpischil and luck

Note: Some optimizations will perform the task with long double (research FLT_EVAL_METHOD) so long double as x86 extended precision results posted too.

Added by Barlop

Some additional points from Eric Postpischil

OP's printf is not exact for 7.21. OP's printf was passed exactly 7.20999999999999996447286321199499070644378662109375, but it printed 7.210…, so it is wrong. It just happens to be wrong in a way that the error in that the OP's printf exactly cancelled out the error that occurred when 7.21 was rounded to binary floating-point.

Eric's correct, printf("%.100f\n", 7.20999999999999996447286321199499070644378662109375); prints 7.210000000

Eric elaborates on how he knows it is 7.20999999999999996447286321199499070644378662109375 which was sent to printf and not some other long number close to that.

Eric commented that he knows that, by using a good C implementation that converts decimal numerals in source code to binary floating-point properly and that prints properly. (Apple’s developer tools on macOS.) And used some Maple (math software) code he wrote to help him with floating-point arithmetic. In the absence of those, he might have had to work it out the long way, just like in elementary school but with more digits.

barlop
  • 12,887
  • 8
  • 80
  • 109
chux - Reinstate Monica
  • 143,097
  • 13
  • 135
  • 256
  • You say my printf is only good up to 16 or so digits, but look at the result I get for 7.21 it's exact, i.e. 7.2100000000000000000000000000000000000000... – barlop Nov 28 '17 at 22:57
  • @Clifford Answer amended. OP's `printf()` is only printing correctly to about 15 decimal digits and then next few are about right. After the 17th or so digit, the functions simply pads with `'0'`. The fact that it printed the exact `7.210000...` is luck. – chux - Reinstate Monica Nov 28 '17 at 23:02
  • @barlop Try more than just 2 examples. Try difficult values like `double x = 1.0 - DBL_EPSILON;` or `1.0 + DBL_EPSILON, DBL_MIN, DBL_MAX, DBL_TRUE_MIN`. – chux - Reinstate Monica Nov 28 '17 at 23:04
  • 4
    @barlop: Your `printf` is not exact for 7.21. Your `printf` was passed exactly 7.20999999999999996447286321199499070644378662109375, but it printed 7.210…, so it is wrong. It just happens to be wrong in a way that the error in your `printf` exactly canceled out the error that occurred when 7.21 was rounded to binary floating-point. – Eric Postpischil Nov 28 '17 at 23:07
  • @EricPostpischil Oops I originally posted `long double` (10 byte) values. 8 byte [binary64](https://en.wikipedia.org/wiki/Double-precision_floating-point_format) ones now included. – chux - Reinstate Monica Nov 28 '17 at 23:13
  • @EricPostpischil you're right `printf("%.100f\n", 7.20999999999999996447286321199499070644378662109375);` prints 7.210000000 – barlop Nov 28 '17 at 23:27
  • @EricPostpischil How can you test that it is `7.20999999999999996447286321199499070644378662109375` which was sent to printf and not some other long number close to that? – barlop Nov 28 '17 at 23:58
  • @barlop: In this case, I used a good C implementation that converts decimal numerals in source code to binary floating-point properly and that prints properly. (Apple’s developer tools on macOS.) I also have some Maple (math software) code I wrote to help me with floating-point arithmetic. In the absence of those, I might have to work it out the long way, just like in elementary school but with more digits. – Eric Postpischil Nov 29 '17 at 00:02
  • @EricPostpischil ok.. I think your comment (with 3 upvotes so far), has hit the nail on the head in answering the question, you're welcome to post it as an answer., with any elaboration eg the elaboration of your recent comment. – barlop Nov 29 '17 at 00:21
  • @barlop: Thanks, but chux can have the points. They can include text from my comment if they wish. This answer includes additional information I think is useful and/or interesting, although perhaps it could be given some additional context to help explain it. – Eric Postpischil Nov 29 '17 at 00:24
  • @EricPostpischil it might be that your implementation is the latest gnu gcc which may give the same result as chux's. Chux is using 6.4.0(albeit cygwin's). Granted you're using a mac osx / BSD derivative implementation, Which version are you using? Do you really think my implementation is not a good one? What Arithmetic is yours doing that is different to the arithmetic that mine is doing? – barlop Nov 29 '17 at 00:34
  • @barlop: On my personal system, I currently have Apple LLVM 8.1.0, based off clang-802.0.42. Apple switched from GCC not too long ago. If your implementation is not printing 7.20999999999999996447286321199499070644378662109375 exactly, then it could be better. There is a notable paper about this, *Correctly Rounded Binary-Decimal and Decimal-Binary Conversions* by David M. Gay, November 30, 1990. I do not know what the current state of the art is. – Eric Postpischil Nov 29 '17 at 00:53
  • Let us [continue this discussion in chat](http://chat.stackoverflow.com/rooms/160040/discussion-between-barlop-and-eric-postpischil). – barlop Nov 29 '17 at 01:16
0

I will keep the accepted answer as accepted.

But I will ad this answer to address the question and the other part of the question about the 839.21

There appears to be a bug in printf of that cygwin gcc implementation. Though R thought it was a ming thing, and there is a case a similar printf bug with mingw which deeper, is a bug with the Windows C/C++ runtime files. Though Chux has found that his cygwin implementation is fine, his cygwin implementation is more current.

I have found an online C compiler which doesn't have the bug so it illustrates things well. https://www.tutorialspoint.com/compile_c_online.php

One of the things my question asked was why 7.21 was displaying accurately and 839.21 was displaying inaccurately, given that they have the same fraction after the decimal point.

printf("%.100f\n", 7.21);
printf("%.100f\n", 839.21);

It's worth looking at a printf that doesn't have the bug.

printf("%f\n",7.21); // 7.210000

printf("%f\n",839.21); // 839.210000

We can see how they are stored in double

printf("%.60f\n",7.21);

7.209999999999999964472863211994990706443786621093750000000000


printf("%.20a\n",7.21);  

0x1.cd70a3d70a3d70000000p+2


printf("%.60f\n",839.21);

839.210000000000036379788070917129516601562500000000000000000000


printf("%.20a\n",839.21);

0x1.a39ae147ae1480000000p+9

Notice that the fraction that is stored in binary, i.e. the number after the binary point, is very different, is very different for 7.21 than for 839.21 This is because in scientific notation, there is only one digit before the point. So even if it were just decimal numbers in scientific notation, the fraction part would be different. 7.21 and 8.3921 If you look at 839 in binary it is 1101000111 So you can see the 1.A there, as the first bit is 1, then for the following nibble you have 1010 which is A. Different to the case of 7.21

So that's why one shouldn't expect the same result in terms of accuracy, storing 7.21 as for 839.21

Chux has mentioned there is an issue in my implementation involving an internal conversion to float in the %a case. So it's worth looking at how they are stored in float, by doing float f; and passing to printf with a format specifier of %.60f and %.20a still using the mentioned working implementation of the printf.

printf("%.60f\n",f); //float f=7.21  
7.210000038146972656250000000000000000000000000000000000000000

printf("%.20a\n",f); //float f=7.21  
0x1.cd70a400000000000000p+2 <-- 7.21 in float, in hex.

printf("%.60f\n",f); //float f=839.21  
839.210021972656250000000000000000000000000000000000000000000000
0x1.a39ae200000000000000p+9 <-- 839.21  in float, in hex

Now if we look at my question, Chux noticed that there is a bug that it converts the number to float in the %a case.

So for example in my early(not latest) cygwin implementation, which is what my question asks about and where a bug is, double x=7.21; printf("%a\n", x); was printing the hex but it was showing the hex of a number stored as float. 0x1.cd70a4p+2

There may be more to the bug, because the %f case may be a bit different, but anyhow, definitely a printf bug, in the early cygwin implementation gcc.exe (rubenvb-4.5.4) 4.5.4 Copyright (C) 2010 Free Software Foundation, Inc. (Whether that affects an early mingw one too I don't know).

There is a related printf bug here with mingw's printf, which is down to the C/C++ runtime files Long double is printed incorrectly with iostreams on MinGW But my issue is resolved in the latest cygwin gcc that chux uses. (I will try the latest cygwin to double check that).

barlop
  • 12,887
  • 8
  • 80
  • 109