I don't think there is any bug here. This is what I've seen from a
little digging:
1) cygwin strtod rounds to even, with about DECIMAL_DIG (==21) digits
precision, as recommended by 7.20.1.3 of WG14/N843. (It acts strange
when the rounding mode is not round to nearest, but since newlib
doesn't provide fesetround(), that's your own problem if you insist on
changing rounding mode by hand, as I understand it).
2) cygwin gcc rounds with *lots* more precision than DECIMAL_DIG,
about 50 digits. It rounds towards zero. These behaviours are a little
different than the recommendation in 6.4.4.2 Ibid., which says
constants should be converted as if by strtod() at runtime, but its
not a big difference.
3) cygwin printf correctly rounds to even, even out to about 41 digits.
4) I have no idea what mingw is doing, but it's different to the
above. Gcc constructs the same double precision constants as on cygwin
but strtod() is different and seems to have less precision, and
printf() seems to work with about 16 digits precision. At a glance I'd
say it rounds to zero, with exactly the precision of
On 05/07/05, Dave Korn wrote:
> Stepping through the code at the weekend, I followed the 0.105 case as far
> down as ldtoa_r, where I observed that although there was code that would
> round 0.105 up to 0.11 if asked for only two sig.figs, the loop that
> generates successive digits was coming up with the sequence
> '0.10499999999.....', and because the rounding only looks at one digit
> beyond the requested s.f., it 'correctly' rounds 0.104 down to 0.10.
The closest double to 0.105 is actually
0.10499999999999999611421941381195210851728916168212890625, and this
is what printf() is (correctly) handed, and what it (correctly) rounds
to 0.10.
> I say 'correctly' in quotes, because it's the correct way to round 0.104,
> but it's not in general correct to round a number by truncating it and then
> rounding the truncated version.
But I don't think that "truncating it and then rounding the truncated
version" is what printf() does. If we give it 0.125 and 0.375, which
are both representible exactly, then it rounds to 0.12 and 0.38, ie
one goes up and one goes down, to round to even. You can see it's not
simply truncating if you give it 0.1250000000000001, which gets
rounded to 0.13 not 0.12 as it would for truncation.
To summarize: cygwin does the right thing
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <float.h>
#ifdef __MINGW32__
#include <fenv.h>
#endif
// show the hex and decimal interpretations of the double
union U{double d;unsigned char c[8];};
void display(union U x)
{
for(int b=7;b>=0;b--)printf("%02x",x.c[b]);
printf(" %.21f",x.d);
}
int main(int argc, const char**argv)
{
#ifdef __MINGW32__
if (argc >= 2)
{
int nrm, n;
n = sscanf (argv[1], "%i", &nrm);
if (n == 1) switch (nrm)
{
case 0:
case 0x400:
case 0x800:
case 0xc00:
fesetround (nrm);
}
}
printf ("Rmode $%08x\n", fegetround ());
#endif
printf("DECIMAL_DIG=%d\n",DECIMAL_DIG);
#define PASTE(x) {.s=#x,.d1=strtod(#x,0),.d2=x}
struct S {const char*s;union U d1,d2;}q[]=
{
// test the rounding with the last digits...
// first right on the boundary, to test round-to-even
PASTE(1.0000000000000000000000000000000000000000000000000000000), // 1
PASTE(1.0000000000000001110223024625156540423631668090820312500), // 1 +
2^-53
PASTE(1.0000000000000002220446049250313080847263336181640625000), // 1 +
2^-52
PASTE(1.0000000000000003330669073875469621270895004272460937500), // 1 +
2^-52 +2^-53
{0},
// now just above the boundary
PASTE(1.00000000000000000001),
PASTE(1.00000000000000011103),
PASTE(1.00000000000000022205),
PASTE(1.00000000000000033307),
{0},
// now just below the boundary
PASTE(0.99999999999999999999),
PASTE(1.00000000000000011101),
PASTE(1.00000000000000022204),
PASTE(1.00000000000000033301),
{0},
PASTE(0.105),
};
printf(" 1.23456789012345678901234567890\n"); // for counting
sig figs
for(int a=0;a<sizeof(q)/sizeof(struct S);a++)
{
if(!q[a].s){printf("\n");continue;}
display(q[a].d1);
printf(" ");
display(q[a].d2);
printf("\n");
}
printf("\n");
union U w[]={0.105,0.115,0.125,0.135,0.145,0.155,0.375};
for(int a=0;a<sizeof(w)/sizeof(union U);a++)
{
display(w[a]);
printf("\n");
}
printf("%.40f\n",0.105);
// here's some numbers with an exact representation in double
printf("%.20f %.2f\n",.125,.125);
printf("%.20f %.2f\n",.375,.375);
printf("%.20f %.13f\n",0.95367431640625000000,0.95367431640625000000);
printf("%.20f %.17f\n",3.5096855163574218750,3.5096855163574218750); // 1 +
657899/2^18
printf("%.20f %.15f\n",-1.5035324096679687500,-1.5035324096679687500);
// how many digits precision does gcc use when making double constants?
display((union U){.d=1.000000000000000333066907387546962127089500427246});
}
// Useful Mathematica stuff...
/*
Table[{x/1000.,
BaseForm[FromDigits[RealDigits[#2, 2] // First // Rest, 2], 16],
N[#2*2^(-53 + #[[2]]),
41]} &[#,
Ceiling[FromDigits[#[[1]], 2]/
2]] &[RealDigits[x/1000, 2, 54]], {x, 105, 155, 10}]
*/
--
Unsubscribe info: http://cygwin.com/ml/#unsubscribe-simple
Problem reports: http://cygwin.com/problems.html
Documentation: http://cygwin.com/docs.html
FAQ: http://cygwin.com/faq/
