This is the mail archive of the
newlib@sourceware.org
mailing list for the newlib project.
Re: RFC/RFA: strtod() magnitude problem
- From: nick clifton <nickc at redhat dot com>
- To: Christian Bruel <christian dot bruel at st dot com>, Corinna Vinschen <vinschen at redhat dot com>
- Cc: newlib at sourceware dot org
- Date: Tue, 23 Apr 2013 17:21:15 +0100
- Subject: Re: RFC/RFA: strtod() magnitude problem
- References: <87vc7d3a4g dot fsf at redhat dot com> <20130423120946 dot GB26397 at calimero dot vinschen dot de>
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;