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.

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