Re: Precision of doubles and stdio
Jim Easton wrote:
Dear Mr. Bagnara,
Roberto Bagnara wrote: ...
does this on Linux/i686
$ a.out 70.9 70.900000000000005684341886080801486968994140625
and does the following under Cygwin on the same machine:
...
$ ./a.exe 70.9 70.90000000000000568434188608080148696899414
Why? Is there a way to reconcile the two behaviors? Notice that I know about the x87 and its vaguaries: nonetheless I wonder why such a scanf immediately followed by a printf shows a difference between Cygwin and Linux.
With all due respect, why would you want to? With double you are guaranteed only 16 or so digits - the rest is noise. Frankly I am amazed that it agrees as far as it does.
Dear Jim,
what is and what is not noise depends on the application. In our applications we systematically use controlled rounding on IEEE 754 floating point numbers. In the end, what we obtain (in memory) are definite (i.e., provably correct) lower or upper bounds of some quantities.
Call `x' such a quantity, and suppose we have that our computed upper bound for `x' is the IEEE 754 Double Precision number
0x4051b9999999999a,
that is (exactly!),
70.900000000000005684341886080801486968994140625.
If that number is (wrongly!) printed as
70.90000000000000568434188608080148696899414
then we lose correctness, since x <= 70.900000000000005684341886080801486968994140625 does not imply x <= 70.90000000000000568434188608080148696899414. So, the final "0625" is not "noise" in our applications: it is what may make the difference between a correct statement and an incorrect one.
Notice also that any IEEE 754 number can be (exactly!) printed with a finite number of decimal digits. Finally, notice that writing an algorithm to print them correctly is not rocket science. Hence my astonishment when Tim showed me that the C standard decided, instead, to allow blatant violations of the principle of least astonishment :-) All the best,
Roberto
[ X-Posted and Followups set to newlib list; this is almost certainly not a cygwin specific issue. To recap:- ]
On 03 March 2006 21:44, Roberto Bagnara wrote:
Hi there,
the following little program
#include <stdio.h>
int main() { double d; scanf("%lf", &d); printf("%.1000g\n", d); return 0; }
does this on Linux/i686
$ gcc -W -Wall in.c $ a.out 70.9 70.900000000000005684341886080801486968994140625
and does the following under Cygwin on the same machine:
roberto@quark /tmp $ gcc -W -Wall in.c
roberto@quark /tmp $ ./a.exe 70.9 70.90000000000000568434188608080148696899414
On 03 March 2006 22:21, Tim Prince wrote:
If you haven't gone out of your way to install similar printf() support libraries on cygwin and linux, they will definitely not be the same. My past reading of various relevant documents convinced me that digits beyond the 17th in formatting of doubles are not required by any standard to be consistent between implementations. They have no useful function, as 17 digits are sufficient to determine uniquely the corresponding binary value in IEEE 754 format.
Jim Easton wrote:
With all due respect, why would you want to? With double you are guaranteed only 16 or so digits - the rest is noise. Frankly I am amazed that it agrees as far as it does.
On 06 March 2006 09:30, Roberto Bagnara wrote:
In our applications we systematically use controlled rounding on IEEE 754 floating point numbers. In the end, what we obtain (in memory) are definite (i.e., provably correct) lower or upper bounds of some quantities.
Call `x' such a quantity, and suppose we have that our computed upper bound for `x' is the IEEE 754 Double Precision number
0x4051b9999999999a,that is (exactly!),
70.900000000000005684341886080801486968994140625.If that number is (wrongly!) printed as
70.90000000000000568434188608080148696899414then we lose correctness, since x <= 70.900000000000005684341886080801486968994140625 does not imply x <= 70.90000000000000568434188608080148696899414. So, the final "0625" is not "noise" in our applications: it is what may make the difference between a correct statement and an incorrect one.
Notice also that any IEEE 754 number can be (exactly!) printed with a finite number of decimal digits. Finally, notice that writing an algorithm to print them correctly is not rocket science. Hence my astonishment when Tim showed me that the C standard decided, instead, to allow blatant violations of the principle of least astonishment :-)
I agree with you that the number is of course precisely representable. Even if Tim isn't onto a red herring (which I think he probably is) and the standard really does give us leeway in this case, we could still do better anyway. Newlib has had a few glitches and corner cases in its fp support before and probably will do again; we can look into this. I notice, for instance, that newlib doesn't define DECIMAL_DIG or the related symbols AFAICT.
Are the newlib-based cygwin system and the glibc-based linux system using the same mode bits in the fpu control register? That '0625' looks awfully like a 1-lsb error, maybe they are using different rounding or guard bits or something like that?
cheers, DaveK
On 06 March 2006 13:33, Dave Korn wrote:
On 03 March 2006 21:44, Roberto Bagnara wrote:
does this on Linux/i686
$ gcc -W -Wall in.c $ a.out 70.9 70.900000000000005684341886080801486968994140625
and does the following under Cygwin on the same machine:
roberto@quark /tmp $ gcc -W -Wall in.c
roberto@quark /tmp $ ./a.exe 70.9 70.90000000000000568434188608080148696899414
On 03 March 2006 22:21, Tim Prince wrote:
If you haven't gone out of your way to install similar printf() support libraries on cygwin and linux, they will definitely not be the same. My past reading of various relevant documents convinced me that digits beyond the 17th in formatting of doubles are not required by any standard to be consistent between implementations. They have no useful function, as 17 digits are sufficient to determine uniquely the corresponding binary value in IEEE 754 format.
I don't see why we still shouldn't print out well-defined numbers precisely.
Sure, if you /assume/ that the number in question was originally converted from an ASCII representation, there's no point converting it back with any more precision in a way that would falsely suggest it was possible to discriminate between numbers that have the same representation once scanned, but if we /don't/ assume the fp number in question was obtained by converting ascii -> float, then there's no reason to arbitrarily truncate the displayed value at some point just because it would convert back to the same fp value if we re-scanf'ed it. We should assume the fp number could just as easily be a precise representation of the result of a calculation.
Finally, I would suggest that in any case, if the user has explicitly requested such a long field length, there's even less reason to truncate the output early.
Notice also that any IEEE 754 number can be (exactly!) printed with a finite number of decimal digits. Finally, notice that writing an algorithm to print them correctly is not rocket science.
I agree with you that the number is of course precisely representable.
before and probably will do again; we can look into this. I notice, for instance, that newlib doesn't define DECIMAL_DIG or the related symbols AFAICT.
Well, it may not have DECIMAL_DIG, but it has an internal symbol called NDEC defined in .../src/newlib/libc/stdlib/ldtoa.c that seems to serve more or less the same purpose:
-----------------------------------snip----------------------------------- /* Number of 16 bit words in external x type format */ #define NE 10
/* Number of 16 bit words in internal format */ #define NI (NE+3)
/* Array offset to exponent */ #define E 1
/* Array offset to high guard word */ #define M 2
/* Number of bits of precision */ #define NBITS ((NI-4)*16)
/* Maximum number of decimal digits in ASCII conversion * = NBITS*log10(2) */ #define NDEC (NBITS*8/27)
/* The exponent of 1.0 */ #define EXONE (0x3fff)
/* Maximum exponent digits - base 10 */ #define MAX_EXP_DIGITS 5 -----------------------------------snip-----------------------------------
Roberto, just as an experiment I turned NDEC up to 64, and I now get the exact same output on Cygwin that you get from your linux version of the code. I can't see any problems that should arise from this, but I don't understand the reasoning behind the NBITS*log10(2) formula, or to be precise, I don't understand why that is or isn't enough. I suspect it has something to do with the fact that the number of ascii digits needed to represent the conversion of NBITS bits from binary to decimal isn't the same thing as the number of digits required to represent the reciprocal of (an NBITS-sized number), but that's not worded very precisely.... what I'm getting at is that b100.0 can be represented in 1 decimal sig.fig., but b0.001 takes 3 decimal sig.figs; the calculation is valid for NBITS binary places above the decimal point, but not NBITS binary places below it.
cheers, DaveK
participants (2)
-
Dave Korn -
Roberto Bagnara