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.
This commit is contained in:
das 2008-09-03 07:23:57 +00:00
parent 88c5e9e192
commit 2732388653
13 changed files with 259 additions and 98 deletions

View File

@ -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 the decimal-point character to be taken from the current locale; otherwise
it is '.'. 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 " Please send comments to David M. Gay (dmg at acm dot org, with " at "
changed at "@" and " dot " changed to "."). changed at "@" and " dot " changed to ".").

View File

@ -66,7 +66,6 @@ THIS SOFTWARE.
*/ */
#ifdef Honor_FLT_ROUNDS #ifdef Honor_FLT_ROUNDS
#define Rounding rounding
#undef Check_FLT_ROUNDS #undef Check_FLT_ROUNDS
#define Check_FLT_ROUNDS #define Check_FLT_ROUNDS
#else #else
@ -127,12 +126,22 @@ dtoa
Bigint *b, *b1, *delta, *mlo, *mhi, *S; Bigint *b, *b1, *delta, *mlo, *mhi, *S;
double d2, ds, eps; double d2, ds, eps;
char *s, *s0; char *s, *s0;
#ifdef Honor_FLT_ROUNDS
int rounding;
#endif
#ifdef SET_INEXACT #ifdef SET_INEXACT
int inexact, oldinexact; int inexact, oldinexact;
#endif #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 #ifndef MULTIPLE_THREADS
if (dtoa_result) { if (dtoa_result) {
@ -178,12 +187,12 @@ dtoa
inexact = 1; inexact = 1;
#endif #endif
#ifdef Honor_FLT_ROUNDS #ifdef Honor_FLT_ROUNDS
if ((rounding = Flt_Rounds) >= 2) { if (Rounding >= 2) {
if (*sign) if (*sign)
rounding = rounding == 2 ? 0 : 2; Rounding = Rounding == 2 ? 0 : 2;
else else
if (rounding != 2) if (Rounding != 2)
rounding = 0; Rounding = 0;
} }
#endif #endif
@ -316,7 +325,7 @@ dtoa
s = s0 = rv_alloc(i); s = s0 = rv_alloc(i);
#ifdef Honor_FLT_ROUNDS #ifdef Honor_FLT_ROUNDS
if (mode > 1 && rounding != 1) if (mode > 1 && Rounding != 1)
leftright = 0; leftright = 0;
#endif #endif
@ -453,7 +462,7 @@ dtoa
if (i == ilim) { if (i == ilim) {
#ifdef Honor_FLT_ROUNDS #ifdef Honor_FLT_ROUNDS
if (mode > 1) if (mode > 1)
switch(rounding) { switch(Rounding) {
case 0: goto ret1; case 0: goto ret1;
case 2: goto bump_up; case 2: goto bump_up;
} }
@ -521,7 +530,7 @@ dtoa
spec_case = 0; spec_case = 0;
if ((mode < 2 || leftright) if ((mode < 2 || leftright)
#ifdef Honor_FLT_ROUNDS #ifdef Honor_FLT_ROUNDS
&& rounding == 1 && Rounding == 1
#endif #endif
) { ) {
if (!word1(d) && !(word0(d) & Bndry_mask) if (!word1(d) && !(word0(d) & Bndry_mask)
@ -614,7 +623,7 @@ dtoa
#ifndef ROUND_BIASED #ifndef ROUND_BIASED
if (j1 == 0 && mode != 1 && !(word1(d) & 1) if (j1 == 0 && mode != 1 && !(word1(d) & 1)
#ifdef Honor_FLT_ROUNDS #ifdef Honor_FLT_ROUNDS
&& rounding >= 1 && Rounding >= 1
#endif #endif
) { ) {
if (dig == '9') if (dig == '9')
@ -642,7 +651,7 @@ dtoa
} }
#ifdef Honor_FLT_ROUNDS #ifdef Honor_FLT_ROUNDS
if (mode > 1) if (mode > 1)
switch(rounding) { switch(Rounding) {
case 0: goto accept_dig; case 0: goto accept_dig;
case 2: goto keep_dig; case 2: goto keep_dig;
} }
@ -660,7 +669,7 @@ dtoa
} }
if (j1 > 0) { if (j1 > 0) {
#ifdef Honor_FLT_ROUNDS #ifdef Honor_FLT_ROUNDS
if (!rounding) if (!Rounding)
goto accept_dig; goto accept_dig;
#endif #endif
if (dig == '9') { /* possible if i == 1 */ if (dig == '9') { /* possible if i == 1 */
@ -703,7 +712,7 @@ dtoa
/* Round off last digit */ /* Round off last digit */
#ifdef Honor_FLT_ROUNDS #ifdef Honor_FLT_ROUNDS
switch(rounding) { switch(Rounding) {
case 0: goto trimzeros; case 0: goto trimzeros;
case 2: goto roundoff; case 2: goto roundoff;
} }

View File

@ -74,9 +74,9 @@ typedef unsigned short UShort;
/* The following may be or-ed into one of the above values. */ /* The following may be or-ed into one of the above values. */
STRTOG_Neg = 0x08, STRTOG_Neg = 0x08, /* does not affect STRTOG_Inexlo or STRTOG_Inexhi */
STRTOG_Inexlo = 0x10, STRTOG_Inexlo = 0x10, /* returned result rounded toward zero */
STRTOG_Inexhi = 0x20, STRTOG_Inexhi = 0x20, /* returned result rounded away from zero */
STRTOG_Inexact = 0x30, STRTOG_Inexact = 0x30,
STRTOG_Underflow= 0x40, STRTOG_Underflow= 0x40,
STRTOG_Overflow = 0x80 STRTOG_Overflow = 0x80

View File

@ -132,13 +132,16 @@ THIS SOFTWARE.
* Infinity and NaN (case insensitively). * Infinity and NaN (case insensitively).
* When INFNAN_CHECK is #defined and No_Hex_NaN is not #defined, * When INFNAN_CHECK is #defined and No_Hex_NaN is not #defined,
* strtodg also accepts (case insensitively) strings of the form * strtodg also accepts (case insensitively) strings of the form
* NaN(x), where x is a string of hexadecimal digits and spaces; * NaN(x), where x is a string of hexadecimal digits (optionally
* if there is only one string of hexadecimal digits, it is taken * preceded by 0x or 0X) and spaces; if there is only one string
* for the fraction bits of the resulting NaN; if there are two or * of hexadecimal digits, it is taken for the fraction bits of the
* more strings of hexadecimal digits, each string is assigned * resulting NaN; if there are two or more strings of hexadecimal
* to the next available sequence of 32-bit words of fractions * digits, each string is assigned to the next available sequence
* bits (starting with the most significant), right-aligned in * of 32-bit words of fractions bits (starting with the most
* each sequence. * 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 * #define MULTIPLE_THREADS if the system offers preemptively scheduled
* multiple threads. In this case, you must provide (or suitably * multiple threads. In this case, you must provide (or suitably
* #define) two locks, acquired by ACQUIRE_DTOA_LOCK(n) and freed * #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. * dtoa. You may do so whether or not MULTIPLE_THREADS is #defined.
* #define IMPRECISE_INEXACT if you do not care about the setting of * #define IMPRECISE_INEXACT if you do not care about the setting of
* the STRTOG_Inexact bits in the special case of doing IEEE double * 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). * dtoa.c).
* #define NO_HEX_FP to disable recognition of C9x's hexadecimal * #define NO_HEX_FP to disable recognition of C9x's hexadecimal
* floating-point constants. * floating-point constants.
@ -204,6 +207,7 @@ extern Char *MALLOC ANSI((size_t));
#define INFNAN_CHECK #define INFNAN_CHECK
#define USE_LOCALE #define USE_LOCALE
#define Honor_FLT_ROUNDS #define Honor_FLT_ROUNDS
#define Trust_FLT_ROUNDS
#undef IEEE_Arith #undef IEEE_Arith
#undef Avoid_Underflow #undef Avoid_Underflow
@ -597,7 +601,7 @@ extern void memcpy_D2A ANSI((void*, const void*, size_t));
extern int cmp ANSI((Bigint*, Bigint*)); extern int cmp ANSI((Bigint*, Bigint*));
extern void copybits ANSI((ULong*, int, Bigint*)); extern void copybits ANSI((ULong*, int, Bigint*));
extern Bigint *d2b ANSI((double, int*, int*)); extern Bigint *d2b ANSI((double, int*, int*));
extern int decrement ANSI((Bigint*)); extern void decrement ANSI((Bigint*));
extern Bigint *diff ANSI((Bigint*, Bigint*)); extern Bigint *diff ANSI((Bigint*, Bigint*));
extern char *dtoa ANSI((double d, int mode, int ndigits, extern char *dtoa ANSI((double d, int mode, int ndigits,
int *decpt, int *sign, char **rve)); int *decpt, int *sign, char **rve));

View File

@ -45,7 +45,7 @@ gethex( CONST char **sp, FPI *fpi, Long *exp, Bigint **bp, int sign)
{ {
Bigint *b; Bigint *b;
CONST unsigned char *decpt, *s0, *s, *s1; 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; ULong L, lostbits, *x;
Long e, e1; Long e, e1;
#ifdef USE_LOCALE #ifdef USE_LOCALE
@ -56,6 +56,7 @@ gethex( CONST char **sp, FPI *fpi, Long *exp, Bigint **bp, int sign)
if (!hexdig['0']) if (!hexdig['0'])
hexdig_init_D2A(); hexdig_init_D2A();
*bp = 0;
havedig = 0; havedig = 0;
s0 = *(CONST unsigned char **)sp + 2; s0 = *(CONST unsigned char **)sp + 2;
while(s0[havedig] == '0') 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); e = -(((Long)(s-decpt)) << 2);
pcheck: pcheck:
s1 = s; s1 = s;
big = esign = 0;
switch(*s) { switch(*s) {
case 'p': case 'p':
case 'P': case 'P':
esign = 0;
switch(*++s) { switch(*++s) {
case '-': case '-':
esign = 1; esign = 1;
@ -106,17 +107,65 @@ gethex( CONST char **sp, FPI *fpi, Long *exp, Bigint **bp, int sign)
break; break;
} }
e1 = n - 0x10; 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; e1 = 10*e1 + n - 0x10;
}
if (esign) if (esign)
e1 = -e1; e1 = -e1;
e += e1; e += e1;
} }
*sp = (char*)s; *sp = (char*)s;
if (zret) {
if (!havedig) if (!havedig)
*sp = s0 - 1; *sp = s0 - 1;
if (zret)
return STRTOG_Zero; 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; n = s1 - s0 - 1;
for(k = 0; n > 7; n >>= 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; k = n - 1;
if (x[k>>kshift] & 1 << (k & kmask)) { if (x[k>>kshift] & 1 << (k & kmask)) {
lostbits = 2; lostbits = 2;
if (k > 1 && any_on(b,k-1)) if (k > 0 && any_on(b,k))
lostbits = 3; lostbits = 3;
} }
} }
@ -165,7 +214,10 @@ gethex( CONST char **sp, FPI *fpi, Long *exp, Bigint **bp, int sign)
if (e > fpi->emax) { if (e > fpi->emax) {
ovfl: ovfl:
Bfree(b); Bfree(b);
*bp = 0; ovfl1:
#ifndef NO_ERRNO
errno = ERANGE;
#endif
return STRTOG_Infinite | STRTOG_Overflow | STRTOG_Inexhi; return STRTOG_Infinite | STRTOG_Overflow | STRTOG_Inexhi;
} }
irv = STRTOG_Normal; irv = STRTOG_Normal;
@ -185,15 +237,22 @@ gethex( CONST char **sp, FPI *fpi, Long *exp, Bigint **bp, int sign)
case FPI_Round_down: case FPI_Round_down:
if (sign) { if (sign) {
one_bit: one_bit:
*exp = fpi->emin;
x[0] = b->wds = 1; x[0] = b->wds = 1;
dret:
*bp = b; *bp = b;
*exp = fpi->emin;
#ifndef NO_ERRNO
errno = ERANGE;
#endif
return STRTOG_Denormal | STRTOG_Inexhi return STRTOG_Denormal | STRTOG_Inexhi
| STRTOG_Underflow; | STRTOG_Underflow;
} }
} }
Bfree(b); Bfree(b);
*bp = 0; retz:
#ifndef NO_ERRNO
errno = ERANGE;
#endif
return STRTOG_Zero | STRTOG_Inexlo | STRTOG_Underflow; return STRTOG_Zero | STRTOG_Inexlo | STRTOG_Underflow;
} }
k = n - 1; k = n - 1;

View File

@ -61,14 +61,18 @@ strtoIg(CONST char *s00, char **se, FPI *fpi, Long *exp, Bigint **B, int *rvp)
if (rv & STRTOG_Inexlo) { if (rv & STRTOG_Inexlo) {
swap = 0; swap = 0;
b1 = increment(b1); b1 = increment(b1);
if (fpi->sudden_underflow if ((rv & STRTOG_Retmask) == STRTOG_Zero) {
&& (rv & STRTOG_Retmask) == STRTOG_Zero) { if (fpi->sudden_underflow) {
b1->x[0] = 0; b1->x[0] = 0;
b1->x[nw1] = 1L << nb11; b1->x[nw1] = 1L << nb11;
rv1 += STRTOG_Normal - STRTOG_Zero; rv1 += STRTOG_Normal - STRTOG_Zero;
rv1 &= ~STRTOG_Underflow; rv1 &= ~STRTOG_Underflow;
goto swapcheck; goto swapcheck;
} }
rv1 &= STRTOG_Inexlo | STRTOG_Underflow | STRTOG_Zero;
rv1 |= STRTOG_Inexhi | STRTOG_Denormal;
goto swapcheck;
}
if (b1->wds > nw if (b1->wds > nw
|| nb1 && b1->x[nw1] & 1L << nb1) { || nb1 && b1->x[nw1] & 1L << nb1) {
if (++e1 > fpi->emax) if (++e1 > fpi->emax)

View File

@ -44,16 +44,15 @@ THIS SOFTWARE.
#ifndef NO_IEEE_Scale #ifndef NO_IEEE_Scale
#define Avoid_Underflow #define Avoid_Underflow
#undef tinytens #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. */ /* 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, static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,
9007199254740992.e-256 9007199254740992.*9007199254740992.e-256
}; };
#endif #endif
#endif #endif
#ifdef Honor_FLT_ROUNDS #ifdef Honor_FLT_ROUNDS
#define Rounding rounding
#undef Check_FLT_ROUNDS #undef Check_FLT_ROUNDS
#define Check_FLT_ROUNDS #define Check_FLT_ROUNDS
#else #else
@ -81,9 +80,19 @@ strtod
#ifdef SET_INEXACT #ifdef SET_INEXACT
int inexact, oldinexact; int inexact, oldinexact;
#endif #endif
#ifdef Honor_FLT_ROUNDS #ifdef Honor_FLT_ROUNDS /*{*/
int rounding; int Rounding;
#endif #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; sign = nz0 = nz = decpt = 0;
dval(rv) = 0.; dval(rv) = 0.;
@ -109,7 +118,7 @@ strtod
} }
break2: break2:
if (*s == '0') { 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 }; static FPI fpi = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI };
Long exp; Long exp;
@ -118,16 +127,20 @@ strtod
case 'x': case 'x':
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; FPI fpi1 = fpi;
#ifdef Honor_FLT_ROUNDS /*{{*/
fpi1.rounding = Rounding;
#else /*}{*/
switch(fegetround()) { switch(fegetround()) {
case FE_TOWARDZERO: fpi1.rounding = 0; break; case FE_TOWARDZERO: fpi1.rounding = 0; break;
case FE_UPWARD: fpi1.rounding = 2; break; case FE_UPWARD: fpi1.rounding = 2; break;
case FE_DOWNWARD: fpi1.rounding = 3; case FE_DOWNWARD: fpi1.rounding = 3;
} }
#else #endif /*}}*/
#else /*}{*/
#define fpi1 fpi #define fpi1 fpi
#endif #endif /*}}*/
switch((i = gethex(&s, &fpi1, &exp, &bb, sign)) & STRTOG_Retmask) { switch((i = gethex(&s, &fpi1, &exp, &bb, sign)) & STRTOG_Retmask) {
case STRTOG_NoNumber: case STRTOG_NoNumber:
s = s00; s = s00;
@ -381,12 +394,12 @@ strtod
scale = 0; scale = 0;
#endif #endif
#ifdef Honor_FLT_ROUNDS #ifdef Honor_FLT_ROUNDS
if ((rounding = Flt_Rounds) >= 2) { if (Rounding >= 2) {
if (sign) if (sign)
rounding = rounding == 2 ? 0 : 2; Rounding = Rounding == 2 ? 0 : 2;
else else
if (rounding != 2) if (Rounding != 2)
rounding = 0; Rounding = 0;
} }
#endif #endif
#endif /*IEEE_Arith*/ #endif /*IEEE_Arith*/
@ -405,7 +418,7 @@ strtod
/* Can't trust HUGE_VAL */ /* Can't trust HUGE_VAL */
#ifdef IEEE_Arith #ifdef IEEE_Arith
#ifdef Honor_FLT_ROUNDS #ifdef Honor_FLT_ROUNDS
switch(rounding) { switch(Rounding) {
case 0: /* toward 0 */ case 0: /* toward 0 */
case 3: /* toward -infinity */ case 3: /* toward -infinity */
word0(rv) = Big0; word0(rv) = Big0;
@ -536,7 +549,7 @@ strtod
bd2 -= bbe; bd2 -= bbe;
bs2 = bb2; bs2 = bb2;
#ifdef Honor_FLT_ROUNDS #ifdef Honor_FLT_ROUNDS
if (rounding != 1) if (Rounding != 1)
bs2++; bs2++;
#endif #endif
#ifdef Avoid_Underflow #ifdef Avoid_Underflow
@ -594,7 +607,7 @@ strtod
delta->sign = 0; delta->sign = 0;
i = cmp(delta, bs); i = cmp(delta, bs);
#ifdef Honor_FLT_ROUNDS #ifdef Honor_FLT_ROUNDS
if (rounding != 1) { if (Rounding != 1) {
if (i < 0) { if (i < 0) {
/* Error is less than an ulp */ /* Error is less than an ulp */
if (!delta->x[0] && delta->wds <= 1) { if (!delta->x[0] && delta->wds <= 1) {
@ -604,7 +617,7 @@ strtod
#endif #endif
break; break;
} }
if (rounding) { if (Rounding) {
if (dsign) { if (dsign) {
adj = 1.; adj = 1.;
goto apply_adj; goto apply_adj;
@ -650,10 +663,10 @@ strtod
if (adj < 1.) if (adj < 1.)
adj = 1.; adj = 1.;
if (adj <= 0x7ffffffe) { if (adj <= 0x7ffffffe) {
/* adj = rounding ? ceil(adj) : floor(adj); */ /* adj = Rounding ? ceil(adj) : floor(adj); */
y = adj; y = adj;
if (y != adj) { if (y != adj) {
if (!((rounding>>1) ^ dsign)) if (!((Rounding>>1) ^ dsign))
y++; y++;
adj = y; adj = y;
} }
@ -676,8 +689,11 @@ strtod
#endif /*Sudden_Underflow*/ #endif /*Sudden_Underflow*/
#endif /*Avoid_Underflow*/ #endif /*Avoid_Underflow*/
adj *= ulp(dval(rv)); adj *= ulp(dval(rv));
if (dsign) if (dsign) {
if (word0(rv) == Big0 && word1(rv) == Big1)
goto ovfl;
dval(rv) += adj; dval(rv) += adj;
}
else else
dval(rv) -= adj; dval(rv) -= adj;
goto cont; goto cont;
@ -770,7 +786,7 @@ strtod
} }
#endif /*Avoid_Underflow*/ #endif /*Avoid_Underflow*/
L = (word0(rv) & Exp_mask) - Exp_msk1; L = (word0(rv) & Exp_mask) - Exp_msk1;
#endif /*Sudden_Underflow}*/ #endif /*Sudden_Underflow}}*/
word0(rv) = L | Bndry_mask1; word0(rv) = L | Bndry_mask1;
word1(rv) = 0xffffffff; word1(rv) = 0xffffffff;
#ifdef IBM #ifdef IBM

View File

@ -89,7 +89,7 @@ increment(Bigint *b)
return b; return b;
} }
int void
#ifdef KR_headers #ifdef KR_headers
decrement(b) Bigint *b; decrement(b) Bigint *b;
#else #else
@ -119,7 +119,6 @@ decrement(Bigint *b)
*x++ = y & 0xffff; *x++ = y & 0xffff;
} while(borrow && x < xe); } while(borrow && x < xe);
#endif #endif
return STRTOG_Inexlo;
} }
static int static int
@ -206,9 +205,9 @@ rvOK
goto ret; goto ret;
} }
switch(rd) { switch(rd) {
case 1: case 1: /* round down (toward -Infinity) */
goto trunc; goto trunc;
case 2: case 2: /* round up (toward +Infinity) */
break; break;
default: /* round near */ default: /* round near */
k = bdif - 1; k = bdif - 1;
@ -330,7 +329,7 @@ strtodg
CONST char *s, *s0, *s1; CONST char *s, *s0, *s1;
double adj, adj0, rv, tol; double adj, adj0, rv, tol;
Long L; Long L;
ULong y, z; ULong *b, *be, y, z;
Bigint *ab, *bb, *bb1, *bd, *bd0, *bs, *delta, *rvb, *rvb0; Bigint *ab, *bb, *bb1, *bd, *bd0, *bs, *delta, *rvb, *rvb0;
irv = STRTOG_Zero; irv = STRTOG_Zero;
@ -822,10 +821,8 @@ strtodg
break; break;
if (dsign) { if (dsign) {
rvb = increment(rvb); rvb = increment(rvb);
if ( (j = rvbits & kmask) !=0) j = kmask & (ULbits - (rvbits & kmask));
j = ULbits - j; if (hi0bits(rvb->x[rvb->wds - 1]) != j)
if (hi0bits(rvb->x[(rvb->wds - 1) >> kshift])
!= j)
rvbits++; rvbits++;
irv = STRTOG_Normal | STRTOG_Inexhi; irv = STRTOG_Normal | STRTOG_Inexhi;
} }
@ -978,6 +975,29 @@ strtodg
Bfree(bd0); Bfree(bd0);
Bfree(delta); Bfree(delta);
if (rve > fpi->emax) { 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: huge:
rvb->wds = 0; rvb->wds = 0;
irv = STRTOG_Infinite | STRTOG_Overflow | STRTOG_Inexhi; irv = STRTOG_Infinite | STRTOG_Overflow | STRTOG_Inexhi;
@ -992,12 +1012,19 @@ strtodg
if (sudden_underflow) { if (sudden_underflow) {
rvb->wds = 0; rvb->wds = 0;
irv = STRTOG_Underflow | STRTOG_Inexlo; irv = STRTOG_Underflow | STRTOG_Inexlo;
#ifndef NO_ERRNO
errno = ERANGE;
#endif
} }
else { else {
irv = (irv & ~STRTOG_Retmask) | irv = (irv & ~STRTOG_Retmask) |
(rvb->wds > 0 ? STRTOG_Denormal : STRTOG_Zero); (rvb->wds > 0 ? STRTOG_Denormal : STRTOG_Zero);
if (irv & STRTOG_Inexact) if (irv & STRTOG_Inexact) {
irv |= STRTOG_Underflow; irv |= STRTOG_Underflow;
#ifndef NO_ERRNO
errno = ERANGE;
#endif
}
} }
} }
if (se) if (se)

View File

@ -55,7 +55,12 @@ logic (on double and double-double conversions).
Program strtodt tests strtod on some hard cases (in file testnos3) Program strtodt tests strtod on some hard cases (in file testnos3)
posted by Fred Tydeman to comp.arch.arithmetic on 26 Feb. 1996. 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, These are simple test programs, not meant for exhaustive testing,
but for manually testing "interesting" cases. Paxson's testbase but for manually testing "interesting" cases. Paxson's testbase
is good for more exhaustive testing, in part with random inputs. is good for more exhaustive testing, in part with random inputs.
See ftp://ftp.ee.lbl.gov/testbase-report.ps.Z .

View File

@ -124,7 +124,9 @@ strtof consumes 9 bytes and returns 0 = #0
g_ffmt(0) gives 1 bytes: "0" g_ffmt(0) gives 1 bytes: "0"
strtoIf returns 80, consuming 9 bytes. 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 Input: 1.23e-320
@ -132,7 +134,9 @@ strtof consumes 9 bytes and returns 0 = #0
g_ffmt(0) gives 1 bytes: "0" g_ffmt(0) gives 1 bytes: "0"
strtoIf returns 80, consuming 9 bytes. 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 Input: 1.23e-20
@ -160,7 +164,9 @@ strtof consumes 15 bytes and returns 0 = #0
g_ffmt(0) gives 1 bytes: "0" g_ffmt(0) gives 1 bytes: "0"
strtoIf returns 80, consuming 15 bytes. 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 Input: 1.234567890123456789
@ -188,7 +194,9 @@ strtof consumes 25 bytes and returns 0 = #0
g_ffmt(0) gives 1 bytes: "0" g_ffmt(0) gives 1 bytes: "0"
strtoIf returns 80, consuming 25 bytes. 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 Input: 1.234567890123456789e-321
@ -196,7 +204,9 @@ strtof consumes 25 bytes and returns 0 = #0
g_ffmt(0) gives 1 bytes: "0" g_ffmt(0) gives 1 bytes: "0"
strtoIf returns 80, consuming 25 bytes. 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 Input: 1e23
@ -224,7 +234,9 @@ strtof consumes 23 bytes and returns 0 = #0
g_ffmt(0) gives 1 bytes: "0" g_ffmt(0) gives 1 bytes: "0"
strtoIf returns 80, consuming 23 bytes. 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 Input: 9.025971879324147880346310405869e-277
@ -232,7 +244,9 @@ strtof consumes 37 bytes and returns 0 = #0
g_ffmt(0) gives 1 bytes: "0" g_ffmt(0) gives 1 bytes: "0"
strtoIf returns 80, consuming 37 bytes. 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 Input: 9.025971879324147880346310405868e-277
@ -240,7 +254,9 @@ strtof consumes 37 bytes and returns 0 = #0
g_ffmt(0) gives 1 bytes: "0" g_ffmt(0) gives 1 bytes: "0"
strtoIf returns 80, consuming 37 bytes. 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 Input: 2.2250738585072014e-308
@ -248,7 +264,9 @@ strtof consumes 23 bytes and returns 0 = #0
g_ffmt(0) gives 1 bytes: "0" g_ffmt(0) gives 1 bytes: "0"
strtoIf returns 80, consuming 23 bytes. 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 Input: 2.2250738585072013e-308
@ -256,7 +274,9 @@ strtof consumes 23 bytes and returns 0 = #0
g_ffmt(0) gives 1 bytes: "0" g_ffmt(0) gives 1 bytes: "0"
strtoIf returns 80, consuming 23 bytes. 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) Rounding mode for strtor... changed from 1 (nearest) to 0 (toward zero)

View File

@ -44,7 +44,14 @@ getround(int r, char *s)
{ {
int i; 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) { if (i >= 0 && i < 4) {
printf("Rounding mode for strtor... "); printf("Rounding mode for strtor... ");
if (i == r) if (i == r)

View File

@ -1,11 +1,11 @@
README e6ebdc91 2429 README e86e9133 2692
Qtest.c e8353ffc 5046 Qtest.c e8353ffc 5046
dItest.c e33800ce 2371 dItest.c e33800ce 2371
ddtest.c f9d06e7b 4984 ddtest.c f9d06e7b 4984
dtest.c ee533ac3 4078 dtest.c ee533ac3 4078
dt.c 7eeda57 6384 dt.c 7eeda57 6384
ftest.c ec8a6654 3999 ftest.c ec8a6654 3999
getround.c f810968 2011 getround.c 161043dd 2143
strtoIdSI.c 7bfb88b 49 strtoIdSI.c 7bfb88b 49
strtoIddSI.c 72e8852 50 strtoIddSI.c 72e8852 50
strtodISI.c ed08b740 49 strtodISI.c ed08b740 49
@ -25,7 +25,7 @@ ddsi.out 1f94bbe2 10251
dd.out e262456e 40923 dd.out e262456e 40923
dtst.out e284ac98 23711 dtst.out e284ac98 23711
d.out f271efc9 28131 d.out f271efc9 28131
f.out 4b0bd51 21207 f.out 9013e91 21537
x.ou0 1402f834 25372 x.ou0 1402f834 25372
xL.ou0 faa3a741 26363 xL.ou0 faa3a741 26363
x.ou1 f1af5a00 34581 x.ou1 f1af5a00 34581

View File

@ -1,7 +1,7 @@
README f2477cff 14190 README 1d9fab5a 14790
arithchk.c ebbe5bc7 4075 arithchk.c ebbe5bc7 4075
dmisc.c c8daa18 4682 dmisc.c c8daa18 4682
dtoa.c 6a7b6fe 16876 dtoa.c 14f5b9e1 17138
g_Qfmt.c f791d807 2839 g_Qfmt.c f791d807 2839
g__fmt.c 14dca85 2504 g__fmt.c 14dca85 2504
g_ddfmt.c 10eae12a 3695 g_ddfmt.c 10eae12a 3695
@ -10,12 +10,12 @@ g_ffmt.c fb83cfb5 2429
g_xLfmt.c f216a096 2686 g_xLfmt.c f216a096 2686
g_xfmt.c ed824bf3 2775 g_xfmt.c ed824bf3 2775
gdtoa.c e29409a6 16988 gdtoa.c e29409a6 16988
gdtoa.h f208c204 4780 gdtoa.h fd46e927 4920
gdtoaimp.h e3c2a970 19441 gdtoaimp.h e753264b 19648
gethex.c dba1616 5201 gethex.c f42e6bdf 6247
gmisc.c 1859d016 2084 gmisc.c 1859d016 2084
hd_init.c efdbe921 1797 hd_init.c efdbe921 1797
hexnan.c f7ea38f9 2958 hexnan.c 13f362b7 3462
makefile f890b12 2932 makefile f890b12 2932
misc.c 1757f7fc 14252 misc.c 1757f7fc 14252
qnan.c efd33d64 3417 qnan.c efd33d64 3417
@ -24,12 +24,12 @@ strtoIQ.c 1809dfcf 1939
strtoId.c f41ddac2 1931 strtoId.c f41ddac2 1931
strtoIdd.c f13e3bc3 2105 strtoIdd.c f13e3bc3 2105
strtoIf.c f12c6af4 1875 strtoIf.c f12c6af4 1875
strtoIg.c ef30d392 3454 strtoIg.c fd8c1bcf 3589
strtoIx.c e50f716d 1960 strtoIx.c e50f716d 1960
strtoIxL.c ea0b821b 1931 strtoIxL.c ea0b821b 1931
strtod.c eec1df60 20532 strtod.c f6943e52 21001
strtodI.c 1c2440ce 3915 strtodI.c 1c2440ce 3915
strtodg.c f6c3dd52 19911 strtodg.c 1c1bb220 20499
strtodnrp.c af895e9 2538 strtodnrp.c af895e9 2538
strtof.c 1c5192d3 2073 strtof.c 1c5192d3 2073
strtopQ.c f116d4f0 2563 strtopQ.c f116d4f0 2563