Subtle conversion bug in strtod (hence atof)

[1-914-784-7648]

Floating-point conversions with glibc printf() and scanf() are remarkably
good, given that they use gmp high-precision arithmetic internally.

Nevertheless, there is a rounding error in decimalstring->BFP conversion
(e.g. strtod()) for values that are equal to, or very close to, decimal
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;

I found the source above last year on sourceforge, but I verified
the behaviour on recent Linux systems.  The sourceforge id was:

  cvs: linux-vax/glibc/stdlib/strtod.c
  Revision: 1.1.1.1 (vendor branch),
  Fri Feb  9 18:04:23 2001 UTC (4 years, 6 months ago) by airlied
  Branch: MAIN, r2_2_3

Sincerely,
  Michel Hack, IBM Research
Sent: 2006-11-06 21:13:19 UTC