Bug in floating-point addition: is anyone else seeing this?

[Henrique Dante de Almeida]

On May 22, 1:41 am, Henrique Dante de Almeida <hda... at gmail.com>
wrote:
>
> > >  Notice that 1e16-1 doesn't exist in IEEE double precision:
> > >  1e16-2 == 0x1.1c37937e07fffp+53
> > >  1e16 == 0x1.1c37937e08p+53
>
> > >  (that is, the hex representation ends with "7fff", then goes to
> > > "8000").
>
> > >  So, it's just rounding. It could go up, to 1e16, or down, to 1e16-2.
> > > This is not a bug, it's a feature.
>
> >  I didn't answer your question. :-/
>
> >  Adding a small number to 1e16-2 should be rounded to nearest (1e16-2)
> > by default. So that's strange.
>
> >  The following code compiled with gcc 4.2 (without optimization) gives
> > the same result:
>
> > #include <stdio.h>
>
> > int main (void)
> > {
> >         double a;
>
> >         while(1) {
> >                 scanf("%lg", &a);
> >                 printf("%a\n", a);
> >                 printf("%a\n", a + 0.999);
> >                 printf("%a\n", a + 0.9999);
> >         }
>
> > }
>
>  However, compiling it with "-mfpmath=sse -msse2" it works. (it
> doesn't work with -msse either).

 Finally (and the answer is obvious). 387 breaks the standards and
doesn't use IEEE double precision when requested to do so.

 It reads the 64-bit double and converts it to a 80-bit long double.
In this case, 1e16-2 + 0.9999 == 1e16-1. When requested by the printf
call, this 80-bit number (1e16-1) is converted to a double, which
happens to be 1e16.