Gay’s strtod() Returns Zero For Inputs Just Above 2^-1075

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.