# 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 ((z = aadj) <= 0)
z = 1;
}
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 ((z = aadj) <= 0)
z = 1;
}
if (rv.d == 0.)
#ifdef NO_STRTOD_BIGCOMP
goto undfl;
#else
{
req_bigcomp = 1;
break;
}
#endif
...

```

The fixed dtoa.c is on ampl.com. 1. water qian says:

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

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

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

4. water qian says:

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;

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

6. water qian says:

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.

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

8. water qian says:

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.

9. water qian says:

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;