[Bug libc/4943] Inconsistent rounding behaviour for sprintf and IEEE doubles

["paul at inet dot co dot za" <sourceware-bugzilla@sourceware.org>]

------- Additional Comments From paul at inet dot co dot za  2007-09-22 20:19 -------
(In reply to comment #24)

> And again, I suggest that you read
>   http://www2.hursley.ibm.com/decimal/
> (in particular the FAQ).

I have read this document, thanks.  FP is accurate enough for my purposes in all
calculations.  Once the result needs to be displayed though, I have an entirely
different situation.  Yes, of course I have work-arounds, as shown in my
original example.  Why are these work-arounds necessary though?

> Paul, you really make every one (including me) think that you do not 
> understand floating-point arithmetic. 

This is the strangest thing.  I probably understand it better than most, which
is really scary when you think about it.  Since the late 80's I could quote the
IEEE storage format without a second thought.

Perhaps it is not common to have people query these things.  Everyone just seems
to accept that "floating point is inaccurate" and are happy to leave it at that.
 When someone comes along and rocks the boat people get all upset.

> Double-precision numbers are binary numbers, not decimal numbers. Saying 
> "15-digit precision" is misleading; 

I thought I was using the correct terminology when I said "15 decimals of
precision" or phrased differently, "15 decimal digits of precision"...

> you should say "53-bit precision" since their representation is in binary.

well, not quite 53 bits, rather 52.  The extra bit is implied simply because the
FPU does not store the value 0.  When you break it down, though, it's not even
52 bits, as you cannot ignore calculation and storage error.  The largest
mantissa is (decimal) 9007199254740991 which is 16 digits.  We know that there
is always an error associated with storage, and therefore we cannot consider
more than 15 decimals (decimal digits) to be accurate.

> Now, if you have so much confidence in Microsoft's C library, I suggest that 
> you try the following program:

Erg.  Confidence?  I think my words were "Microsoft does not seem to have this
problem"...

> #include <stdio.h>
> int main (void)
> {
>   volatile double x = 2597.525;
>   volatile double y = 5000.525;
>   printf ("%.2f\n%.2f\n", x, y);
>   return 0;
> }
> 
> (please do not change the numbers), and see if you get "consistent" results for
> x and for y.

Ok, Vincent, you are right.  Well done.  In this case, even my "work-around" fails.

MS stores the value
  5000.525
as
  5000.524999999999600

Ignoring the error, the value 
  5000.52499999999 
being less than
  5000.525
rounds to
  5000.52
(closest)

Linux produces the identical result.

Should the FPU be ignoring the error when promoting 5000.5249999999996 to
5000.52499999999 to account for precision?  Or should it not instead have
promoted the value to 5000.52500000000?  I cannot answer that. Looks like a true
case of FP inaccuracy to me.

This is very different from the case of 2597.625 where we had an exactly stored
value being rounded towards the even value simply because it was exactly stored
and on a boundary.  

Seems no-one can get this right.  I will see what Microsoft have to say about
this.  They don't have the "round to nearest even on boundaries" argument to
throw my way.  I will certainly inform you as to the outcome.

Their result is certainly inconsistent with your earlier example:

#include <stdio.h>
int main (void)
{
  double x;
  for (x = -8.5; x <= 8.5; x += 1.0)
    printf ("%4g %.0f\n", x, x);
  return 0;
}

[Microsoft output]

-8.5 -9
-7.5 -8
-6.5 -7
-5.5 -6
-4.5 -5
-3.5 -4
-2.5 -3
-1.5 -2
-0.5 -1
 0.5 1
 1.5 2
 2.5 3
 3.5 4
 4.5 5
 5.5 6
 6.5 7
 7.5 8
 8.5 9


-- 


http://sourceware.org/bugzilla/show_bug.cgi?id=4943

------- You are receiving this mail because: -------
You are on the CC list for the bug, or are watching someone who is.