[Bug libc/3479] New: Incorrect rounding in strtod()

["hack at watson dot ibm dot com" <sourceware-bugzilla at sourceware dot org>]

There is a subtle rounding error in decimal_string->BFP conversion
(e.g. strtod()) for values that are equal to, or very close to, 
half-way points, i.e. sufficiently-long prefixes of the exact decimal
representation of a binary half-way point.
 
An example is strtod(" 3.518437208883201171875E+013 ") which truncates
to the odd 0x42c0000000000001 when in fact default IEEE rounding should
have rounded up to 0x42c0000000000002 (nearest-ties-to-even).
 
In fact, strtod only looks at the first 20 digits (well, it scans the
others for validity), so that strtod(" 3.518437208883201171999E+013 "
also rounds down, even though it is definitely above the half-way point.
Indeed, with an extended fraction this would be
 JL16'+3.518437208883201171999E+013' =  42C00000 00000001 800A66E1 358F83EC
(This is the output of the IBM Research experimental assembler Phantasm,
for which I wrote the conversion routines ten years ago.)
The first number to round up correctly is "3.518437208883201172000E+013".
* (I hope a fixed-width font is used here)  ....:....1....:....2....:....3
 
Here is the broken part of strtod.c -- the comment is simply wrong in
the case of exact half-way numbers, where all digits are significant
in order to determine which way to round!
 
     /* For the fractional part we need not process too many digits.  One
        decimal digits gives us log_2(10) ~ 3.32 bits.  If we now compute
                         ceil(BITS / 3) =: N
        digits we should have enough bits for the result.  The remaining
        decimal digits give us the information that more bits are following.
        This can be used while rounding.  (One added as a safety margin.)  */
     if (dig_no - int_no > (MANT_DIG - bits + 2) / 3 + 1)
       {
         dig_no = int_no + (MANT_DIG - bits + 2) / 3 + 1;
         more_bits = 1;
       }
     else
       more_bits = 0;


 The basic idea of limiting the number of decimal digits that need to
participate in conversion arithmetic is a good one, btw; it's just that
the limit picked by strtod() is too small and does not take the exponent
into account.
 
As many digits as it takes to represent a binary halfway point exactly
must be processed.  The worst case is 752 digits for double, for numbers
near Nmin/2.  If the binary exponent is known at this point (I think it
is, based on quick code browsing), the max can be computed.  I enclose
the comments from my own conversion program:
 
*
*     When the normalised binary exponent B exceeds the frac-
*     tion length F (in bits), the length of an exact decimal
*     represention is .3 B digits, otherwise (e.g. negative B)
*     it is (F + .7 B) digits (or (F-B) + .3 B).  We need to
*     detect exact half-way points, so we need to consider F+1
*     instead of the F derived from the format.
*
 
(Here .3 and .7 are short for log10(2) and log10(5) of course.)
 
Michel.

-- 
           Summary: Incorrect rounding in strtod()
           Product: glibc
           Version: unspecified
            Status: NEW
          Severity: normal
          Priority: P2
         Component: libc
        AssignedTo: drepper at redhat dot com
        ReportedBy: hack at watson dot ibm dot com
                CC: glibc-bugs at sources dot redhat dot com


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

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