While running some of GCC’s string to double conversion testcases I discovered a bug in David Gay’s strtod(): it converts some very small subnormal numbers incorrectly. Unlike numbers 2^{-1075} or smaller, which should convert to zero under round-to-nearest/ties-to-even rounding, numbers between 2^{-1075} and 2^{-1074} should convert to 2^{-1074}, the smallest number representable in double-precision binary floating-point. strtod() correctly converts the former to 0, but it incorrectly converts the latter to 0 as well.

(**Update 11/25/13**: This bug has been fixed.)

## Example

Here’s an example of a number between 2^{-1074} and 2^{-1075} that strtod() gets wrong:

2.4703282292062327208828439643411068618252990130716238221279284125033775363510437593264991818081799618989828234772285886546332835517796989819938739800539093906315035659515570226392290858392449105184435931802849936536152500319370457678249219365623669863658480757001585769269903706311928279558551332927834338409351978015531246597263579574622766465272827220056374006485499977096599470454020828166226237857393450736339007967761930577506740176324673600968951340535537458516661134223766678604162159680461914467291840300530057530849048765391711386591646239524912623653881879636239373280423891018672348497668235089863388587925628302755995657524455507255189313690836254779186948667994968324049705821028513185451396213837722826145437693412532098591327667236328125001e-324

(This is the testcase named ‘d1c’ in GCC’s string to double test suite.)

This number is *very* close to 2^{-1075}; it is just a tiny bit bigger. Its first 1 bit is at bit position 1075, the rounding bit for subnormal numbers. That 1 bit is followed by 2506 zeros, and then by a non-zero portion. As such, it should round up to 2^{-1074}, or 0x0.0000000000001p-1022 in hexadecimal floating-point constant form.

## The Bug

I didn’t try to find the cause of the bug; I just reported it David Gay. ~~I will update this article when he responds.~~

### Update 11/25/13

As usual, Dave quickly provided a fix; this is the description of it in the dtoa.c changes file:

20131122 dtoa.c: fix a possible glitch with strtod in deciding whether a decimal value less in absolute value than the smallest denormal should be rounded to zero. (The fix is to force use of bigcomp() when at the bottom of the exponent range and we otherwise would have returned 0.)

It is a simple fix; these two highlighted lines

#ifdef Avoid_Underflow if (bc.scale && y <= 2*P*Exp_msk1) { if (aadj <= 0x7fffffff) { if ((z = aadj) <= 0) z = 1; aadj = z; aadj1 = bc.dsign ? aadj : -aadj; } dval(&aadj2) = aadj1; word0(&aadj2) += (2*P+1)*Exp_msk1 - y; aadj1 = dval(&aadj2); adj.d = aadj1 * ulp(&rv); dval(&rv) += adj.d; if (rv.d == 0.) #ifdef NO_STRTOD_BIGCOMP goto undfl; #else { if (bc.nd > nd) bc.dsign = 1; break; } #endif ...

were changed to this single line (highlighted):

#ifdef Avoid_Underflow if (bc.scale && y <= 2*P*Exp_msk1) { if (aadj <= 0x7fffffff) { if ((z = aadj) <= 0) z = 1; aadj = z; aadj1 = bc.dsign ? aadj : -aadj; } dval(&aadj2) = aadj1; word0(&aadj2) += (2*P+1)*Exp_msk1 - y; aadj1 = dval(&aadj2); adj.d = aadj1 * ulp(&rv); dval(&rv) += adj.d; if (rv.d == 0.) #ifdef NO_STRTOD_BIGCOMP goto undfl; #else { req_bigcomp = 1; break; } #endif ...

The fixed dtoa.c is on ampl.com.

It is a bug from bigcomp. #define NO_STRTOD_BIGCOMP will give correct result.

Water,

Yes, good catch. But interestingly it’s not bigcomp per se, but the code added for bigcomp that decides whether to call bigcomp.

i = 1 because strtod() thinks the approximation is more than 1/2 ulp above the decimal input, so req_bigcomp is never set. Since the approximation is scaled this made me think it has to do with the interaction of the “Avoid_Underflow” logic. I quickly confirmed this by defining the “NO_IEEE_Scale” flag. This makes the bug disappear as well.

11/25/13: Article updated with description of fix.

When I try to read Gay’s code, I don’t think it is bug free. It is hard to read and such code always has bugs.

I find one bug with #define IBM. The number 12345678912345678 is parsed wrong when calculate rv. see following code.

if (k > DBL_DIG)

oldinexact = get_inexact();

#endif

dval(&rv) = tens[k – 9] * dval(&rv) + z;

Water,

I hear what you’re saying.

I never tried any arithmetic other than IEEE. You can report the bug to Dave Gay directly — dmg at acm dot org. He responds quickly to bug reports.

Thanks, Rick. Actually it is bug in the process of parsing decimal strings to rough value rv. Nothing concerning IBM doubles. Anyway, I don’t have access to IBM machine to unable to check it in running.

It looks like that piece of code is under the ‘SET_INEXACT’ flag and not the ‘IBM’ flag. The comments say SET_INEXACT applies to IEEE arithmetic.

Rick, don’t mind. I think the bug is like this.

On IBM machine, DBL_DIG is 16. for input “12345678912345678”.

Line 2573:

for(nd = nf = 0; (c = *s) >= ‘0’ && c <= '9'; nd++, s++)

if (nd < 9)

y = 10*y + c – '0';

else if (nd < 16)

z = 10*z + c – '0';

y has 9 digits "123456789", z has 7 digits "1234567".

Then 2726: dval(&rv) = tens[k – 9] * dval(&rv) + z;

k should be 16. but it is 17 = DBL_DIG + 1.

Yes, I see it now.

Rick. I can say a little more to make you understand why I find this bug. I am very suspicious about the line:

if (nd < 9)

y = 10*y + c – '0';

else if (nd = ‘0’ && c <= '9'; nd++, s++)

if (nd < 9)

y = 10*y + c – '0';

else if (nd < DBL_DIG + 2)

z = 10*z + c – '0';

place 2: for(i = 1; i < nz; i++)

if (nd++ < 9)

y *= 10;

else if (nd <= DBL_DIG + 2)

z *= 10;

if (nd++ < 9)

y = 10*y + c;

else if (nd <= DBL_DIG + 2)

z = 10*z + c;

place 3: if (!nd0)

nd0 = nd;

k = nd < DBL_DIG + 2 ? nd : DBL_DIG + 2;

Let me know what you think about this. It can reduce guess error to about 6ULP maximum, instead of original 10ULP.

It seems blog post has limits of chars. I will send email instead.