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.

Dingbat
RSS feed icon
RSS e-mail icon

11 comments

  1. Water,

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

    #ifndef NO_STRTOD_BIGCOMP /*{*/
    		if (bc.nd > nd && i < = 0) {
    			if (bc.dsign) {
    				/* Must use bigcomp(). */
    				req_bigcomp = 1;
    				break;
    				}
    

    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.

  2. 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;

  3. 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.

  4. 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.

  5. 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.

  6. 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.

  7. 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.

Leave a Reply

Your email address will not be published. Required fields are marked *

(Cookies must be enabled to leave a comment...it reduces spam.)