This is the mail archive of the
newlib@sourceware.org
mailing list for the newlib project.
Re: RFC/RFA: strtod() magnitude problem
- From: Sebastian Huber <sebastian dot huber at embedded-brains dot de>
- To: newlib at sourceware dot org
- Date: Wed, 24 Apr 2013 09:49:40 +0200
- Subject: Re: RFC/RFA: strtod() magnitude problem
- References: <87vc7d3a4g dot fsf at redhat dot com> <20130423120946 dot GB26397 at calimero dot vinschen dot de> <5176B4FB dot 8060308 at redhat dot com>
On 04/23/2013 06:21 PM, nick clifton wrote:
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 ?
Attached is a patch from FreeBSD that introduced this change. I looks pretty
intentional.
--
Sebastian Huber, embedded brains GmbH
Address : Dornierstr. 4, D-82178 Puchheim, Germany
Phone : +49 89 189 47 41-16
Fax : +49 89 189 47 41-09
E-Mail : sebastian.huber@embedded-brains.de
PGP : Public key available on request.
Diese Nachricht ist keine geschÃftliche Mitteilung im Sinne des EHUG.
>From 2732388653bc1d73e3e22df8de7d745845b3f37f Mon Sep 17 00:00:00 2001
From: das <das@FreeBSD.org>
Date: Wed, 3 Sep 2008 07:23:57 +0000
Subject: [PATCH] Merge gdtoa 20080831. This fixes several bugs, including an
infinite loop pointed out by cognet@ that occurs when
calling strtod() with a string representing a number
between DBL_MAX and 2*DBL_MAX, when the rounding mode is
anything other than the default.
---
contrib/gdtoa/README | 10 +++++
contrib/gdtoa/dtoa.c | 39 ++++++++++++--------
contrib/gdtoa/gdtoa.h | 6 ++--
contrib/gdtoa/gdtoaimp.h | 22 +++++++-----
contrib/gdtoa/gethex.c | 79 +++++++++++++++++++++++++++++++++++-----
contrib/gdtoa/strtoIg.c | 16 +++++---
contrib/gdtoa/strtod.c | 60 +++++++++++++++++++-----------
contrib/gdtoa/strtodg.c | 47 +++++++++++++++++++-----
contrib/gdtoa/test/README | 5 +++
contrib/gdtoa/test/f.out | 40 +++++++++++++++-----
contrib/gdtoa/test/getround.c | 9 ++++-
contrib/gdtoa/test/xsum0.out | 6 ++--
contrib/gdtoa/xsum0.out | 18 +++++-----
13 files changed, 259 insertions(+), 98 deletions(-)
diff --git a/contrib/gdtoa/README b/contrib/gdtoa/README
index cf1947a..95a40aa 100644
--- a/contrib/gdtoa/README
+++ b/contrib/gdtoa/README
@@ -332,5 +332,15 @@ Compiling g__fmt.c, strtod.c, and strtodg.c with -DUSE_LOCALE causes
the decimal-point character to be taken from the current locale; otherwise
it is '.'.
+Source files dtoa.c and strtod.c in this directory are derived from
+netlib's "dtoa.c from fp" and are meant to function equivalently.
+When compiled with Honor_FLT_ROUNDS #defined (on systems that provide
+FLT_ROUNDS and fegetround() as specified in the C99 standard), they
+honor the current rounding mode. Because FLT_ROUNDS is buggy on some
+(Linux) systems -- not reflecting calls on fesetround(), as the C99
+standard says it should -- when Honor_FLT_ROUNDS is #defined, the
+current rounding mode is obtained from fegetround() rather than from
+FLT_ROUNDS, unless Trust_FLT_ROUNDS is also #defined.
+
Please send comments to David M. Gay (dmg at acm dot org, with " at "
changed at "@" and " dot " changed to ".").
diff --git a/contrib/gdtoa/dtoa.c b/contrib/gdtoa/dtoa.c
index e808cc1..48fdf5e 100644
--- a/contrib/gdtoa/dtoa.c
+++ b/contrib/gdtoa/dtoa.c
@@ -66,7 +66,6 @@ THIS SOFTWARE.
*/
#ifdef Honor_FLT_ROUNDS
-#define Rounding rounding
#undef Check_FLT_ROUNDS
#define Check_FLT_ROUNDS
#else
@@ -127,12 +126,22 @@ dtoa
Bigint *b, *b1, *delta, *mlo, *mhi, *S;
double d2, ds, eps;
char *s, *s0;
-#ifdef Honor_FLT_ROUNDS
- int rounding;
-#endif
#ifdef SET_INEXACT
int inexact, oldinexact;
#endif
+#ifdef Honor_FLT_ROUNDS /*{*/
+ int Rounding;
+#ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */
+ Rounding = Flt_Rounds;
+#else /*}{*/
+ Rounding = 1;
+ switch(fegetround()) {
+ case FE_TOWARDZERO: Rounding = 0; break;
+ case FE_UPWARD: Rounding = 2; break;
+ case FE_DOWNWARD: Rounding = 3;
+ }
+#endif /*}}*/
+#endif /*}*/
#ifndef MULTIPLE_THREADS
if (dtoa_result) {
@@ -178,12 +187,12 @@ dtoa
inexact = 1;
#endif
#ifdef Honor_FLT_ROUNDS
- if ((rounding = Flt_Rounds) >= 2) {
+ if (Rounding >= 2) {
if (*sign)
- rounding = rounding == 2 ? 0 : 2;
+ Rounding = Rounding == 2 ? 0 : 2;
else
- if (rounding != 2)
- rounding = 0;
+ if (Rounding != 2)
+ Rounding = 0;
}
#endif
@@ -316,7 +325,7 @@ dtoa
s = s0 = rv_alloc(i);
#ifdef Honor_FLT_ROUNDS
- if (mode > 1 && rounding != 1)
+ if (mode > 1 && Rounding != 1)
leftright = 0;
#endif
@@ -453,7 +462,7 @@ dtoa
if (i == ilim) {
#ifdef Honor_FLT_ROUNDS
if (mode > 1)
- switch(rounding) {
+ switch(Rounding) {
case 0: goto ret1;
case 2: goto bump_up;
}
@@ -521,7 +530,7 @@ dtoa
spec_case = 0;
if ((mode < 2 || leftright)
#ifdef Honor_FLT_ROUNDS
- && rounding == 1
+ && Rounding == 1
#endif
) {
if (!word1(d) && !(word0(d) & Bndry_mask)
@@ -614,7 +623,7 @@ dtoa
#ifndef ROUND_BIASED
if (j1 == 0 && mode != 1 && !(word1(d) & 1)
#ifdef Honor_FLT_ROUNDS
- && rounding >= 1
+ && Rounding >= 1
#endif
) {
if (dig == '9')
@@ -642,7 +651,7 @@ dtoa
}
#ifdef Honor_FLT_ROUNDS
if (mode > 1)
- switch(rounding) {
+ switch(Rounding) {
case 0: goto accept_dig;
case 2: goto keep_dig;
}
@@ -660,7 +669,7 @@ dtoa
}
if (j1 > 0) {
#ifdef Honor_FLT_ROUNDS
- if (!rounding)
+ if (!Rounding)
goto accept_dig;
#endif
if (dig == '9') { /* possible if i == 1 */
@@ -703,7 +712,7 @@ dtoa
/* Round off last digit */
#ifdef Honor_FLT_ROUNDS
- switch(rounding) {
+ switch(Rounding) {
case 0: goto trimzeros;
case 2: goto roundoff;
}
diff --git a/contrib/gdtoa/gdtoa.h b/contrib/gdtoa/gdtoa.h
index ee6a9e5..ed0e020 100644
--- a/contrib/gdtoa/gdtoa.h
+++ b/contrib/gdtoa/gdtoa.h
@@ -74,9 +74,9 @@ typedef unsigned short UShort;
/* The following may be or-ed into one of the above values. */
- STRTOG_Neg = 0x08,
- STRTOG_Inexlo = 0x10,
- STRTOG_Inexhi = 0x20,
+ STRTOG_Neg = 0x08, /* does not affect STRTOG_Inexlo or STRTOG_Inexhi */
+ STRTOG_Inexlo = 0x10, /* returned result rounded toward zero */
+ STRTOG_Inexhi = 0x20, /* returned result rounded away from zero */
STRTOG_Inexact = 0x30,
STRTOG_Underflow= 0x40,
STRTOG_Overflow = 0x80
diff --git a/contrib/gdtoa/gdtoaimp.h b/contrib/gdtoa/gdtoaimp.h
index 0790b5e..916e156 100644
--- a/contrib/gdtoa/gdtoaimp.h
+++ b/contrib/gdtoa/gdtoaimp.h
@@ -132,13 +132,16 @@ THIS SOFTWARE.
* Infinity and NaN (case insensitively).
* When INFNAN_CHECK is #defined and No_Hex_NaN is not #defined,
* strtodg also accepts (case insensitively) strings of the form
- * NaN(x), where x is a string of hexadecimal digits and spaces;
- * if there is only one string of hexadecimal digits, it is taken
- * for the fraction bits of the resulting NaN; if there are two or
- * more strings of hexadecimal digits, each string is assigned
- * to the next available sequence of 32-bit words of fractions
- * bits (starting with the most significant), right-aligned in
- * each sequence.
+ * NaN(x), where x is a string of hexadecimal digits (optionally
+ * preceded by 0x or 0X) and spaces; if there is only one string
+ * of hexadecimal digits, it is taken for the fraction bits of the
+ * resulting NaN; if there are two or more strings of hexadecimal
+ * digits, each string is assigned to the next available sequence
+ * of 32-bit words of fractions bits (starting with the most
+ * significant), right-aligned in each sequence.
+ * Unless GDTOA_NON_PEDANTIC_NANCHECK is #defined, input "NaN(...)"
+ * is consumed even when ... has the wrong form (in which case the
+ * "(...)" is consumed but ignored).
* #define MULTIPLE_THREADS if the system offers preemptively scheduled
* multiple threads. In this case, you must provide (or suitably
* #define) two locks, acquired by ACQUIRE_DTOA_LOCK(n) and freed
@@ -150,7 +153,7 @@ THIS SOFTWARE.
* dtoa. You may do so whether or not MULTIPLE_THREADS is #defined.
* #define IMPRECISE_INEXACT if you do not care about the setting of
* the STRTOG_Inexact bits in the special case of doing IEEE double
- * precision conversions (which could also be done by the strtog in
+ * precision conversions (which could also be done by the strtod in
* dtoa.c).
* #define NO_HEX_FP to disable recognition of C9x's hexadecimal
* floating-point constants.
@@ -204,6 +207,7 @@ extern Char *MALLOC ANSI((size_t));
#define INFNAN_CHECK
#define USE_LOCALE
#define Honor_FLT_ROUNDS
+#define Trust_FLT_ROUNDS
#undef IEEE_Arith
#undef Avoid_Underflow
@@ -597,7 +601,7 @@ extern void memcpy_D2A ANSI((void*, const void*, size_t));
extern int cmp ANSI((Bigint*, Bigint*));
extern void copybits ANSI((ULong*, int, Bigint*));
extern Bigint *d2b ANSI((double, int*, int*));
- extern int decrement ANSI((Bigint*));
+ extern void decrement ANSI((Bigint*));
extern Bigint *diff ANSI((Bigint*, Bigint*));
extern char *dtoa ANSI((double d, int mode, int ndigits,
int *decpt, int *sign, char **rve));
diff --git a/contrib/gdtoa/gethex.c b/contrib/gdtoa/gethex.c
index f130fd5..b3f7c50 100644
--- a/contrib/gdtoa/gethex.c
+++ b/contrib/gdtoa/gethex.c
@@ -45,7 +45,7 @@ gethex( CONST char **sp, FPI *fpi, Long *exp, Bigint **bp, int sign)
{
Bigint *b;
CONST unsigned char *decpt, *s0, *s, *s1;
- int esign, havedig, irv, k, n, nbits, up, zret;
+ int big, esign, havedig, irv, j, k, n, n0, nbits, up, zret;
ULong L, lostbits, *x;
Long e, e1;
#ifdef USE_LOCALE
@@ -56,6 +56,7 @@ gethex( CONST char **sp, FPI *fpi, Long *exp, Bigint **bp, int sign)
if (!hexdig['0'])
hexdig_init_D2A();
+ *bp = 0;
havedig = 0;
s0 = *(CONST unsigned char **)sp + 2;
while(s0[havedig] == '0')
@@ -90,10 +91,10 @@ gethex( CONST char **sp, FPI *fpi, Long *exp, Bigint **bp, int sign)
e = -(((Long)(s-decpt)) << 2);
pcheck:
s1 = s;
+ big = esign = 0;
switch(*s) {
case 'p':
case 'P':
- esign = 0;
switch(*++s) {
case '-':
esign = 1;
@@ -106,17 +107,65 @@ gethex( CONST char **sp, FPI *fpi, Long *exp, Bigint **bp, int sign)
break;
}
e1 = n - 0x10;
- while((n = hexdig[*++s]) !=0 && n <= 0x19)
+ while((n = hexdig[*++s]) !=0 && n <= 0x19) {
+ if (e1 & 0xf8000000)
+ big = 1;
e1 = 10*e1 + n - 0x10;
+ }
if (esign)
e1 = -e1;
e += e1;
}
*sp = (char*)s;
- if (zret) {
- if (!havedig)
- *sp = s0 - 1;
+ if (!havedig)
+ *sp = s0 - 1;
+ if (zret)
return STRTOG_Zero;
+ if (big) {
+ if (esign) {
+ switch(fpi->rounding) {
+ case FPI_Round_up:
+ if (sign)
+ break;
+ goto ret_tiny;
+ case FPI_Round_down:
+ if (!sign)
+ break;
+ goto ret_tiny;
+ }
+ goto retz;
+ ret_tiny:
+ b = Balloc(0);
+ b->wds = 1;
+ b->x[0] = 1;
+ goto dret;
+ }
+ switch(fpi->rounding) {
+ case FPI_Round_near:
+ goto ovfl1;
+ case FPI_Round_up:
+ if (!sign)
+ goto ovfl1;
+ goto ret_big;
+ case FPI_Round_down:
+ if (sign)
+ goto ovfl1;
+ goto ret_big;
+ }
+ ret_big:
+ nbits = fpi->nbits;
+ n0 = n = nbits >> kshift;
+ if (nbits & kmask)
+ ++n;
+ for(j = n, k = 0; j >>= 1; ++k);
+ *bp = b = Balloc(k);
+ b->wds = n;
+ for(j = 0; j < n0; ++j)
+ b->x[j] = ALL_ON;
+ if (n > n0)
+ b->x[j] = ULbits >> (ULbits - (nbits & kmask));
+ *exp = fpi->emin;
+ return STRTOG_Normal | STRTOG_Inexlo;
}
n = s1 - s0 - 1;
for(k = 0; n > 7; n >>= 1)
@@ -149,7 +198,7 @@ gethex( CONST char **sp, FPI *fpi, Long *exp, Bigint **bp, int sign)
k = n - 1;
if (x[k>>kshift] & 1 << (k & kmask)) {
lostbits = 2;
- if (k > 1 && any_on(b,k-1))
+ if (k > 0 && any_on(b,k))
lostbits = 3;
}
}
@@ -165,7 +214,10 @@ gethex( CONST char **sp, FPI *fpi, Long *exp, Bigint **bp, int sign)
if (e > fpi->emax) {
ovfl:
Bfree(b);
- *bp = 0;
+ ovfl1:
+#ifndef NO_ERRNO
+ errno = ERANGE;
+#endif
return STRTOG_Infinite | STRTOG_Overflow | STRTOG_Inexhi;
}
irv = STRTOG_Normal;
@@ -185,15 +237,22 @@ gethex( CONST char **sp, FPI *fpi, Long *exp, Bigint **bp, int sign)
case FPI_Round_down:
if (sign) {
one_bit:
- *exp = fpi->emin;
x[0] = b->wds = 1;
+ dret:
*bp = b;
+ *exp = fpi->emin;
+#ifndef NO_ERRNO
+ errno = ERANGE;
+#endif
return STRTOG_Denormal | STRTOG_Inexhi
| STRTOG_Underflow;
}
}
Bfree(b);
- *bp = 0;
+ retz:
+#ifndef NO_ERRNO
+ errno = ERANGE;
+#endif
return STRTOG_Zero | STRTOG_Inexlo | STRTOG_Underflow;
}
k = n - 1;
diff --git a/contrib/gdtoa/strtoIg.c b/contrib/gdtoa/strtoIg.c
index ec46a3e..90c88db 100644
--- a/contrib/gdtoa/strtoIg.c
+++ b/contrib/gdtoa/strtoIg.c
@@ -61,12 +61,16 @@ strtoIg(CONST char *s00, char **se, FPI *fpi, Long *exp, Bigint **B, int *rvp)
if (rv & STRTOG_Inexlo) {
swap = 0;
b1 = increment(b1);
- if (fpi->sudden_underflow
- && (rv & STRTOG_Retmask) == STRTOG_Zero) {
- b1->x[0] = 0;
- b1->x[nw1] = 1L << nb11;
- rv1 += STRTOG_Normal - STRTOG_Zero;
- rv1 &= ~STRTOG_Underflow;
+ if ((rv & STRTOG_Retmask) == STRTOG_Zero) {
+ if (fpi->sudden_underflow) {
+ b1->x[0] = 0;
+ b1->x[nw1] = 1L << nb11;
+ rv1 += STRTOG_Normal - STRTOG_Zero;
+ rv1 &= ~STRTOG_Underflow;
+ goto swapcheck;
+ }
+ rv1 &= STRTOG_Inexlo | STRTOG_Underflow | STRTOG_Zero;
+ rv1 |= STRTOG_Inexhi | STRTOG_Denormal;
goto swapcheck;
}
if (b1->wds > nw
diff --git a/contrib/gdtoa/strtod.c b/contrib/gdtoa/strtod.c
index 3703f94..69567f5 100644
--- a/contrib/gdtoa/strtod.c
+++ b/contrib/gdtoa/strtod.c
@@ -44,16 +44,15 @@ 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
#ifdef Honor_FLT_ROUNDS
-#define Rounding rounding
#undef Check_FLT_ROUNDS
#define Check_FLT_ROUNDS
#else
@@ -81,9 +80,19 @@ strtod
#ifdef SET_INEXACT
int inexact, oldinexact;
#endif
-#ifdef Honor_FLT_ROUNDS
- int rounding;
-#endif
+#ifdef Honor_FLT_ROUNDS /*{*/
+ int Rounding;
+#ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */
+ Rounding = Flt_Rounds;
+#else /*}{*/
+ Rounding = 1;
+ switch(fegetround()) {
+ case FE_TOWARDZERO: Rounding = 0; break;
+ case FE_UPWARD: Rounding = 2; break;
+ case FE_DOWNWARD: Rounding = 3;
+ }
+#endif /*}}*/
+#endif /*}*/
sign = nz0 = nz = decpt = 0;
dval(rv) = 0.;
@@ -109,7 +118,7 @@ strtod
}
break2:
if (*s == '0') {
-#ifndef NO_HEX_FP
+#ifndef NO_HEX_FP /*{{*/
{
static FPI fpi = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI };
Long exp;
@@ -118,16 +127,20 @@ strtod
case 'x':
case 'X':
{
-#if defined(FE_DOWNWARD) && defined(FE_TONEAREST) && defined(FE_TOWARDZERO) && defined(FE_UPWARD)
+#if defined(FE_DOWNWARD) && defined(FE_TONEAREST) && defined(FE_TOWARDZERO) && defined(FE_UPWARD) /*{{*/
FPI fpi1 = fpi;
+#ifdef Honor_FLT_ROUNDS /*{{*/
+ fpi1.rounding = Rounding;
+#else /*}{*/
switch(fegetround()) {
case FE_TOWARDZERO: fpi1.rounding = 0; break;
case FE_UPWARD: fpi1.rounding = 2; break;
case FE_DOWNWARD: fpi1.rounding = 3;
}
-#else
+#endif /*}}*/
+#else /*}{*/
#define fpi1 fpi
-#endif
+#endif /*}}*/
switch((i = gethex(&s, &fpi1, &exp, &bb, sign)) & STRTOG_Retmask) {
case STRTOG_NoNumber:
s = s00;
@@ -381,12 +394,12 @@ strtod
scale = 0;
#endif
#ifdef Honor_FLT_ROUNDS
- if ((rounding = Flt_Rounds) >= 2) {
+ if (Rounding >= 2) {
if (sign)
- rounding = rounding == 2 ? 0 : 2;
+ Rounding = Rounding == 2 ? 0 : 2;
else
- if (rounding != 2)
- rounding = 0;
+ if (Rounding != 2)
+ Rounding = 0;
}
#endif
#endif /*IEEE_Arith*/
@@ -405,7 +418,7 @@ strtod
/* Can't trust HUGE_VAL */
#ifdef IEEE_Arith
#ifdef Honor_FLT_ROUNDS
- switch(rounding) {
+ switch(Rounding) {
case 0: /* toward 0 */
case 3: /* toward -infinity */
word0(rv) = Big0;
@@ -536,7 +549,7 @@ strtod
bd2 -= bbe;
bs2 = bb2;
#ifdef Honor_FLT_ROUNDS
- if (rounding != 1)
+ if (Rounding != 1)
bs2++;
#endif
#ifdef Avoid_Underflow
@@ -594,7 +607,7 @@ strtod
delta->sign = 0;
i = cmp(delta, bs);
#ifdef Honor_FLT_ROUNDS
- if (rounding != 1) {
+ if (Rounding != 1) {
if (i < 0) {
/* Error is less than an ulp */
if (!delta->x[0] && delta->wds <= 1) {
@@ -604,7 +617,7 @@ strtod
#endif
break;
}
- if (rounding) {
+ if (Rounding) {
if (dsign) {
adj = 1.;
goto apply_adj;
@@ -650,10 +663,10 @@ strtod
if (adj < 1.)
adj = 1.;
if (adj <= 0x7ffffffe) {
- /* adj = rounding ? ceil(adj) : floor(adj); */
+ /* adj = Rounding ? ceil(adj) : floor(adj); */
y = adj;
if (y != adj) {
- if (!((rounding>>1) ^ dsign))
+ if (!((Rounding>>1) ^ dsign))
y++;
adj = y;
}
@@ -676,8 +689,11 @@ strtod
#endif /*Sudden_Underflow*/
#endif /*Avoid_Underflow*/
adj *= ulp(dval(rv));
- if (dsign)
+ if (dsign) {
+ if (word0(rv) == Big0 && word1(rv) == Big1)
+ goto ovfl;
dval(rv) += adj;
+ }
else
dval(rv) -= adj;
goto cont;
@@ -770,7 +786,7 @@ strtod
}
#endif /*Avoid_Underflow*/
L = (word0(rv) & Exp_mask) - Exp_msk1;
-#endif /*Sudden_Underflow}*/
+#endif /*Sudden_Underflow}}*/
word0(rv) = L | Bndry_mask1;
word1(rv) = 0xffffffff;
#ifdef IBM
diff --git a/contrib/gdtoa/strtodg.c b/contrib/gdtoa/strtodg.c
index cbdf4aa..9b73d95 100644
--- a/contrib/gdtoa/strtodg.c
+++ b/contrib/gdtoa/strtodg.c
@@ -89,7 +89,7 @@ increment(Bigint *b)
return b;
}
- int
+ void
#ifdef KR_headers
decrement(b) Bigint *b;
#else
@@ -119,7 +119,6 @@ decrement(Bigint *b)
*x++ = y & 0xffff;
} while(borrow && x < xe);
#endif
- return STRTOG_Inexlo;
}
static int
@@ -206,9 +205,9 @@ rvOK
goto ret;
}
switch(rd) {
- case 1:
+ case 1: /* round down (toward -Infinity) */
goto trunc;
- case 2:
+ case 2: /* round up (toward +Infinity) */
break;
default: /* round near */
k = bdif - 1;
@@ -330,7 +329,7 @@ strtodg
CONST char *s, *s0, *s1;
double adj, adj0, rv, tol;
Long L;
- ULong y, z;
+ ULong *b, *be, y, z;
Bigint *ab, *bb, *bb1, *bd, *bd0, *bs, *delta, *rvb, *rvb0;
irv = STRTOG_Zero;
@@ -822,10 +821,8 @@ strtodg
break;
if (dsign) {
rvb = increment(rvb);
- if ( (j = rvbits & kmask) !=0)
- j = ULbits - j;
- if (hi0bits(rvb->x[(rvb->wds - 1) >> kshift])
- != j)
+ j = kmask & (ULbits - (rvbits & kmask));
+ if (hi0bits(rvb->x[rvb->wds - 1]) != j)
rvbits++;
irv = STRTOG_Normal | STRTOG_Inexhi;
}
@@ -978,6 +975,29 @@ strtodg
Bfree(bd0);
Bfree(delta);
if (rve > fpi->emax) {
+ switch(fpi->rounding & 3) {
+ case FPI_Round_near:
+ goto huge;
+ case FPI_Round_up:
+ if (!sign)
+ goto huge;
+ break;
+ case FPI_Round_down:
+ if (sign)
+ goto huge;
+ }
+ /* Round to largest representable magnitude */
+ Bfree(rvb);
+ rvb = 0;
+ irv = STRTOG_Normal | STRTOG_Inexlo;
+ *exp = fpi->emax;
+ b = bits;
+ be = b + (fpi->nbits >> 5) + 1;
+ while(b < be)
+ *b++ = -1;
+ if ((j = fpi->nbits & 0x1f))
+ *--be >>= (32 - j);
+ goto ret;
huge:
rvb->wds = 0;
irv = STRTOG_Infinite | STRTOG_Overflow | STRTOG_Inexhi;
@@ -992,12 +1012,19 @@ strtodg
if (sudden_underflow) {
rvb->wds = 0;
irv = STRTOG_Underflow | STRTOG_Inexlo;
+#ifndef NO_ERRNO
+ errno = ERANGE;
+#endif
}
else {
irv = (irv & ~STRTOG_Retmask) |
(rvb->wds > 0 ? STRTOG_Denormal : STRTOG_Zero);
- if (irv & STRTOG_Inexact)
+ if (irv & STRTOG_Inexact) {
irv |= STRTOG_Underflow;
+#ifndef NO_ERRNO
+ errno = ERANGE;
+#endif
+ }
}
}
if (se)
diff --git a/contrib/gdtoa/test/README b/contrib/gdtoa/test/README
index 79b1c91..685fe13 100644
--- a/contrib/gdtoa/test/README
+++ b/contrib/gdtoa/test/README
@@ -55,7 +55,12 @@ logic (on double and double-double conversions).
Program strtodt tests strtod on some hard cases (in file testnos3)
posted by Fred Tydeman to comp.arch.arithmetic on 26 Feb. 1996.
+To get correct results on Intel (x86) systems, the rounding precision
+must be set to 53 bits. This can be done, e.g., by invoking
+fpinit_ASL(), whose source appears in
+http://www.netlib.org/ampl/solvers/fpinit.c .
These are simple test programs, not meant for exhaustive testing,
but for manually testing "interesting" cases. Paxson's testbase
is good for more exhaustive testing, in part with random inputs.
+See ftp://ftp.ee.lbl.gov/testbase-report.ps.Z .
diff --git a/contrib/gdtoa/test/f.out b/contrib/gdtoa/test/f.out
index ca8d6eb..d82e4e6 100644
--- a/contrib/gdtoa/test/f.out
+++ b/contrib/gdtoa/test/f.out
@@ -124,7 +124,9 @@ strtof consumes 9 bytes and returns 0 = #0
g_ffmt(0) gives 1 bytes: "0"
strtoIf returns 80, consuming 9 bytes.
-fI[0] == fI[1] == strtof
+fI[0] = #0 = 0
+fI[1] = #1 = 1.4012985e-45
+fI[0] == strtof
Input: 1.23e-320
@@ -132,7 +134,9 @@ strtof consumes 9 bytes and returns 0 = #0
g_ffmt(0) gives 1 bytes: "0"
strtoIf returns 80, consuming 9 bytes.
-fI[0] == fI[1] == strtof
+fI[0] = #0 = 0
+fI[1] = #1 = 1.4012985e-45
+fI[0] == strtof
Input: 1.23e-20
@@ -160,7 +164,9 @@ strtof consumes 15 bytes and returns 0 = #0
g_ffmt(0) gives 1 bytes: "0"
strtoIf returns 80, consuming 15 bytes.
-fI[0] == fI[1] == strtof
+fI[0] = #0 = 0
+fI[1] = #1 = 1.4012985e-45
+fI[0] == strtof
Input: 1.234567890123456789
@@ -188,7 +194,9 @@ strtof consumes 25 bytes and returns 0 = #0
g_ffmt(0) gives 1 bytes: "0"
strtoIf returns 80, consuming 25 bytes.
-fI[0] == fI[1] == strtof
+fI[0] = #0 = 0
+fI[1] = #1 = 1.4012985e-45
+fI[0] == strtof
Input: 1.234567890123456789e-321
@@ -196,7 +204,9 @@ strtof consumes 25 bytes and returns 0 = #0
g_ffmt(0) gives 1 bytes: "0"
strtoIf returns 80, consuming 25 bytes.
-fI[0] == fI[1] == strtof
+fI[0] = #0 = 0
+fI[1] = #1 = 1.4012985e-45
+fI[0] == strtof
Input: 1e23
@@ -224,7 +234,9 @@ strtof consumes 23 bytes and returns 0 = #0
g_ffmt(0) gives 1 bytes: "0"
strtoIf returns 80, consuming 23 bytes.
-fI[0] == fI[1] == strtof
+fI[0] = #0 = 0
+fI[1] = #1 = 1.4012985e-45
+fI[0] == strtof
Input: 9.025971879324147880346310405869e-277
@@ -232,7 +244,9 @@ strtof consumes 37 bytes and returns 0 = #0
g_ffmt(0) gives 1 bytes: "0"
strtoIf returns 80, consuming 37 bytes.
-fI[0] == fI[1] == strtof
+fI[0] = #0 = 0
+fI[1] = #1 = 1.4012985e-45
+fI[0] == strtof
Input: 9.025971879324147880346310405868e-277
@@ -240,7 +254,9 @@ strtof consumes 37 bytes and returns 0 = #0
g_ffmt(0) gives 1 bytes: "0"
strtoIf returns 80, consuming 37 bytes.
-fI[0] == fI[1] == strtof
+fI[0] = #0 = 0
+fI[1] = #1 = 1.4012985e-45
+fI[0] == strtof
Input: 2.2250738585072014e-308
@@ -248,7 +264,9 @@ strtof consumes 23 bytes and returns 0 = #0
g_ffmt(0) gives 1 bytes: "0"
strtoIf returns 80, consuming 23 bytes.
-fI[0] == fI[1] == strtof
+fI[0] = #0 = 0
+fI[1] = #1 = 1.4012985e-45
+fI[0] == strtof
Input: 2.2250738585072013e-308
@@ -256,7 +274,9 @@ strtof consumes 23 bytes and returns 0 = #0
g_ffmt(0) gives 1 bytes: "0"
strtoIf returns 80, consuming 23 bytes.
-fI[0] == fI[1] == strtof
+fI[0] = #0 = 0
+fI[1] = #1 = 1.4012985e-45
+fI[0] == strtof
Rounding mode for strtor... changed from 1 (nearest) to 0 (toward zero)
diff --git a/contrib/gdtoa/test/getround.c b/contrib/gdtoa/test/getround.c
index 1754fb3..786bcc3 100644
--- a/contrib/gdtoa/test/getround.c
+++ b/contrib/gdtoa/test/getround.c
@@ -44,7 +44,14 @@ getround(int r, char *s)
{
int i;
- i = atoi(s+1);
+ while(*++s <= ' ') {
+ if (!*s) {
+ printf("Current round mode for strtor... is %d (%s).\n",
+ r, dir[r]);
+ return r;
+ }
+ }
+ i = atoi(s);
if (i >= 0 && i < 4) {
printf("Rounding mode for strtor... ");
if (i == r)
diff --git a/contrib/gdtoa/test/xsum0.out b/contrib/gdtoa/test/xsum0.out
index 67bbe69..9837ac8 100644
--- a/contrib/gdtoa/test/xsum0.out
+++ b/contrib/gdtoa/test/xsum0.out
@@ -1,11 +1,11 @@
-README e6ebdc91 2429
+README e86e9133 2692
Qtest.c e8353ffc 5046
dItest.c e33800ce 2371
ddtest.c f9d06e7b 4984
dtest.c ee533ac3 4078
dt.c 7eeda57 6384
ftest.c ec8a6654 3999
-getround.c f810968 2011
+getround.c 161043dd 2143
strtoIdSI.c 7bfb88b 49
strtoIddSI.c 72e8852 50
strtodISI.c ed08b740 49
@@ -25,7 +25,7 @@ ddsi.out 1f94bbe2 10251
dd.out e262456e 40923
dtst.out e284ac98 23711
d.out f271efc9 28131
-f.out 4b0bd51 21207
+f.out 9013e91 21537
x.ou0 1402f834 25372
xL.ou0 faa3a741 26363
x.ou1 f1af5a00 34581
diff --git a/contrib/gdtoa/xsum0.out b/contrib/gdtoa/xsum0.out
index 701b55d..5fc30d9 100644
--- a/contrib/gdtoa/xsum0.out
+++ b/contrib/gdtoa/xsum0.out
@@ -1,7 +1,7 @@
-README f2477cff 14190
+README 1d9fab5a 14790
arithchk.c ebbe5bc7 4075
dmisc.c c8daa18 4682
-dtoa.c 6a7b6fe 16876
+dtoa.c 14f5b9e1 17138
g_Qfmt.c f791d807 2839
g__fmt.c 14dca85 2504
g_ddfmt.c 10eae12a 3695
@@ -10,12 +10,12 @@ g_ffmt.c fb83cfb5 2429
g_xLfmt.c f216a096 2686
g_xfmt.c ed824bf3 2775
gdtoa.c e29409a6 16988
-gdtoa.h f208c204 4780
-gdtoaimp.h e3c2a970 19441
-gethex.c dba1616 5201
+gdtoa.h fd46e927 4920
+gdtoaimp.h e753264b 19648
+gethex.c f42e6bdf 6247
gmisc.c 1859d016 2084
hd_init.c efdbe921 1797
-hexnan.c f7ea38f9 2958
+hexnan.c 13f362b7 3462
makefile f890b12 2932
misc.c 1757f7fc 14252
qnan.c efd33d64 3417
@@ -24,12 +24,12 @@ strtoIQ.c 1809dfcf 1939
strtoId.c f41ddac2 1931
strtoIdd.c f13e3bc3 2105
strtoIf.c f12c6af4 1875
-strtoIg.c ef30d392 3454
+strtoIg.c fd8c1bcf 3589
strtoIx.c e50f716d 1960
strtoIxL.c ea0b821b 1931
-strtod.c eec1df60 20532
+strtod.c f6943e52 21001
strtodI.c 1c2440ce 3915
-strtodg.c f6c3dd52 19911
+strtodg.c 1c1bb220 20499
strtodnrp.c af895e9 2538
strtof.c 1c5192d3 2073
strtopQ.c f116d4f0 2563
--
1.7.7