Incorrect Directed Conversions in David Gay’s strtod()

For correctly rounded decimal to floating-point conversions, many open source projects rely on David Gay’s strtod() function. In the default rounding mode, IEEE 754 round-to-nearest, this function is known to give correct results (notwithstanding recent bugs, which have been fixed). However, in the less frequently used IEEE 754 directed rounding modes — round toward positive infinity, round toward negative infinity, and round toward zero — strtod() gives incorrectly rounded results for some inputs.

(Update 7/11/10: dtoa.c is now fixed — see below.)

Discovering the Problem

As a follow up to my article “Visual C++ and GLIBC strtod() Ignore Rounding Mode,” I asked Mark Dickinson — the developer who maintains Python’s copy of dtoa.c — why Python does not enable what I call directed conversions. Out of that discussion came Mark’s conjecture: the directed rounding code was not correct. He soon followed up with an example to prove it.

An Automated Way to Find Examples

I wanted to find more examples, but I wanted a C program to do it for me. I needed a way to generate examples whose correct outcomes were known (in my evaluations of round-to-nearest conversions, David Gay’s strtod() served this purpose). Inspired by Mark’s example, this is the algorithm I devised:

1. Generate a random decimal string that represents a number that is not an integer and is not dyadic.
2. Convert the decimal string to a double (I used David Gay’s strtod(), but the built-in strtod() would have been fine — correct rounding does not matter here).
3. Convert the double to a decimal string, giving the exact decimal representation of a floating-point number (I used David Gay’s g_fmt() function, since I’m running on Windows and printf() in Visual C++ can’t print all the decimal digits of a floating-point number).
4. Change the last digit of the new decimal string — which by the way is always 5 (this article indirectly explains why) — to either 4 or 6, depending on whether you want a number just below or just above the floating-point number.

Except for Mark’s example, the examples below were found by my C program implementation of this algorithm. I ran the program against a copy of dtoa.c most recently changed on “Tue Feb 2 23:05:34 MST 2010,” which I downloaded on May 27, 2010.

Examples: Decimal Numbers Just Above a Floating-Point Number

This section shows two examples that David Gay’s strtod() (hereafter referred to only as strtod()) rounds incorrectly; they are very close to, but just above (in absolute value), a floating-point number.

Example 1

This is Mark’s example, and it shows a positive number rounding to nearest in round toward positive infinity mode:

1.100000000000000088817841970012523233890533447265626

It converts to this binary number, which has a non-terminating binary fraction portion:

1.00011001100110011001100110011001100110011001100110100000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001…

Rounded toward positive infinity the result should be

1.0001100110011001100110011001100110011001100110011011

In hexadecimal “%a” notation, this equals 0x1.199999999999bp+0.

strtod() converts it incorrectly to 0x1.199999999999ap+0, which is one ULP below the correctly rounded value. 0x1.199999999999ap+0 is the floating-point number equal to decimal 1.100000000000000088817841970012523233890533447265625, which differs from the input decimal value only in the last decimal digit (5 vs 6).

Example 2

This example shows a negative number rounding to nearest in round toward negative infinity mode:

-0.91276999999999997026378650843980722129344940185546876

It converts to this non-terminating binary number:

-0.111010011010101101001011011100101100010100011001011110000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001…

Rounded toward negative infinity the result should be -0x1.d35696e58a330p-1; strtod() converts it incorrectly to -0x1.d35696e58a32fp-1.

Examples: Decimal Numbers Just Below a Floating-Point Number

This section shows two examples that strtod() rounds incorrectly; they are very close to, but just below (in absolute value), a floating-point number.

Example 3

This example shows a negative number rounding to nearest in round toward positive infinity mode:

-266.240000000000009094947017729282379150390624

Rounded toward positive infinity it should be -0x1.0a3d70a3d70a3p+8; strtod() rounds it incorrectly to -0x1.0a3d70a3d70a4p+8.

Example 4

This example shows a positive number rounding to nearest in both round toward negative infinity and round toward zero modes:

8.255628858767918002472043289952338102302250764062685473021474535926245152950286865234374e-17

Rounded toward negative infinity and rounded toward zero it should be 0x1.7cb9433617c9bp-54; strtod() rounds it incorrectly — in both cases — to 0x1.7cb9433617c9cp-54.

Discussion

In my article “Visual C++ and GLIBC strtod() Ignore Rounding Mode” I said

The percentage of incorrect conversions for both Visual C++ and glibc ranged between 49.00% and 49.06%, suggesting they were not reacting to the rounding mode at all. (Why it wasn’t almost exactly 50% I don’t know.)

It might not have been 50% because, although I didn’t know it at the time, strtod()’s conversions weren’t 100% correct! (I’ll rerun those tests when dtoa.c is fixed.) (Update: the percentages are still near 49%, even with the fixed dtoa.c, so scratch that theory.)

Do You Use Directed Rounding?

If your project uses dtoa.c, could you let me know whether it enables directed conversions? That is, does it define the “Honor_FLT_ROUNDS” flag?

You should be aware that strtod() reacts to the rounding mode, although incorrectly, when Honor_FLT_ROUNDS is not defined.

Bug Report

I reported this to David Gay; I’ll keep you posted.

Update: dtoa.c is Fixed

David Gay updated dtoa.c with a fix, and it appears to work: all four examples above now round correctly, as did the millions of random examples I checked.