This is the mail archive of the
`glibc-bugs@sourceware.org`
mailing list for the glibc project.

Index Nav: | [Date Index] [Subject Index] [Author Index] [Thread Index] | |
---|---|---|

Message Nav: | [Date Prev] [Date Next] | [Thread Prev] [Thread Next] |

Other format: | [Raw text] |

*From*: "hack at watson dot ibm dot com" <sourceware-bugzilla at sourceware dot org>*To*: glibc-bugs at sources dot redhat dot com*Date*: 7 Nov 2006 17:20:25 -0000*Subject*: [Bug libc/3479] New: Incorrect rounding in strtod()*Reply-to*: 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.

Index Nav: | [Date Index] [Subject Index] [Author Index] [Thread Index] | |
---|---|---|

Message Nav: | [Date Prev] [Date Next] | [Thread Prev] [Thread Next] |