RFC/RFA: strtod() magnitude problem
Christian Bruel
christian.bruel@st.com
Thu Apr 25 07:19:00 GMT 2013
On 04/23/2013 02:09 PM, Corinna Vinschen wrote:
> 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?
Hello Corinna,
yes, my -m4-single-only tests are OK with your patch.
I just had to guard dword1 that is undef for DOUBLE_IS_32BITS, (see
bellow and stdlib/mprec.h), fixing a "error: lvalue required as left
operand of assignment" compilation error
Apart that, no regressions. Thanks !
Christian
--- strtod.c 2013-04-25 08:56:29.000000000 +0200
+++ strtod.c~ 2013-04-25 08:52:59.000000000 +0200
@@ -158,9 +158,7 @@
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 /*}*/
>
> 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;
>
>
>
More information about the Newlib
mailing list