Re: RFC/RFA: strtod() magnitude problem

[Corinna Vinschen <vinschen at redhat dot com>]

Hi Nick,

On Apr 23 09:49, Nick Clifton wrote:
> Hi Guys,
> 
>   The strtod() function appears to misbehave on targets where
>   DOUBLE_IS_32_BITS is true.  For example:
> 
>     strtod ("123456789.10", &endptr)
> 
>   return will return a value or 1234567.125 rather than 123456792.0, ie
>   two orders of magnitude too small.

There's also another problem with strtod in conjunction with 32 bit
doubles, see, for instance the thread starting at
http://sourceware.org/ml/newlib/2011/msg00178.html and the followup
reversion of the patch here:
http://sourceware.org/ml/newlib/2012/msg00576.html

>   Unfortunately endptr will not be left pointing at "89.10", indicating
>   that the remaining digits had not been consumed, but instead it will
>   be point at the NUL following ".10".  Thus you cannot examine endptr
>   to determine that the conversion has failed.
> 
>   The problem happens because the preprocessor constant DBL_DIG is set
>   to 6, (for these small sized doubles targets), so only the first 7
>   digits are considered for conversion.  Presumably the idea is that any
>   value returned by strtod() must be able to be converted back into a
>   similar string by atod().
> 
>   This problem is made worse by the fact that scanf() uses strtod() to
>   perform ascii to floating point conversions so for example
>   sscanf ("123456789.10") fails as well.
> 
>   I am not sure if the behaviour of strtod() is considered to be correct
>   or not.

It's not.  This certainly requires a change in strtod.  Patching scanf
doesn't sounds like the right solution here.  The problem shouldn't
have happened.

Can you try if the below patch fixes the issue?  It's just a manual
update of strtod (or rather _strtod_r) to the latest implementation
from NetBSD.  There were a few changes which might explain your problem
as well as the one encountered by Christian Bruel.

Christian, can you please try this patch in your scenario, too?

Additionally, when looking closely into the code, there seem to be a few
places where an `#ifndef _DOUBLE_IS_32BITS' might be missing, and that's
not just due to the below patch.  Rather it seems there are a few
unguarded tests for 'dword1(rv)' which look wrong in 32 bit double mode.
Can you please have a look?


Thanks,
Corinna


	* libc/stdlib/strtod.c: Manual update to latest algorithm from
	NetBSD.


Index: libc/stdlib/strtod.c
===================================================================
RCS file: /cvs/src/src/newlib/libc/stdlib/strtod.c,v
retrieving revision 1.19
diff -u -p -r1.19 strtod.c
--- libc/stdlib/strtod.c	19 Dec 2012 10:16:00 -0000	1.19
+++ libc/stdlib/strtod.c	23 Apr 2013 12:05:08 -0000
@@ -128,10 +128,10 @@ THIS SOFTWARE.
 #ifndef NO_IEEE_Scale
 #define Avoid_Underflow
 #undef tinytens
-/* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */
+/* The factor of 2^106 in tinytens[4] helps us avoid setting the underflow */
 /* flag unnecessarily.  It leads to a song and dance at the end of strtod. */
 static _CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,
-		9007199254740992.e-256
+		9007199254740992.*9007199254740992.e-256
 		};
 #endif
 #endif
@@ -144,6 +144,26 @@ static _CONST double tinytens[] = { 1e-1
 #define Rounding Flt_Rounds
 #endif
 
+#ifdef Avoid_Underflow /*{*/
+ static double
+_DEFUN (sulp, (x, scale),
+       	U x _AND
+	int scale)
+{
+        U u;
+        double rv;
+        int i;
+
+        rv = ulp(dval(x));
+        if (!scale || (i = 2*P + 1 - ((dword0(x) & Exp_mask) >> Exp_shift)) <= 0)
+                return rv; /* Is there an example where i <= 0 ? */
+        dword0(u) = Exp_1 + (i << Exp_shift);
+        dword1(u) = 0;
+        return rv * u.d;
+        }
+#endif /*}*/
+
+
 #ifndef NO_HEX_FP
 
 static void
@@ -221,7 +241,10 @@ _DEFUN (_strtod_r, (ptr, s00, se),
 	U aadj1, rv, rv0;
 	Long L;
 	__ULong y, z;
-	_Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
+	_Bigint *bb = NULL, *bb1, *bd = NULL, *bd0, *bs = NULL, *delta = NULL;
+#ifdef Avoid_Underflow
+	__ULong Lsb, Lsb1;
+#endif
 #ifdef SET_INEXACT
 	int inexact, oldinexact;
 #endif
@@ -279,6 +302,8 @@ _DEFUN (_strtod_r, (ptr, s00, se),
 			switch((i = gethex(ptr, &s, &fpi1, &exp, &bb, sign)) & STRTOG_Retmask) {
 			  case STRTOG_NoNumber:
 				s = s00;
+				sign = 0;
+				/* FALLTHROUGH */
 			  case STRTOG_Zero:
 				break;
 			  default:
@@ -299,14 +324,11 @@ _DEFUN (_strtod_r, (ptr, s00, se),
 		}
 	s0 = s;
 	y = z = 0;
-	for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++) {
-		if (nd < DBL_DIG + 1) {
-			if (nd < 9)
-				y = 10*y + c - '0';
-			else
-				z = 10*z + c - '0';
-		}
-        }
+	for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
+		if (nd < 9)
+			y = 10*y + c - '0';
+		else
+			z = 10*z + c - '0';
 	nd0 = nd;
 	if (strncmp (s, _localeconv_r (ptr)->decimal_point,
 		     strlen (_localeconv_r (ptr)->decimal_point)) == 0)
@@ -329,20 +351,15 @@ _DEFUN (_strtod_r, (ptr, s00, se),
 			nz++;
 			if (c -= '0') {
 				nf += nz;
-				for(i = 1; i < nz; i++) {
-					if (nd++ <= DBL_DIG + 1) {
-						if (nd < 10)
-							y *= 10;
-						else
-							z *= 10;
-					}
-				}
-				if (nd++ <= DBL_DIG + 1) {
-					if (nd < 10)
-						y = 10*y + c;
-					else
-						z = 10*z + c;
-				}
+				for(i = 1; i < nz; i++)
+					if (nd++ < 9)
+						y *= 10;
+					else if (nd <= DBL_DIG + 1)
+						z *= 10;
+				if (nd++ < 9)
+					y = 10*y + c;
+				else if (nd <= DBL_DIG + 1)
+					z = 10*z + c;
 				nz = 0;
 				}
 			}
@@ -691,12 +708,20 @@ _DEFUN (_strtod_r, (ptr, s00, se),
 	/* Put digits into bd: true value = bd * 10^e */
 
 	bd0 = s2b(ptr, s0, nd0, nd, y);
+	if (bd0 == NULL)
+		goto ovfl;
 
 	for(;;) {
 		bd = Balloc(ptr,bd0->_k);
+		if (bd == NULL)
+			goto ovfl;
 		Bcopy(bd, bd0);
 		bb = d2b(ptr,dval(rv), &bbe, &bbbits);	/* rv = bb * 2^bbe */
+		if (bb == NULL)
+			goto ovfl;
 		bs = i2b(ptr,1);
+		if (bs == NULL)
+			goto ovfl;
 
 		if (e >= 0) {
 			bb2 = bb5 = 0;
@@ -716,12 +741,19 @@ _DEFUN (_strtod_r, (ptr, s00, se),
 			bs2++;
 #endif
 #ifdef Avoid_Underflow
+		Lsb = LSB;
+		Lsb1 = 0;
 		j = bbe - scale;
 		i = j + bbbits - 1;	/* logb(rv) */
-		if (i < Emin)	/* denormal */
-			j += P - Emin;
-		else
-			j = P + 1 - bbbits;
+		j = P + 1 - bbbits;
+		if (i < Emin) {	/* denormal */
+			i = Emin - i;
+			j -= i;
+			if (i < 32)
+				Lsb <<= i;
+			else
+				Lsb1 = Lsb << (i-32);
+			}
 #else /*Avoid_Underflow*/
 #ifdef Sudden_Underflow
 #ifdef IBM
@@ -753,19 +785,37 @@ _DEFUN (_strtod_r, (ptr, s00, se),
 			}
 		if (bb5 > 0) {
 			bs = pow5mult(ptr, bs, bb5);
+			if (bs == NULL)
+				goto ovfl;
 			bb1 = mult(ptr, bs, bb);
+			if (bb1 == NULL)
+				goto ovfl;
 			Bfree(ptr, bb);
 			bb = bb1;
 			}
-		if (bb2 > 0)
+		if (bb2 > 0) {
 			bb = lshift(ptr, bb, bb2);
-		if (bd5 > 0)
+			if (bb == NULL)
+				goto ovfl;
+			}
+		if (bd5 > 0) {
 			bd = pow5mult(ptr, bd, bd5);
-		if (bd2 > 0)
+			if (bd == NULL)
+				goto ovfl;
+			}
+		if (bd2 > 0) {
 			bd = lshift(ptr, bd, bd2);
-		if (bs2 > 0)
+			if (bd == NULL)
+				goto ovfl;
+			}
+		if (bs2 > 0) {
 			bs = lshift(ptr, bs, bs2);
+			if (bs == NULL)
+				goto ovfl;
+			}
 		delta = diff(ptr, bb, bd);
+		if (delta == NULL)
+			goto ovfl;
 		dsign = delta->_sign;
 		delta->_sign = 0;
 		i = cmp(delta, bs);
@@ -852,7 +902,9 @@ _DEFUN (_strtod_r, (ptr, s00, se),
 #endif /*Sudden_Underflow*/
 #endif /*Avoid_Underflow*/
 			adj *= ulp(dval(rv));
-			if (dsign)
+			if (dsign) {
+				if (dword0(rv) == Big0 && dword1(rv) == Big1)
+					goto ovfl;
 				dval(rv) += adj;
 			else
 				dval(rv) -= adj;
@@ -902,6 +954,8 @@ _DEFUN (_strtod_r, (ptr, s00, se),
 #endif
 						   0xffffffff)) {
 					/*boundary case -- increment exponent*/
+					if (dword0(rv) == Big0 && dword1(rv) == Big1)
+						goto ovfl;
 					dword0(rv) = (dword0(rv) & Exp_mask)
 						+ Exp_msk1
 #ifdef IBM
@@ -960,14 +1014,31 @@ _DEFUN (_strtod_r, (ptr, s00, se),
 #endif
 				}
 #ifndef ROUND_BIASED
+#ifdef Avoid_Underflow
+			if (Lsb1) {
+				if (!(dword0(rv) & Lsb1))
+					break;
+				}
+			else if (!(dword1(rv) & Lsb))
+				break;
+#else
 			if (!(dword1(rv) & LSB))
 				break;
 #endif
+#endif
 			if (dsign)
+#ifdef Avoid_Underflow
+				dval(rv) += sulp(rv, scale);
+#else
 				dval(rv) += ulp(dval(rv));
+#endif
 #ifndef ROUND_BIASED
 			else {
+#ifdef Avoid_Underflow
+				dval(rv) -= sulp(rv, scale);
+#else
 				dval(rv) -= ulp(dval(rv));
+#endif
 #ifndef Sudden_Underflow
 				if (!dval(rv))
 					goto undfl;
@@ -1044,7 +1115,7 @@ _DEFUN (_strtod_r, (ptr, s00, se),
 #ifdef Avoid_Underflow
 			if (scale && y <= 2*P*Exp_msk1) {
 				if (aadj <= 0x7fffffff) {
-					if ((z = aadj) <= 0)
+					if ((z = aadj) == 0)
 						z = 1;
 					aadj = z;
 					dval(aadj1) = dsign ? aadj : -aadj;



-- 
Corinna Vinschen
Cygwin Maintainer
Red Hat