Re: RFC/RFA: strtod() magnitude problem

[nick clifton <nickc at redhat dot com>]

Hi Corinna,

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

It does - although I had to tweak it slightly.  (Tweaked version 
attached).  There was one place I found where the dword1() macro was 
being used unsafely.  The attached patch includes a fix for this. There 
are other places where dword1() is used without checking 
DOUBLE_IS_32BITS, but these places are safe.

One point - this part of your patch:

>   static _CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,
> -		9007199254740992.e-256
> +		9007199254740992.*9007199254740992.e-256
>   		};

appears to be a typo.  I have removed it in my revised patch, but maybe 
it needs fixing ?

Cheers
  Nick



Index: libc/stdlib/strtod.c
===================================================================
RCS file: /cvs/cvsfiles/devo/newlib/libc/stdlib/strtod.c,v
retrieving revision 1.31
diff -c -3 -p -r1.31 strtod.c
*** libc/stdlib/strtod.c	27 Jan 2013 12:46:12 -0000	1.31
--- libc/stdlib/strtod.c	23 Apr 2013 16:10:18 -0000
*************** THIS SOFTWARE.
*** 128,138 ****
  #ifndef NO_IEEE_Scale
  #define Avoid_Underflow
  #undef tinytens
! /* The factor of 2^53 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
! 		};
  #endif
  #endif
  
--- 128,137 ----
  #ifndef NO_IEEE_Scale
  #define Avoid_Underflow
  #undef tinytens
! /* 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 };
  #endif
  #endif
  
*************** static _CONST double tinytens[] = { 1e-1
*** 144,149 ****
--- 143,170 ----
  #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);
+ #ifndef _DOUBLE_IS_32BITS
+         dword1(u) = 0;
+ #endif
+         return rv * u.d;
+         }
+ #endif /*}*/
+ 
+ 
  #ifndef NO_HEX_FP
  
  static void
*************** _DEFUN (_strtod_r, (ptr, s00, se),
*** 221,227 ****
  	U aadj1, rv, rv0;
  	Long L;
  	__ULong y, z;
! 	_Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
  #ifdef SET_INEXACT
  	int inexact, oldinexact;
  #endif
--- 242,251 ----
  	U aadj1, rv, rv0;
  	Long L;
  	__ULong y, z;
! 	_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
*************** _DEFUN (_strtod_r, (ptr, s00, se),
*** 279,284 ****
--- 303,310 ----
  			switch((i = gethex(ptr, &s, &fpi1, &exp, &bb, sign)) & STRTOG_Retmask) {
  			  case STRTOG_NoNumber:
  				s = s00;
+ 				sign = 0;
+ 				/* FALLTHROUGH */
  			  case STRTOG_Zero:
  				break;
  			  default:
*************** _DEFUN (_strtod_r, (ptr, s00, se),
*** 299,312 ****
  		}
  	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';
! 		}
!         }
  	nd0 = nd;
  	if (strncmp (s, _localeconv_r (ptr)->decimal_point,
  		     strlen (_localeconv_r (ptr)->decimal_point)) == 0)
--- 325,335 ----
  		}
  	s0 = s;
  	y = z = 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)
*************** _DEFUN (_strtod_r, (ptr, s00, se),
*** 329,348 ****
  			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;
! 				}
  				nz = 0;
  				}
  			}
--- 352,366 ----
  			nz++;
  			if (c -= '0') {
  				nf += nz;
! 				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;
  				}
  			}
*************** _DEFUN (_strtod_r, (ptr, s00, se),
*** 691,702 ****
--- 709,728 ----
  	/* 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;
*************** _DEFUN (_strtod_r, (ptr, s00, se),
*** 716,727 ****
  			bs2++;
  #endif
  #ifdef Avoid_Underflow
  		j = bbe - scale;
  		i = j + bbbits - 1;	/* logb(rv) */
! 		if (i < Emin)	/* denormal */
! 			j += P - Emin;
! 		else
! 			j = P + 1 - bbbits;
  #else /*Avoid_Underflow*/
  #ifdef Sudden_Underflow
  #ifdef IBM
--- 742,760 ----
  			bs2++;
  #endif
  #ifdef Avoid_Underflow
+ 		Lsb = LSB;
+ 		Lsb1 = 0;
  		j = bbe - scale;
  		i = j + bbbits - 1;	/* logb(rv) */
! 		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
*************** _DEFUN (_strtod_r, (ptr, s00, se),
*** 753,771 ****
  			}
  		if (bb5 > 0) {
  			bs = pow5mult(ptr, bs, bb5);
  			bb1 = mult(ptr, bs, bb);
  			Bfree(ptr, bb);
  			bb = bb1;
  			}
! 		if (bb2 > 0)
  			bb = lshift(ptr, bb, bb2);
! 		if (bd5 > 0)
  			bd = pow5mult(ptr, bd, bd5);
! 		if (bd2 > 0)
  			bd = lshift(ptr, bd, bd2);
! 		if (bs2 > 0)
  			bs = lshift(ptr, bs, bs2);
  		delta = diff(ptr, bb, bd);
  		dsign = delta->_sign;
  		delta->_sign = 0;
  		i = cmp(delta, bs);
--- 786,822 ----
  			}
  		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) {
  			bb = lshift(ptr, bb, bb2);
! 			if (bb == NULL)
! 				goto ovfl;
! 			}
! 		if (bd5 > 0) {
  			bd = pow5mult(ptr, bd, bd5);
! 			if (bd == NULL)
! 				goto ovfl;
! 			}
! 		if (bd2 > 0) {
  			bd = lshift(ptr, bd, bd2);
! 			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);
*************** _DEFUN (_strtod_r, (ptr, s00, se),
*** 789,795 ****
  				else if (!dsign) {
  					adj = -1.;
  					if (!dword1(rv)
! 					 && !(dword0(rv) & Frac_mask)) {
  						y = dword0(rv) & Exp_mask;
  #ifdef Avoid_Underflow
  						if (!scale || y > 2*P*Exp_msk1)
--- 840,846 ----
  				else if (!dsign) {
  					adj = -1.;
  					if (!dword1(rv)
! 					    && !(dword0(rv) & Frac_mask)) {
  						y = dword0(rv) & Exp_mask;
  #ifdef Avoid_Underflow
  						if (!scale || y > 2*P*Exp_msk1)
*************** _DEFUN (_strtod_r, (ptr, s00, se),
*** 852,858 ****
  #endif /*Sudden_Underflow*/
  #endif /*Avoid_Underflow*/
  			adj *= ulp(dval(rv));
! 			if (dsign)
  				dval(rv) += adj;
  			else
  				dval(rv) -= adj;
--- 903,911 ----
  #endif /*Sudden_Underflow*/
  #endif /*Avoid_Underflow*/
  			adj *= ulp(dval(rv));
! 			if (dsign) {
! 				if (dword0(rv) == Big0 && dword1(rv) == Big1)
! 					goto ovfl;
  				dval(rv) += adj;
  			else
  				dval(rv) -= adj;
*************** _DEFUN (_strtod_r, (ptr, s00, se),
*** 902,907 ****
--- 955,962 ----
  #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
*************** _DEFUN (_strtod_r, (ptr, s00, se),
*** 960,973 ****
--- 1015,1045 ----
  #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;
*************** _DEFUN (_strtod_r, (ptr, s00, se),
*** 1044,1050 ****
  #ifdef Avoid_Underflow
  			if (scale && y <= 2*P*Exp_msk1) {
  				if (aadj <= 0x7fffffff) {
! 					if ((z = aadj) <= 0)
  						z = 1;
  					aadj = z;
  					dval(aadj1) = dsign ? aadj : -aadj;
--- 1116,1122 ----
  #ifdef Avoid_Underflow
  			if (scale && y <= 2*P*Exp_msk1) {
  				if (aadj <= 0x7fffffff) {
! 					if ((z = aadj) == 0)
  						z = 1;
  					aadj = z;
  					dval(aadj1) = dsign ? aadj : -aadj;