- Updated gdtoa in order to fix clang warnings.

gdtoa/README

@@ -56,7 +56,9 @@ two letters:
 		whose sum is the desired value
 
 For decimal -> binary conversions, there are three families of
-helper routines: one for round-nearest:
+helper routines: one for round-nearest (or the current rounding
+mode on IEEE-arithmetic systems that provide the C99 fegetround()
+function, if compiled with -DHonor_FLT_ROUNDS):
 
 	strtof
 	strtod
@@ -150,12 +152,14 @@ suffer double rounding due to use of extended-precision registers.
 For some conversions this variant of strtod is less efficient than the
 one in strtod.c when the latter is run with 53-bit rounding precision.
 
-The values that the strto* routines return for NaNs are determined by
+When float or double are involved, the values that the strto* routines
-gd_qnan.h, which the makefile generates by running the program whose
+return for NaNs are determined by gd_qnan.h, which the makefile
-source is qnan.c.  Note that the rules for distinguishing signaling
+generates by running the program whose source is qnan.c.  For other
-from quiet NaNs are system-dependent.  For cross-compilation, you need
+types, default NaN values are specified in g__fmt.c and may need
-to determine arith.h and gd_qnan.h suitably, e.g., using the
+adjusting.  Note that the rules for distinguishing signaling from
-arithmetic of the target machine.
+quiet NaNs are system-dependent.  For cross-compilation, you need to
+determine arith.h and gd_qnan.h suitably, e.g., using the arithmetic
+of the target machine.
 
 C99's hexadecimal floating-point constants are recognized by the
 strto* routines (but this feature has not yet been heavily tested).
@@ -170,10 +174,11 @@ 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.
+each sequence.  Strings of hexadecimal digits may be preceded by "0x"
+or "0X".
 
-For binary -> decimal conversions, I've provided just one family
+For binary -> decimal conversions, I've provided a family of helper
-of helper routines:
+routines:
 
 	g_ffmt
 	g_dfmt
@@ -181,6 +186,12 @@ of helper routines:
 	g_xfmt
 	g_xLfmt
 	g_Qfmt
+	g_ffmt_p
+	g_dfmt_p
+	g_ddfmt_p
+	g_xfmt_p
+	g_xLfmt_p
+	g_Qfmt_p
 
 which do a "%g" style conversion either to a specified number of decimal
 places (if their ndig argument is positive), or to the shortest
@@ -191,6 +202,36 @@ in the buffer, if the buffer was long enough, or 0.  Other forms of
 conversion are easily done with the help of gdtoa(), such as %e or %f
 style and conversions with direction of rounding specified (so that, if
 desired, the decimal value is either >= or <= the binary value).
+On IEEE-arithmetic systems that provide the C99 fegetround() function,
+if compiled with -DHonor_FLT_ROUNDS, these routines honor the current
+rounding mode.  For pedants, the ...fmt_p() routines are similar to the
+...fmt() routines, but have an additional final int argument, nik,
+that for conversions of Infinity or NaN, determines whether upper,
+lower, or mixed case is used, whether (...) is added to NaN values,
+and whether the sign of a NaN is reported or suppressed:
+
+	nik = ic + 6*(nb + 3*ns),
+
+where ic with 0 <= ic < 6 controls the rendering of Infinity and NaN:
+
+	0 ==> Infinity or NaN
+	1 ==> infinity or nan
+	2 ==> INFINITY or NAN
+	3 ==> Inf or NaN
+	4 ==> inf or nan
+	5 ==> INF or NAN
+
+nb with 0 <= nb < 3 determines whether NaN values are rendered
+as NaN(...):
+
+	0 ==> no
+	1 ==> yes
+	2 ==> no for default NaN values; yes otherwise
+
+ns = 0 or 1 determines whether the sign of NaN values reported:
+
+	0 ==> distinguish NaN and -NaN
+	1 ==> report both as NaN
 
 For an example of more general conversions based on dtoa(), see
 netlib's "printf.c from ampl/solvers".
@@ -332,5 +373,28 @@ 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.
+
+Compile with -DUSE_LOCALE to use the current locale; otherwise
+decimal points are assumed to be '.'.  With -DUSE_LOCALE, unless
+you also compile with -DNO_LOCALE_CACHE, the details about the
+current "decimal point" character string are cached and assumed not
+to change during the program's execution.
+
+On machines with a 64-bit long double and perhaps a 113-bit "quad"
+type, you can invoke "make Printf" to add Printf (and variants, such
+as Fprintf) to gdtoa.a.  These are analogs, declared in stdio1.h, of
+printf and fprintf, etc. in which %La, %Le, %Lf, and %Lg are for long
+double and (if appropriate) %Lqa, %Lqe, %Lqf, and %Lqg are for quad
+precision printing.
+
 Please send comments to	David M. Gay (dmg at acm dot org, with " at "
 changed at "@" and " dot " changed to ".").

gdtoa/arithchk.c

@@ -106,7 +106,7 @@ ccheck()
 	long Cray1;
 
 	/* Cray1 = 4617762693716115456 -- without overflow on non-Crays */
-	Cray1 = printf(emptyfmt) < 0 ? 0 : 4617762;
+	Cray1 = printf("%s", emptyfmt) < 0 ? 0 : 4617762;
 	if (printf(emptyfmt, Cray1) >= 0)
 		Cray1 = 1000000*Cray1 + 693716;
 	if (printf(emptyfmt, Cray1) >= 0)

gdtoa/dmisc.c

@@ -46,7 +46,7 @@ rv_alloc(int i)
 
 	j = sizeof(ULong);
 	for(k = 0;
-		sizeof(Bigint) - sizeof(ULong) - sizeof(int) + j <= (size_t)i;
+		sizeof(Bigint) - sizeof(ULong) - sizeof(int) + j <= (size_t)(i);
 		j <<= 1)
 			k++;
 	r = (int*)Balloc(k);

gdtoa/dtoa.c

@@ -30,7 +30,6 @@ THIS SOFTWARE.
  * with " at " changed at "@" and " dot " changed to ".").	*/
 
 #include "gdtoaimp.h"
-#include <limits.h>
 
 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
  *
@@ -67,7 +66,6 @@ THIS SOFTWARE.
  */
 
 #ifdef Honor_FLT_ROUNDS
-#define Rounding rounding
 #undef Check_FLT_ROUNDS
 #define Check_FLT_ROUNDS
 #else
@@ -77,17 +75,17 @@ THIS SOFTWARE.
  char *
 dtoa
 #ifdef KR_headers
-	(d, mode, ndigits, decpt, sign, rve)
+	(d0, mode, ndigits, decpt, sign, rve)
-	double d; int mode, ndigits, *decpt, *sign; char **rve;
+	double d0; int mode, ndigits, *decpt, *sign; char **rve;
 #else
-	(double _d, int mode, int ndigits, int *decpt, int *sign, char **rve)
+	(double d0, int mode, int ndigits, int *decpt, int *sign, char **rve)
 #endif
 {
  /*	Arguments ndigits, decpt, sign are similar to those
 	of ecvt and fcvt; trailing zeros are suppressed from
 	the returned string.  If not null, *rve is set to point
 	to the end of the return value.  If d is +-Infinity or NaN,
-	then *decpt is set to INT_MAX.
+	then *decpt is set to 9999.
 
 	mode:
 		0 ==> shortest string that yields d when read in
@@ -129,12 +127,22 @@ dtoa
 	U d, d2, eps;
 	double ds;
 	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) {
@@ -142,36 +150,35 @@ dtoa
 		dtoa_result = 0;
 		}
 #endif
-
+	d.d = d0;
-	dval(d) = _d;
+	if (word0(&d) & Sign_bit) {
-	if (word0(d) & Sign_bit) {
 		/* set sign for everything, including 0's and NaNs */
 		*sign = 1;
-		word0(d) &= ~Sign_bit;	/* clear sign bit */
+		word0(&d) &= ~Sign_bit;	/* clear sign bit */
 		}
 	else
 		*sign = 0;
 
 #if defined(IEEE_Arith) + defined(VAX)
 #ifdef IEEE_Arith
-	if ((word0(d) & Exp_mask) == Exp_mask)
+	if ((word0(&d) & Exp_mask) == Exp_mask)
 #else
-	if (word0(d)  == 0x8000)
+	if (word0(&d)  == 0x8000)
 #endif
 		{
 		/* Infinity or NaN */
-		*decpt = INT_MAX;
+		*decpt = 9999;
 #ifdef IEEE_Arith
-		if (!word1(d) && !(word0(d) & 0xfffff))
+		if (!word1(&d) && !(word0(&d) & 0xfffff))
 			return nrv_alloc("Infinity", rve, 8);
 #endif
 		return nrv_alloc("NaN", rve, 3);
 		}
 #endif
 #ifdef IBM
-	dval(d) += 0; /* normalize */
+	dval(&d) += 0; /* normalize */
 #endif
-	if (!dval(d)) {
+	if (!dval(&d)) {
 		*decpt = 1;
 		return nrv_alloc("0", rve, 1);
 		}
@@ -181,35 +188,35 @@ 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)
+			if (Rounding != 2)
-				rounding = 0;
+				Rounding = 0;
 		}
 #endif
 
-	b = d2b(dval(d), &be, &bbits);
+	b = d2b(dval(&d), &be, &bbits);
 #ifdef Sudden_Underflow
-	i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
+	i = (int)(word0(&d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
 #else
-	if (( i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1)) )!=0) {
+	if (( i = (int)(word0(&d) >> Exp_shift1 & (Exp_mask>>Exp_shift1)) )!=0) {
 #endif
-		dval(d2) = dval(d);
+		dval(&d2) = dval(&d);
-		word0(d2) &= Frac_mask1;
+		word0(&d2) &= Frac_mask1;
-		word0(d2) |= Exp_11;
+		word0(&d2) |= Exp_11;
 #ifdef IBM
-		if (( j = 11 - hi0bits(word0(d2) & Frac_mask) )!=0)
+		if (( j = 11 - hi0bits(word0(&d2) & Frac_mask) )!=0)
-			dval(d2) /= 1 << j;
+			dval(&d2) /= 1 << j;
 #endif
 
 		/* log(x)	~=~ log(1.5) + (x-1.5)/1.5
 		 * log10(x)	 =  log(x) / log(10)
 		 *		~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
-		 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
+		 * log10(&d) = (i-Bias)*log(2)/log(10) + log10(&d2)
 		 *
-		 * This suggests computing an approximation k to log10(d) by
+		 * This suggests computing an approximation k to log10(&d) by
 		 *
 		 * k = (i - Bias)*0.301029995663981
 		 *	+ ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
@@ -238,21 +245,21 @@ dtoa
 		/* d is denormalized */
 
 		i = bbits + be + (Bias + (P-1) - 1);
-		x = i > 32  ? word0(d) << 64 - i | word1(d) >> i - 32
+		x = i > 32  ? word0(&d) << (64 - i) | word1(&d) >> (i - 32)
-			    : word1(d) << 32 - i;
+			    : word1(&d) << (32 - i);
-		dval(d2) = x;
+		dval(&d2) = x;
-		word0(d2) -= 31*Exp_msk1; /* adjust exponent */
+		word0(&d2) -= 31*Exp_msk1; /* adjust exponent */
 		i -= (Bias + (P-1) - 1) + 1;
 		denorm = 1;
 		}
 #endif
-	ds = (dval(d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
+	ds = (dval(&d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
 	k = (int)ds;
 	if (ds < 0. && ds != k)
 		k--;	/* want k = floor(ds) */
 	k_check = 1;
 	if (k >= 0 && k <= Ten_pmax) {
-		if (dval(d) < tens[k])
+		if (dval(&d) < tens[k])
 			k--;
 		k_check = 0;
 		}
@@ -291,10 +298,11 @@ dtoa
 		try_quick = 0;
 		}
 	leftright = 1;
+	ilim = ilim1 = -1;	/* Values for cases 0 and 1; done here to */
+				/* silence erroneous "gcc -Wall" warning. */
 	switch(mode) {
 		case 0:
 		case 1:
-			ilim = ilim1 = -1;
 			i = 18;
 			ndigits = 0;
 			break;
@@ -319,7 +327,7 @@ dtoa
 	s = s0 = rv_alloc(i);
 
 #ifdef Honor_FLT_ROUNDS
-	if (mode > 1 && rounding != 1)
+	if (mode > 1 && Rounding != 1)
 		leftright = 0;
 #endif
 
@@ -328,7 +336,7 @@ dtoa
 		/* Try to get by with floating-point arithmetic. */
 
 		i = 0;
-		dval(d2) = dval(d);
+		dval(&d2) = dval(&d);
 		k0 = k;
 		ilim0 = ilim;
 		ieps = 2; /* conservative */
@@ -338,7 +346,7 @@ dtoa
 			if (j & Bletch) {
 				/* prevent overflows */
 				j &= Bletch - 1;
-				dval(d) /= bigtens[n_bigtens-1];
+				dval(&d) /= bigtens[n_bigtens-1];
 				ieps++;
 				}
 			for(; j; j >>= 1, i++)
@@ -346,32 +354,32 @@ dtoa
 					ieps++;
 					ds *= bigtens[i];
 					}
-			dval(d) /= ds;
+			dval(&d) /= ds;
 			}
 		else if (( j1 = -k )!=0) {
-			dval(d) *= tens[j1 & 0xf];
+			dval(&d) *= tens[j1 & 0xf];
 			for(j = j1 >> 4; j; j >>= 1, i++)
 				if (j & 1) {
 					ieps++;
-					dval(d) *= bigtens[i];
+					dval(&d) *= bigtens[i];
 					}
 			}
-		if (k_check && dval(d) < 1. && ilim > 0) {
+		if (k_check && dval(&d) < 1. && ilim > 0) {
 			if (ilim1 <= 0)
 				goto fast_failed;
 			ilim = ilim1;
 			k--;
-			dval(d) *= 10.;
+			dval(&d) *= 10.;
 			ieps++;
 			}
-		dval(eps) = ieps*dval(d) + 7.;
+		dval(&eps) = ieps*dval(&d) + 7.;
-		word0(eps) -= (P-1)*Exp_msk1;
+		word0(&eps) -= (P-1)*Exp_msk1;
 		if (ilim == 0) {
 			S = mhi = 0;
-			dval(d) -= 5.;
+			dval(&d) -= 5.;
-			if (dval(d) > dval(eps))
+			if (dval(&d) > dval(&eps))
 				goto one_digit;
-			if (dval(d) < -dval(eps))
+			if (dval(&d) < -dval(&eps))
 				goto no_digits;
 			goto fast_failed;
 			}
@@ -380,34 +388,34 @@ dtoa
 			/* Use Steele & White method of only
 			 * generating digits needed.
 			 */
-			dval(eps) = 0.5/tens[ilim-1] - dval(eps);
+			dval(&eps) = 0.5/tens[ilim-1] - dval(&eps);
 			for(i = 0;;) {
-				L = (Long)dval(d);
+				L = (Long)dval(&d);
-				dval(d) -= L;
+				dval(&d) -= L;
 				*s++ = '0' + (int)L;
-				if (dval(d) < dval(eps))
+				if (dval(&d) < dval(&eps))
 					goto ret1;
-				if (1. - dval(d) < dval(eps))
+				if (1. - dval(&d) < dval(&eps))
 					goto bump_up;
 				if (++i >= ilim)
 					break;
-				dval(eps) *= 10.;
+				dval(&eps) *= 10.;
-				dval(d) *= 10.;
+				dval(&d) *= 10.;
 				}
 			}
 		else {
 #endif
 			/* Generate ilim digits, then fix them up. */
-			dval(eps) *= tens[ilim-1];
+			dval(&eps) *= tens[ilim-1];
-			for(i = 1;; i++, dval(d) *= 10.) {
+			for(i = 1;; i++, dval(&d) *= 10.) {
-				L = (Long)(dval(d));
+				L = (Long)(dval(&d));
-				if (!(dval(d) -= L))
+				if (!(dval(&d) -= L))
 					ilim = i;
 				*s++ = '0' + (int)L;
 				if (i == ilim) {
-					if (dval(d) > 0.5 + dval(eps))
+					if (dval(&d) > 0.5 + dval(&eps))
 						goto bump_up;
-					else if (dval(d) < 0.5 - dval(eps)) {
+					else if (dval(&d) < 0.5 - dval(&eps)) {
 						while(*--s == '0');
 						s++;
 						goto ret1;
@@ -420,7 +428,7 @@ dtoa
 #endif
  fast_failed:
 		s = s0;
-		dval(d) = dval(d2);
+		dval(&d) = dval(&d2);
 		k = k0;
 		ilim = ilim0;
 		}
@@ -432,22 +440,22 @@ dtoa
 		ds = tens[k];
 		if (ndigits < 0 && ilim <= 0) {
 			S = mhi = 0;
-			if (ilim < 0 || dval(d) <= 5*ds)
+			if (ilim < 0 || dval(&d) <= 5*ds)
 				goto no_digits;
 			goto one_digit;
 			}
-		for(i = 1;; i++, dval(d) *= 10.) {
+		for(i = 1;; i++, dval(&d) *= 10.) {
-			L = (Long)(dval(d) / ds);
+			L = (Long)(dval(&d) / ds);
-			dval(d) -= L*ds;
+			dval(&d) -= L*ds;
 #ifdef Check_FLT_ROUNDS
 			/* If FLT_ROUNDS == 2, L will usually be high by 1 */
-			if (dval(d) < 0) {
+			if (dval(&d) < 0) {
 				L--;
-				dval(d) += ds;
+				dval(&d) += ds;
 				}
 #endif
 			*s++ = '0' + (int)L;
-			if (!dval(d)) {
+			if (!dval(&d)) {
 #ifdef SET_INEXACT
 				inexact = 0;
 #endif
@@ -456,13 +464,18 @@ dtoa
 			if (i == ilim) {
 #ifdef Honor_FLT_ROUNDS
 				if (mode > 1)
-				switch(rounding) {
+				switch(Rounding) {
 				  case 0: goto ret1;
 				  case 2: goto bump_up;
 				  }
 #endif
-				dval(d) += dval(d);
+				dval(&d) += dval(&d);
-				if (dval(d) > ds || dval(d) == ds && L & 1) {
+#ifdef ROUND_BIASED
+				if (dval(&d) >= ds)
+#else
+				if (dval(&d) > ds || (dval(&d) == ds && L & 1))
+#endif
+					{
  bump_up:
 					while(*--s == '9')
 						if (s == s0) {
@@ -524,12 +537,12 @@ dtoa
 	spec_case = 0;
 	if ((mode < 2 || leftright)
 #ifdef Honor_FLT_ROUNDS
-			&& rounding == 1
+			&& Rounding == 1
 #endif
 				) {
-		if (!word1(d) && !(word0(d) & Bndry_mask)
+		if (!word1(&d) && !(word0(&d) & Bndry_mask)
 #ifndef Sudden_Underflow
-		 && word0(d) & (Exp_mask & ~Exp_msk1)
+		 && word0(&d) & (Exp_mask & ~Exp_msk1)
 #endif
 				) {
 			/* The special case */
@@ -615,9 +628,9 @@ dtoa
 			j1 = delta->sign ? 1 : cmp(b, delta);
 			Bfree(delta);
 #ifndef ROUND_BIASED
-			if (j1 == 0 && mode != 1 && !(word1(d) & 1)
+			if (j1 == 0 && mode != 1 && !(word1(&d) & 1)
 #ifdef Honor_FLT_ROUNDS
-				&& rounding >= 1
+				&& Rounding >= 1
 #endif
 								   ) {
 				if (dig == '9')
@@ -632,11 +645,11 @@ dtoa
 				goto ret;
 				}
 #endif
-			if (j < 0 || j == 0 && mode != 1
+			if (j < 0 || (j == 0 && mode != 1
 #ifndef ROUND_BIASED
-							&& !(word1(d) & 1)
+							&& !(word1(&d) & 1)
 #endif
-					) {
+					)) {
 				if (!b->x[0] && b->wds <= 1) {
 #ifdef SET_INEXACT
 					inexact = 0;
@@ -645,7 +658,7 @@ dtoa
 					}
 #ifdef Honor_FLT_ROUNDS
 				if (mode > 1)
-				 switch(rounding) {
+				 switch(Rounding) {
 				  case 0: goto accept_dig;
 				  case 2: goto keep_dig;
 				  }
@@ -653,7 +666,11 @@ dtoa
 				if (j1 > 0) {
 					b = lshift(b, 1);
 					j1 = cmp(b, S);
-					if ((j1 > 0 || j1 == 0 && dig & 1)
+#ifdef ROUND_BIASED
+					if (j1 >= 0 /*)*/
+#else
+					if ((j1 > 0 || (j1 == 0 && dig & 1))
+#endif
 					&& dig++ == '9')
 						goto round_9_up;
 					}
@@ -663,7 +680,7 @@ dtoa
 				}
 			if (j1 > 0) {
 #ifdef Honor_FLT_ROUNDS
-				if (!rounding)
+				if (!Rounding)
 					goto accept_dig;
 #endif
 				if (dig == '9') { /* possible if i == 1 */
@@ -706,14 +723,19 @@ dtoa
 	/* Round off last digit */
 
 #ifdef Honor_FLT_ROUNDS
-	switch(rounding) {
+	switch(Rounding) {
 	  case 0: goto trimzeros;
 	  case 2: goto roundoff;
 	  }
 #endif
 	b = lshift(b, 1);
 	j = cmp(b, S);
-	if (j > 0 || j == 0 && dig & 1) {
+#ifdef ROUND_BIASED
+	if (j >= 0)
+#else
+	if (j > 0 || (j == 0 && dig & 1))
+#endif
+		{
  roundoff:
 		while(*--s == '9')
 			if (s == s0) {
@@ -724,7 +746,9 @@ dtoa
 		++*s++;
 		}
 	else {
+#ifdef Honor_FLT_ROUNDS
  trimzeros:
+#endif
 		while(*--s == '0');
 		s++;
 		}
@@ -739,9 +763,9 @@ dtoa
 #ifdef SET_INEXACT
 	if (inexact) {
 		if (!oldinexact) {
-			word0(d) = Exp_1 + (70 << Exp_shift);
+			word0(&d) = Exp_1 + (70 << Exp_shift);
-			word1(d) = 0;
+			word1(&d) = 0;
-			dval(d) += 1.;
+			dval(&d) += 1.;
 			}
 		}
 	else if (!oldinexact)

gdtoa/g_Qfmt.c

@@ -51,19 +51,24 @@ THIS SOFTWARE.
 
  char*
 #ifdef KR_headers
-g_Qfmt(buf, V, ndig, bufsize) char *buf; char *V; int ndig; unsigned bufsize;
+g_Qfmt(buf, V, ndig, bufsize) char *buf; char *V; int ndig; size_t bufsize;
 #else
-g_Qfmt(char *buf, void *V, int ndig, unsigned bufsize)
+g_Qfmt(char *buf, void *V, int ndig, size_t bufsize)
 #endif
 {
-	static CONST FPI fpi = { 113, 1-16383-113+1, 32766 - 16383 - 113 + 1, 1, 0 };
+	static FPI fpi0 = { 113, 1-16383-113+1, 32766 - 16383 - 113 + 1, 1, 0, Int_max };
 	char *b, *s, *se;
 	ULong bits[4], *L, sign;
 	int decpt, ex, i, mode;
+#ifdef Honor_FLT_ROUNDS
+#include "gdtoa_fltrnds.h"
+#else
+#define fpi &fpi0
+#endif
 
 	if (ndig < 0)
 		ndig = 0;
-	if (bufsize < (unsigned)(ndig + 10))
+	if (bufsize < (size_t)(ndig + 10))
 		return 0;
 
 	L = (ULong*)V;
@@ -109,6 +114,6 @@ g_Qfmt(char *buf, void *V, int ndig, unsigned bufsize)
 			return 0;
 		mode = 0;
 		}
-	s = gdtoa(&fpi, ex, bits, &i, mode, ndig, &decpt, &se);
+	s = gdtoa(fpi, ex, bits, &i, mode, ndig, &decpt, &se);
-	return g__fmt(buf, s, se, decpt, sign);
+	return g__fmt(buf, s, se, decpt, sign, bufsize);
 	}

gdtoa/g__fmt.c

@@ -35,26 +35,77 @@ THIS SOFTWARE.
 #include "locale.h"
 #endif
 
+#ifndef ldus_QNAN0
+#define ldus_QNAN0 0x7fff
+#endif
+#ifndef ldus_QNAN1
+#define ldus_QNAN1 0xc000
+#endif
+#ifndef ldus_QNAN2
+#define ldus_QNAN2 0
+#endif
+#ifndef ldus_QNAN3
+#define ldus_QNAN3 0
+#endif
+#ifndef ldus_QNAN4
+#define ldus_QNAN4 0
+#endif
+
+ const char *InfName[6] = { "Infinity", "infinity", "INFINITY", "Inf", "inf", "INF" };
+ const char *NanName[3] = { "NaN", "nan", "NAN" };
+ ULong NanDflt_Q_D2A[4] = { 0xffffffff, 0xffffffff, 0xffffffff, 0x7fffffff };
+ ULong NanDflt_d_D2A[2] = { d_QNAN1, d_QNAN0 };
+ ULong NanDflt_f_D2A[1] = { f_QNAN };
+ ULong NanDflt_xL_D2A[3] = { 1, 0x80000000, 0x7fff0000 };
+ UShort NanDflt_ldus_D2A[5] = { ldus_QNAN4, ldus_QNAN3, ldus_QNAN2, ldus_QNAN1, ldus_QNAN0 };
+
  char *
 #ifdef KR_headers
-g__fmt(b, s, se, decpt, sign) char *b; char *s; char *se; int decpt; ULong sign;
+g__fmt(b, s, se, decpt, sign, blen) char *b; char *s; char *se; int decpt; ULong sign; size_t blen;
 #else
-g__fmt(char *b, char *s, char *se, int decpt, ULong sign)
+g__fmt(char *b, char *s, char *se, int decpt, ULong sign, size_t blen)
 #endif
 {
 	int i, j, k;
-	char *s0 = s;
+	char *be, *s0;
+	size_t len;
 #ifdef USE_LOCALE
-	char decimalpoint = *localeconv()->decimal_point;
+#ifdef NO_LOCALE_CACHE
+	char *decimalpoint = localeconv()->decimal_point;
+	size_t dlen = strlen(decimalpoint);
 #else
-#define decimalpoint '.'
+	char *decimalpoint;
+	static char *decimalpoint_cache;
+	static size_t dlen;
+	if (!(s0 = decimalpoint_cache)) {
+		s0 = localeconv()->decimal_point;
+		dlen = strlen(s0);
+		if ((decimalpoint_cache = (char*)MALLOC(strlen(s0) + 1))) {
+			strcpy(decimalpoint_cache, s0);
+			s0 = decimalpoint_cache;
+			}
+		}
+	decimalpoint = s0;
 #endif
+#else
+#define dlen 0
+#endif
+	s0 = s;
+	len = (se-s) + dlen + 6; /* 6 = sign + e+dd + trailing null */
+	if (blen < len)
+		goto ret0;
+	be = b + blen - 1;
 	if (sign)
 		*b++ = '-';
 	if (decpt <= -4 || decpt > se - s + 5) {
 		*b++ = *s++;
 		if (*s) {
-			*b++ = decimalpoint;
+#ifdef USE_LOCALE
+			while((*b = *decimalpoint++))
+				++b;
+#else
+			*b++ = '.';
+#endif
 			while((*b = *s++) !=0)
 				b++;
 			}
@@ -69,6 +120,8 @@ g__fmt(char *b, char *s, char *se, int decpt, ULong sign)
 		for(j = 2, k = 10; 10*k <= decpt; j++, k *= 10){}
 		for(;;) {
 			i = decpt / k;
+			if (b >= be)
+				goto ret0;
 			*b++ = i + '0';
 			if (--j <= 0)
 				break;
@@ -78,7 +131,14 @@ g__fmt(char *b, char *s, char *se, int decpt, ULong sign)
 		*b = 0;
 		}
 	else if (decpt <= 0) {
-		*b++ = decimalpoint;
+#ifdef USE_LOCALE
+		while((*b = *decimalpoint++))
+			++b;
+#else
+		*b++ = '.';
+#endif
+		if (be < b - decpt + (se - s))
+			goto ret0;
 		for(; decpt < 0; decpt++)
 			*b++ = '0';
 		while((*b = *s++) != 0)
@@ -87,13 +147,57 @@ g__fmt(char *b, char *s, char *se, int decpt, ULong sign)
 	else {
 		while((*b = *s++) != 0) {
 			b++;
-			if (--decpt == 0 && *s)
+			if (--decpt == 0 && *s) {
-				*b++ = decimalpoint;
+#ifdef USE_LOCALE
+				while(*b = *decimalpoint++)
+					++b;
+#else
+				*b++ = '.';
+#endif
+				}
+			}
+		if (b + decpt > be) {
+ ret0:
+			b = 0;
+			goto ret;
 			}
 		for(; decpt > 0; decpt--)
 			*b++ = '0';
 		*b = 0;
 		}
+ ret:
 	freedtoa(s0);
 	return b;
  	}
+
+ char *
+add_nanbits_D2A(char *b, size_t blen, ULong *bits, int nb)
+{
+	ULong t;
+	char *rv;
+	int i, j;
+	size_t L;
+	static char Hexdig[16] = "0123456789abcdef";
+
+	while(!bits[--nb])
+		if (!nb)
+			return b;
+	L = 8*nb + 3;
+	t = bits[nb];
+	do ++L; while((t >>= 4));
+	if (L > blen)
+		return b;
+	b += L;
+	*--b = 0;
+	rv = b;
+	*--b = /*(*/ ')';
+	for(i = 0; i < nb; ++i) {
+		t = bits[i];
+		for(j = 0; j < 8; ++j, t >>= 4)
+			*--b = Hexdig[t & 0xf];
+		}
+	t = bits[nb];
+	do *--b = Hexdig[t & 0xf]; while(t >>= 4);
+	*--b = '('; /*)*/
+	return rv;
+	}

gdtoa/g_ddfmt.c

@@ -33,9 +33,9 @@ THIS SOFTWARE.
 
  char *
 #ifdef KR_headers
-g_ddfmt(buf, dd, ndig, bufsize) char *buf; double *dd; int ndig; unsigned bufsize;
+g_ddfmt(buf, dd0, ndig, bufsize) char *buf; double *dd0; int ndig; size_t bufsize;
 #else
-g_ddfmt(char *buf, double *dd, int ndig, unsigned bufsize)
+g_ddfmt(char *buf, double *dd0, int ndig, size_t bufsize)
 #endif
 {
 	FPI fpi;
@@ -43,12 +43,28 @@ g_ddfmt(char *buf, double *dd, int ndig, unsigned bufsize)
 	ULong *L, bits0[4], *bits, *zx;
 	int bx, by, decpt, ex, ey, i, j, mode;
 	Bigint *x, *y, *z;
-	double ddx[2];
+	U *dd, ddx[2];
+#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 /*}}*/
+#else /*}{*/
+#define Rounding FPI_Round_near
+#endif /*}}*/
 
-	if (bufsize < 10 || bufsize < (unsigned)(ndig + 8))
+	if (bufsize < 10 || bufsize < (size_t)(ndig + 8))
 		return 0;
 
-	L = (ULong*)dd;
+	dd = (U*)dd0;
+	L = dd->L;
 	if ((L[_0] & 0x7ff00000L) == 0x7ff00000L) {
 		/* Infinity or NaN */
 		if (L[_0] & 0xfffff || L[_1]) {
@@ -73,7 +89,7 @@ g_ddfmt(char *buf, double *dd, int ndig, unsigned bufsize)
 			goto nanret;
 		goto infret;
 		}
-	if (dd[0] + dd[1] == 0.) {
+	if (dval(&dd[0]) + dval(&dd[1]) == 0.) {
 		b = buf;
 #ifndef IGNORE_ZERO_SIGN
 		if (L[_0] & L[2+_0] & 0x80000000L)
@@ -84,16 +100,16 @@ g_ddfmt(char *buf, double *dd, int ndig, unsigned bufsize)
 		return b;
 		}
 	if ((L[_0] & 0x7ff00000L) < (L[2+_0] & 0x7ff00000L)) {
-		ddx[1] = dd[0];
+		dval(&ddx[1]) = dval(&dd[0]);
-		ddx[0] = dd[1];
+		dval(&ddx[0]) = dval(&dd[1]);
 		dd = ddx;
-		L = (ULong*)dd;
+		L = dd->L;
 		}
-	z = d2b(dd[0], &ex, &bx);
+	z = d2b(dval(&dd[0]), &ex, &bx);
-	if (dd[1] == 0.)
+	if (dval(&dd[1]) == 0.)
 		goto no_y;
 	x = z;
-	y = d2b(dd[1], &ey, &by);
+	y = d2b(dval(&dd[1]), &ey, &by);
 	if ( (i = ex - ey) !=0) {
 		if (i > 0) {
 			x = lshift(x, i);
@@ -136,7 +152,7 @@ g_ddfmt(char *buf, double *dd, int ndig, unsigned bufsize)
 		}
 	mode = 2;
 	if (ndig <= 0) {
-		if (bufsize < (unsigned)((int)(fpi.nbits * .301029995664) + 10)) {
+		if (bufsize < (size_t)((int)(fpi.nbits * .301029995664) + 10)) {
 			Bfree(z);
 			return 0;
 			}
@@ -144,11 +160,12 @@ g_ddfmt(char *buf, double *dd, int ndig, unsigned bufsize)
 		}
 	fpi.emin = 1-1023-53+1;
 	fpi.emax = 2046-1023-106+1;
-	fpi.rounding = FPI_Round_near;
+	fpi.rounding = Rounding;
 	fpi.sudden_underflow = 0;
+	fpi.int_max = Int_max;
 	i = STRTOG_Normal;
 	s = gdtoa(&fpi, ex, bits, &i, mode, ndig, &decpt, &se);
-	b = g__fmt(buf, s, se, decpt, z->sign);
+	b = g__fmt(buf, s, se, decpt, z->sign, bufsize);
 	Bfree(z);
 	return b;
 	}

gdtoa/g_dfmt.c

@@ -33,25 +33,32 @@ THIS SOFTWARE.
 
  char*
 #ifdef KR_headers
-g_dfmt(buf, d, ndig, bufsize) char *buf; double *d; int ndig; unsigned bufsize;
+g_dfmt(buf, d, ndig, bufsize) char *buf; double *d; int ndig; size_t bufsize;
 #else
-g_dfmt(char *buf, double *d, int ndig, unsigned bufsize)
+g_dfmt(char *buf, double *d, int ndig, size_t bufsize)
 #endif
 {
-	static CONST FPI fpi = { 53, 1-1023-53+1, 2046-1023-53+1, 1, 0 };
+	static FPI fpi0 = { 53, 1-1023-53+1, 2046-1023-53+1, 1, 0, Int_max };
 	char *b, *s, *se;
 	ULong bits[2], *L, sign;
 	int decpt, ex, i, mode;
+#ifdef Honor_FLT_ROUNDS
+#include "gdtoa_fltrnds.h"
+#else
+#define fpi &fpi0
+#endif
 
 	if (ndig < 0)
 		ndig = 0;
-	if (bufsize < (unsigned)(ndig + 10))
+	if (bufsize < (size_t)(ndig + 10))
 		return 0;
 
 	L = (ULong*)d;
 	sign = L[_0] & 0x80000000L;
 	if ((L[_0] & 0x7ff00000) == 0x7ff00000) {
 		/* Infinity or NaN */
+		if (bufsize < 10)
+			return 0;
 		if (L[_0] & 0xfffff || L[_1]) {
 			return strcp(buf, "NaN");
 			}
@@ -78,12 +85,11 @@ g_dfmt(char *buf, double *d, int ndig, unsigned bufsize)
 		ex = 1;
 	ex -= 0x3ff + 52;
 	mode = 2;
-	if (ndig <= 0) {
+	if (ndig <= 0)
-		if (bufsize < 25)
-			return 0;
 		mode = 0;
-		}
 	i = STRTOG_Normal;
-	s = gdtoa(&fpi, ex, bits, &i, mode, ndig, &decpt, &se);
+	if (sign)
-	return g__fmt(buf, s, se, decpt, sign);
+		i = STRTOG_Normal | STRTOG_Neg;
+	s = gdtoa(fpi, ex, bits, &i, mode, ndig, &decpt, &se);
+	return g__fmt(buf, s, se, decpt, sign, bufsize);
 	}

gdtoa/g_ffmt.c

@@ -33,19 +33,24 @@ THIS SOFTWARE.
 
  char*
 #ifdef KR_headers
-g_ffmt(buf, f, ndig, bufsize) char *buf; float *f; int ndig; unsigned bufsize;
+g_ffmt(buf, f, ndig, bufsize) char *buf; float *f; int ndig; size_t bufsize;
 #else
-g_ffmt(char *buf, float *f, int ndig, unsigned bufsize)
+g_ffmt(char *buf, float *f, int ndig, size_t bufsize)
 #endif
 {
-	static CONST FPI fpi = { 24, 1-127-24+1,  254-127-24+1, 1, 0 };
+	static FPI fpi0 = { 24, 1-127-24+1,  254-127-24+1, 1, 0, 6 };
 	char *b, *s, *se;
 	ULong bits[1], *L, sign;
 	int decpt, ex, i, mode;
+#ifdef Honor_FLT_ROUNDS
+#include "gdtoa_fltrnds.h"
+#else
+#define fpi &fpi0
+#endif
 
 	if (ndig < 0)
 		ndig = 0;
-	if (bufsize < (unsigned)(ndig + 10))
+	if (bufsize < (size_t)(ndig + 10))
 		return 0;
 
 	L = (ULong*)f;
@@ -83,6 +88,6 @@ g_ffmt(char *buf, float *f, int ndig, unsigned bufsize)
 		mode = 0;
 		}
 	i = STRTOG_Normal;
-	s = gdtoa(&fpi, ex, bits, &i, mode, ndig, &decpt, &se);
+	s = gdtoa(fpi, ex, bits, &i, mode, ndig, &decpt, &se);
-	return g__fmt(buf, s, se, decpt, sign);
+	return g__fmt(buf, s, se, decpt, sign, bufsize);
 	}

gdtoa/g_xLfmt.c

@@ -49,19 +49,24 @@ THIS SOFTWARE.
 
  char*
 #ifdef KR_headers
-g_xLfmt(buf, V, ndig, bufsize) char *buf; char *V; int ndig; unsigned bufsize;
+g_xLfmt(buf, V, ndig, bufsize) char *buf; char *V; int ndig; size_t bufsize;
 #else
-g_xLfmt(char *buf, void *V, int ndig, unsigned bufsize)
+g_xLfmt(char *buf, void *V, int ndig, size_t bufsize)
 #endif
 {
-	static CONST FPI fpi = { 64, 1-16383-64+1, 32766 - 16383 - 64 + 1, 1, 0 };
+	static FPI fpi0 = { 64, 1-16383-64+1, 32766 - 16383 - 64 + 1, 1, 0, Int_max };
 	char *b, *s, *se;
 	ULong bits[2], *L, sign;
 	int decpt, ex, i, mode;
+#ifdef Honor_FLT_ROUNDS
+#include "gdtoa_fltrnds.h"
+#else
+#define fpi &fpi0
+#endif
 
 	if (ndig < 0)
 		ndig = 0;
-	if (bufsize < (unsigned)(ndig + 10))
+	if (bufsize < (size_t)(ndig + 10))
 		return 0;
 
 	L = (ULong*)V;
@@ -103,6 +108,6 @@ g_xLfmt(char *buf, void *V, int ndig, unsigned bufsize)
 			return 0;
 		mode = 0;
 		}
-	s = gdtoa(&fpi, ex, bits, &i, mode, ndig, &decpt, &se);
+	s = gdtoa(fpi, ex, bits, &i, mode, ndig, &decpt, &se);
-	return g__fmt(buf, s, se, decpt, sign);
+	return g__fmt(buf, s, se, decpt, sign, bufsize);
 	}

gdtoa/g_xfmt.c

@@ -53,20 +53,25 @@ THIS SOFTWARE.
 
  char*
 #ifdef KR_headers
-g_xfmt(buf, V, ndig, bufsize) char *buf; char *V; int ndig; unsigned bufsize;
+g_xfmt(buf, V, ndig, bufsize) char *buf; char *V; int ndig; size_t bufsize;
 #else
-g_xfmt(char *buf, void *V, int ndig, unsigned bufsize)
+g_xfmt(char *buf, void *V, int ndig, size_t bufsize)
 #endif
 {
-	static CONST FPI fpi = { 64, 1-16383-64+1, 32766 - 16383 - 64 + 1, 1, 0 };
+	static FPI fpi0 = { 64, 1-16383-64+1, 32766 - 16383 - 64 + 1, 1, 0, Int_max };
 	char *b, *s, *se;
 	ULong bits[2], sign;
 	UShort *L;
 	int decpt, ex, i, mode;
+#ifdef Honor_FLT_ROUNDS
+#include "gdtoa_fltrnds.h"
+#else
+#define fpi &fpi0
+#endif
 
 	if (ndig < 0)
 		ndig = 0;
-	if (bufsize < (unsigned)(ndig + 10))
+	if (bufsize < (size_t)(ndig + 10))
 		return 0;
 
 	L = (UShort *)V;
@@ -76,14 +81,14 @@ g_xfmt(char *buf, void *V, int ndig, unsigned bufsize)
 	if ( (ex = L[_0] & 0x7fff) !=0) {
 		if (ex == 0x7fff) {
 			/* Infinity or NaN */
-			if (bits[0] | bits[1])
+			if (!bits[0] && bits[1]== 0x80000000) {
-				b = strcp(buf, "NaN");
-			else {
 				b = buf;
 				if (sign)
 					*b++ = '-';
 				b = strcp(b, "Infinity");
 				}
+			else
+				b = strcp(buf, "NaN");
 			return b;
 			}
 		i = STRTOG_Normal;
@@ -109,6 +114,6 @@ g_xfmt(char *buf, void *V, int ndig, unsigned bufsize)
 			return 0;
 		mode = 0;
 		}
-	s = gdtoa(&fpi, ex, bits, &i, mode, ndig, &decpt, &se);
+	s = gdtoa(fpi, ex, bits, &i, mode, ndig, &decpt, &se);
-	return g__fmt(buf, s, se, decpt, sign);
+	return g__fmt(buf, s, se, decpt, sign, bufsize);
 	}

gdtoa/gdtoa.c

@@ -30,7 +30,6 @@ THIS SOFTWARE.
  * with " at " changed at "@" and " dot " changed to ".").	*/
 
 #include "gdtoaimp.h"
-#include <limits.h>
 
  static Bigint *
 #ifdef KR_headers
@@ -62,7 +61,7 @@ bitstob(ULong *bits, int nbits, int *bbits)
 		*x++ = (*bits >> 16) & ALL_ON;
 #endif
 		} while(++bits <= be);
-	i = (int)(x - x0);
+	i = x - x0;
 	while(!x0[--i])
 		if (!i) {
 			b->wds = 0;
@@ -116,14 +115,15 @@ gdtoa
 	FPI *fpi; int be; ULong *bits;
 	int *kindp, mode, ndigits, *decpt; char **rve;
 #else
-	(CONST FPI *fpi, int be, ULong *bits, int *kindp, int mode, int ndigits, int *decpt, char **rve)
+	(FPI *fpi, int be, ULong *bits, int *kindp, int mode, int ndigits, int *decpt, char **rve)
 #endif
 {
  /*	Arguments ndigits and decpt are similar to the second and third
 	arguments of ecvt and fcvt; trailing zeros are suppressed from
 	the returned string.  If not null, *rve is set to point
 	to the end of the return value.  If d is +-Infinity or NaN,
-	then *decpt is set to INT_MAX.
+	then *decpt is set to 9999.
+	be = exponent: value = (integer represented by bits) * (2 to the power of be).
 
 	mode:
 		0 ==> shortest string that yields d when read in
@@ -159,8 +159,8 @@ gdtoa
 	Long L;
 	Bigint *b, *b1, *delta, *mlo, *mhi, *mhi1, *S;
 	double d2, ds;
-	U d, eps;
 	char *s, *s0;
+	U d, eps;
 
 #ifndef MULTIPLE_THREADS
 	if (dtoa_result) {
@@ -177,10 +177,10 @@ gdtoa
 	  case STRTOG_Denormal:
 		break;
 	  case STRTOG_Infinite:
-		*decpt = INT_MAX;
+		*decpt = -32768;
 		return nrv_alloc("Infinity", rve, 8);
 	  case STRTOG_NaN:
-		*decpt = INT_MAX;
+		*decpt = -32768;
 		return nrv_alloc("NaN", rve, 3);
 	  default:
 		return 0;
@@ -199,21 +199,21 @@ gdtoa
 		return nrv_alloc("0", rve, 1);
 		}
 
-	dval(d) = b2d(b, &i);
+	dval(&d) = b2d(b, &i);
 	i = be + bbits - 1;
-	word0(d) &= Frac_mask1;
+	word0(&d) &= Frac_mask1;
-	word0(d) |= Exp_11;
+	word0(&d) |= Exp_11;
 #ifdef IBM
-	if ( (j = 11 - hi0bits(word0(d) & Frac_mask)) !=0)
+	if ( (j = 11 - hi0bits(word0(&d) & Frac_mask)) !=0)
-		dval(d) /= 1 << j;
+		dval(&d) /= 1 << j;
 #endif
 
 	/* log(x)	~=~ log(1.5) + (x-1.5)/1.5
 	 * log10(x)	 =  log(x) / log(10)
 	 *		~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
-	 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
+	 * log10(&d) = (i-Bias)*log(2)/log(10) + log10(d2)
 	 *
-	 * This suggests computing an approximation k to log10(d) by
+	 * This suggests computing an approximation k to log10(&d) by
 	 *
 	 * k = (i - Bias)*0.301029995663981
 	 *	+ ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
@@ -233,7 +233,7 @@ gdtoa
 	i <<= 2;
 	i += j;
 #endif
-	ds = (dval(d)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
+	ds = (dval(&d)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
 
 	/* correct assumption about exponent range */
 	if ((j = i) < 0)
@@ -248,13 +248,13 @@ gdtoa
 #ifdef IBM
 	j = be + bbits - 1;
 	if ( (j1 = j & 3) !=0)
-		dval(d) *= 1 << j1;
+		dval(&d) *= 1 << j1;
-	word0(d) += j << Exp_shift - 2 & Exp_mask;
+	word0(&d) += j << Exp_shift - 2 & Exp_mask;
 #else
-	word0(d) += (be + bbits - 1) << Exp_shift;
+	word0(&d) += (be + bbits - 1) << Exp_shift;
 #endif
 	if (k >= 0 && k <= Ten_pmax) {
-		if (dval(d) < tens[k])
+		if (dval(&d) < tens[k])
 			k--;
 		k_check = 0;
 		}
@@ -284,11 +284,14 @@ gdtoa
 		mode -= 4;
 		try_quick = 0;
 		}
+	else if (i >= -4 - Emin || i < Emin)
+		try_quick = 0;
 	leftright = 1;
+	ilim = ilim1 = -1;	/* Values for cases 0 and 1; done here to */
+				/* silence erroneous "gcc -Wall" warning. */
 	switch(mode) {
 		case 0:
 		case 1:
-			ilim = ilim1 = -1;
 			i = (int)(nbits * .30103) + 3;
 			ndigits = 0;
 			break;
@@ -330,10 +333,10 @@ gdtoa
 		/* Try to get by with floating-point arithmetic. */
 
 		i = 0;
-		d2 = dval(d);
+		d2 = dval(&d);
 #ifdef IBM
-		if ( (j = 11 - hi0bits(word0(d) & Frac_mask)) !=0)
+		if ( (j = 11 - hi0bits(word0(&d) & Frac_mask)) !=0)
-			dval(d) /= 1 << j;
+			dval(&d) /= 1 << j;
 #endif
 		k0 = k;
 		ilim0 = ilim;
@@ -344,7 +347,7 @@ gdtoa
 			if (j & Bletch) {
 				/* prevent overflows */
 				j &= Bletch - 1;
-				dval(d) /= bigtens[n_bigtens-1];
+				dval(&d) /= bigtens[n_bigtens-1];
 				ieps++;
 				}
 			for(; j; j >>= 1, i++)
@@ -356,30 +359,30 @@ gdtoa
 		else  {
 			ds = 1.;
 			if ( (j1 = -k) !=0) {
-				dval(d) *= tens[j1 & 0xf];
+				dval(&d) *= tens[j1 & 0xf];
 				for(j = j1 >> 4; j; j >>= 1, i++)
 					if (j & 1) {
 						ieps++;
-						dval(d) *= bigtens[i];
+						dval(&d) *= bigtens[i];
 						}
 				}
 			}
-		if (k_check && dval(d) < 1. && ilim > 0) {
+		if (k_check && dval(&d) < 1. && ilim > 0) {
 			if (ilim1 <= 0)
 				goto fast_failed;
 			ilim = ilim1;
 			k--;
-			dval(d) *= 10.;
+			dval(&d) *= 10.;
 			ieps++;
 			}
-		dval(eps) = ieps*dval(d) + 7.;
+		dval(&eps) = ieps*dval(&d) + 7.;
-		word0(eps) -= (P-1)*Exp_msk1;
+		word0(&eps) -= (P-1)*Exp_msk1;
 		if (ilim == 0) {
 			S = mhi = 0;
-			dval(d) -= 5.;
+			dval(&d) -= 5.;
-			if (dval(d) > dval(eps))
+			if (dval(&d) > dval(&eps))
 				goto one_digit;
-			if (dval(d) < -dval(eps))
+			if (dval(&d) < -dval(&eps))
 				goto no_digits;
 			goto fast_failed;
 			}
@@ -388,42 +391,40 @@ gdtoa
 			/* Use Steele & White method of only
 			 * generating digits needed.
 			 */
-			dval(eps) = ds*0.5/tens[ilim-1] - dval(eps);
+			dval(&eps) = ds*0.5/tens[ilim-1] - dval(&eps);
 			for(i = 0;;) {
-				L = (Long)(dval(d)/ds);
+				L = (Long)(dval(&d)/ds);
-				dval(d) -= L*ds;
+				dval(&d) -= L*ds;
 				*s++ = '0' + (int)L;
-				if (dval(d) < dval(eps)) {
+				if (dval(&d) < dval(&eps)) {
-					if (dval(d))
+					if (dval(&d))
 						inex = STRTOG_Inexlo;
 					goto ret1;
 					}
-				if (ds - dval(d) < dval(eps))
+				if (ds - dval(&d) < dval(&eps))
 					goto bump_up;
 				if (++i >= ilim)
 					break;
-				dval(eps) *= 10.;
+				dval(&eps) *= 10.;
-				dval(d) *= 10.;
+				dval(&d) *= 10.;
 				}
 			}
 		else {
 #endif
 			/* Generate ilim digits, then fix them up. */
-			dval(eps) *= tens[ilim-1];
+			dval(&eps) *= tens[ilim-1];
-			for(i = 1;; i++, dval(d) *= 10.) {
+			for(i = 1;; i++, dval(&d) *= 10.) {
-				if ( (L = (Long)(dval(d)/ds)) !=0)
+				if ( (L = (Long)(dval(&d)/ds)) !=0)
-					dval(d) -= L*ds;
+					dval(&d) -= L*ds;
 				*s++ = '0' + (int)L;
 				if (i == ilim) {
 					ds *= 0.5;
-					if (dval(d) > ds + dval(eps))
+					if (dval(&d) > ds + dval(&eps))
 						goto bump_up;
-					else if (dval(d) < ds - dval(eps)) {
+					else if (dval(&d) < ds - dval(&eps)) {
-						while(*--s == '0'){}
+						if (dval(&d))
-						s++;
-						if (dval(d))
 							inex = STRTOG_Inexlo;
-						goto ret1;
+						goto clear_trailing0;
 						}
 					break;
 					}
@@ -433,34 +434,34 @@ gdtoa
 #endif
  fast_failed:
 		s = s0;
-		dval(d) = d2;
+		dval(&d) = d2;
 		k = k0;
 		ilim = ilim0;
 		}
 
 	/* Do we have a "small" integer? */
 
-	if (be >= 0 && k <= Int_max) {
+	if (be >= 0 && k <= fpi->int_max) {
 		/* Yes. */
 		ds = tens[k];
 		if (ndigits < 0 && ilim <= 0) {
 			S = mhi = 0;
-			if (ilim < 0 || dval(d) <= 5*ds)
+			if (ilim < 0 || dval(&d) <= 5*ds)
 				goto no_digits;
 			goto one_digit;
 			}
-		for(i = 1;; i++, dval(d) *= 10.) {
+		for(i = 1;; i++, dval(&d) *= 10.) {
-			L = (Long)(dval(d) / ds);
+			L = (Long)(dval(&d) / ds);
-			dval(d) -= L*ds;
+			dval(&d) -= L*ds;
 #ifdef Check_FLT_ROUNDS
 			/* If FLT_ROUNDS == 2, L will usually be high by 1 */
-			if (dval(d) < 0) {
+			if (dval(&d) < 0) {
 				L--;
-				dval(d) += ds;
+				dval(&d) += ds;
 				}
 #endif
 			*s++ = '0' + (int)L;
-			if (dval(d) == 0.)
+			if (dval(&d) == 0.)
 				break;
 			if (i == ilim) {
 				if (rdir) {
@@ -469,8 +470,13 @@ gdtoa
 					inex = STRTOG_Inexlo;
 					goto ret1;
 					}
-				dval(d) += dval(d);
+				dval(&d) += dval(&d);
-				if (dval(d) > ds || dval(d) == ds && L & 1) {
+#ifdef ROUND_BIASED
+				if (dval(&d) >= ds)
+#else
+				if (dval(&d) > ds || (dval(&d) == ds && L & 1))
+#endif
+					{
  bump_up:
 					inex = STRTOG_Inexhi;
 					while(*--s == '9')
@@ -481,8 +487,12 @@ gdtoa
 							}
 					++*s++;
 					}
-				else
+				else {
 					inex = STRTOG_Inexlo;
+ clear_trailing0:
+					while(*--s == '0'){}
+					++s;
+					}
 				break;
 				}
 			}
@@ -493,13 +503,15 @@ gdtoa
 	m5 = b5;
 	mhi = mlo = 0;
 	if (leftright) {
-		if (mode < 2) {
 		i = nbits - bbits;
-			if (be - i++ < fpi->emin)
+		if (be - i++ < fpi->emin && mode != 3 && mode != 5) {
 			/* denormal */
 			i = be - fpi->emin + 1;
+			if (mode >= 2 && ilim > 0 && ilim < i)
+				goto small_ilim;
 			}
-		else {
+		else if (mode >= 2) {
+ small_ilim:
 			j = ilim - 1;
 			if (m5 >= j)
 				m5 -= j;
@@ -560,28 +572,11 @@ gdtoa
 	 * and for all and pass them and a shift to quorem, so it
 	 * can do shifts and ors to compute the numerator for q.
 	 */
-#ifdef Pack_32
+	i = ((s5 ? hi0bits(S->x[S->wds-1]) : ULbits - 1) - s2 - 4) & kmask;
-	if ( (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f) !=0)
-		i = 32 - i;
-#else
-	if ( (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf) !=0)
-		i = 16 - i;
-#endif
-	if (i > 4) {
-		i -= 4;
-		b2 += i;
 	m2 += i;
-		s2 += i;
+	if ((b2 += i) > 0)
-		}
-	else if (i < 4) {
-		i += 28;
-		b2 += i;
-		m2 += i;
-		s2 += i;
-		}
-	if (b2 > 0)
 		b = lshift(b, b2);
-	if (s2 > 0)
+	if ((s2 += i) > 0)
 		S = lshift(S, s2);
 	if (k_check) {
 		if (cmp(b,S) < 0) {
@@ -646,11 +641,11 @@ gdtoa
 				goto ret;
 				}
 #endif
-			if (j < 0 || j == 0 && !mode
+			if (j < 0 || (j == 0 && !mode
 #ifndef ROUND_BIASED
 							&& !(bits[0] & 1)
 #endif
-					) {
+					)) {
 				if (rdir && (b->wds > 1 || b->x[0])) {
 					if (rdir == 2) {
 						inex = STRTOG_Inexlo;
@@ -673,7 +668,11 @@ gdtoa
 				if (j1 > 0) {
 					b = lshift(b, 1);
 					j1 = cmp(b, S);
-					if ((j1 > 0 || j1 == 0 && dig & 1)
+#ifdef ROUND_BIASED
+					if (j1 >= 0 /*)*/
+#else
+					if ((j1 > 0 || (j1 == 0 && dig & 1))
+#endif
 					&& dig++ == '9')
 						goto round_9_up;
 					inex = STRTOG_Inexhi;
@@ -718,13 +717,18 @@ gdtoa
 	/* Round off last digit */
 
 	if (rdir) {
-		if (rdir == 2 || b->wds <= 1 && !b->x[0])
+		if (rdir == 2 || (b->wds <= 1 && !b->x[0]))
 			goto chopzeros;
 		goto roundoff;
 		}
 	b = lshift(b, 1);
 	j = cmp(b, S);
-	if (j > 0 || j == 0 && dig & 1) {
+#ifdef ROUND_BIASED
+	if (j >= 0)
+#else
+	if (j > 0 || (j == 0 && dig & 1))
+#endif
+		{
  roundoff:
 		inex = STRTOG_Inexhi;
 		while(*--s == '9')
@@ -740,7 +744,7 @@ gdtoa
 		if (b->wds > 1 || b->x[0])
 			inex = STRTOG_Inexlo;
 		while(*--s == '0'){}
-		s++;
+		++s;
 		}
  ret:
 	Bfree(S);

gdtoa/gdtoa.h

@@ -66,16 +66,13 @@ THIS SOFTWARE.
 #else
 #include "arith.h"
 #endif
+#include <stddef.h> /* for size_t */
 
-/* [RH] On 64-bit Unix systems, long is a 64-bit type. I do not that is
- * is what is desired here, so I sacrifice compatibility with systems
- * that use 16-bit ints (oh no!) and make Long an int instead.
- */
 #ifndef Long
-typedef int Long;
+#define Long int
 #endif
 #ifndef ULong
-typedef unsigned int ULong;
+typedef unsigned Long ULong;
 #endif
 #ifndef UShort
 typedef unsigned short UShort;
@@ -111,9 +108,9 @@ typedef unsigned short UShort;
 
 	/* 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_Underflow= 0x40,
 	STRTOG_Overflow	= 0x80
@@ -126,6 +123,7 @@ FPI {
 	int emax;
 	int rounding;
 	int sudden_underflow;
+	int int_max;
 	} FPI;
 
 enum {	/* FPI.rounding values: same as FLT_ROUNDS */
@@ -135,29 +133,31 @@ enum {	/* FPI.rounding values: same as FLT_ROUNDS */
 	FPI_Round_down = 3
 	};
 
-// Our strtod is not the CRT's strtod.
-#include <stdlib.h>
-#define strtod mystrtod
-
 #ifdef __cplusplus
 extern "C" {
 #endif
 
 extern char* dtoa  ANSI((double d, int mode, int ndigits, int *decpt,
 			int *sign, char **rve));
-extern char* gdtoa ANSI((CONST FPI *fpi, int be, ULong *bits, int *kindp,
+extern char* gdtoa ANSI((FPI *fpi, int be, ULong *bits, int *kindp,
 			int mode, int ndigits, int *decpt, char **rve));
 extern void freedtoa ANSI((char*));
 extern float  strtof ANSI((CONST char *, char **));
 extern double strtod ANSI((CONST char *, char **));
-extern int strtodg ANSI((CONST char*, char**, CONST FPI*, Long*, ULong*));
+extern int strtodg ANSI((CONST char*, char**, FPI*, Long*, ULong*));
 
-extern char*	g_ddfmt  ANSI((char*, double*, int, unsigned));
+extern char*	g_ddfmt   ANSI((char*, double*, int, size_t));
-extern char*	g_dfmt   ANSI((char*, double*, int, unsigned));
+extern char*	g_ddfmt_p ANSI((char*, double*,	int, size_t, int));
-extern char*	g_ffmt   ANSI((char*, float*,  int, unsigned));
+extern char*	g_dfmt    ANSI((char*, double*, int, size_t));
-extern char*	g_Qfmt   ANSI((char*, void*,   int, unsigned));
+extern char*	g_dfmt_p  ANSI((char*, double*,	int, size_t, int));
-extern char*	g_xfmt   ANSI((char*, void*,   int, unsigned));
+extern char*	g_ffmt    ANSI((char*, float*,  int, size_t));
-extern char*	g_xLfmt  ANSI((char*, void*,   int, unsigned));
+extern char*	g_ffmt_p  ANSI((char*, float*,	int, size_t, int));
+extern char*	g_Qfmt    ANSI((char*, void*,   int, size_t));
+extern char*	g_Qfmt_p  ANSI((char*, void*,	int, size_t, int));
+extern char*	g_xfmt    ANSI((char*, void*,   int, size_t));
+extern char*	g_xfmt_p  ANSI((char*, void*,	int, size_t, int));
+extern char*	g_xLfmt   ANSI((char*, void*,   int, size_t));
+extern char*	g_xLfmt_p ANSI((char*, void*,	int, size_t, int));
 
 extern int	strtoId  ANSI((CONST char*, char**, double*, double*));
 extern int	strtoIdd ANSI((CONST char*, char**, double*, double*));

gdtoa/gdtoa_fltrnds.h

@@ -0,0 +1,18 @@
+	FPI *fpi, fpi1;
+	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 /*}}*/
+	fpi = &fpi0;
+	if (Rounding != 1) {
+		fpi1 = fpi0;
+		fpi = &fpi1;
+		fpi1.rounding = Rounding;
+		}

gdtoa/gdtoaimp.h

@@ -89,12 +89,18 @@ THIS SOFTWARE.
  * #define IBM for IBM mainframe-style floating-point arithmetic.
  * #define VAX for VAX-style floating-point arithmetic (D_floating).
  * #define No_leftright to omit left-right logic in fast floating-point
- *	computation of dtoa.
+ *	computation of dtoa and gdtoa.  This will cause modes 4 and 5 to be
+ *	treated the same as modes 2 and 3 for some inputs.
  * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3.
  * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
  *	that use extended-precision instructions to compute rounded
  *	products and quotients) with IBM.
- * #define ROUND_BIASED for IEEE-format with biased rounding.
+ * #define ROUND_BIASED for IEEE-format with biased rounding and arithmetic
+ *	that rounds toward +Infinity.
+ * #define ROUND_BIASED_without_Round_Up for IEEE-format with biased
+ *	rounding when the underlying floating-point arithmetic uses
+ *	unbiased rounding.  This prevent using ordinary floating-point
+ *	arithmetic when the result could be computed with one rounding error.
  * #define Inaccurate_Divide for IEEE-format with correctly rounded
  *	products but inaccurate quotients, e.g., for Intel i860.
  * #define NO_LONG_LONG on machines that do not have a "long long"
@@ -113,7 +119,12 @@ THIS SOFTWARE.
  * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n)
  *	if memory is available and otherwise does something you deem
  *	appropriate.  If MALLOC is undefined, malloc will be invoked
- *	directly -- and assumed always to succeed.
+ *	directly -- and assumed always to succeed.  Similarly, if you
+ *	want something other than the system's free() to be called to
+ *	recycle memory acquired from MALLOC, #define FREE to be the
+ *	name of the alternate routine.  (FREE or free is only called in
+ *	pathological cases, e.g., in a gdtoa call after a gdtoa return in
+ *	mode 3 with thousands of digits requested.)
  * #define Omit_Private_Memory to omit logic (added Jan. 1998) for making
  *	memory allocations from a private pool of memory when possible.
  *	When used, the private pool is PRIVATE_MEM bytes long:  2304 bytes,
@@ -126,8 +137,10 @@ THIS SOFTWARE.
  *	conversions of IEEE doubles in single-threaded executions with
  *	8-byte pointers, PRIVATE_MEM >= 7400 appears to suffice; with
  *	4-byte pointers, PRIVATE_MEM >= 7112 appears adequate.
- * #define INFNAN_CHECK on IEEE systems to cause strtod to check for
+ * #define NO_INFNAN_CHECK if you do not wish to have INFNAN_CHECK
- *	Infinity and NaN (case insensitively).
+ *	#defined automatically on IEEE systems.  On such systems,
+ *	when INFNAN_CHECK is #defined, strtod checks
+ *	for 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 (optionally
@@ -151,7 +164,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.
@@ -220,6 +233,10 @@ THIS SOFTWARE.
 #include "gd_qnan.h"
 #endif
 
+#ifdef Honor_FLT_ROUNDS
+#include <fenv.h>
+#endif
+
 #ifdef DEBUG
 #include "stdio.h"
 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
@@ -229,13 +246,13 @@ THIS SOFTWARE.
 #include "string.h"
 
 #ifdef KR_headers
-#define KR_VOID char
+#define Char char
 #else
-#define KR_VOID void
+#define Char void
 #endif
 
 #ifdef MALLOC
-extern KR_VOID *MALLOC ANSI((size_t));
+extern Char *MALLOC ANSI((size_t));
 #else
 #define MALLOC malloc
 #endif
@@ -307,19 +324,19 @@ extern "C" {
 #endif
 
 #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1
-#error Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
+Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
 #endif
 
 typedef union { double d; ULong L[2]; } U;
 
 #ifdef IEEE_8087
-#define word0(x) x.L[1]
+#define word0(x) (x)->L[1]
-#define word1(x) x.L[0]
+#define word1(x) (x)->L[0]
 #else
-#define word0(x) x.L[0]
+#define word0(x) (x)->L[0]
-#define word1(x) x.L[1]
+#define word1(x) (x)->L[1]
 #endif
-#define dval(x)  x.d
+#define dval(x) (x)->d
 
 /* The following definition of Storeinc is appropriate for MIPS processors.
  * An alternative that might be better on some machines is
@@ -433,6 +450,11 @@ typedef union { double d; ULong L[2]; } U;
 
 #ifndef IEEE_Arith
 #define ROUND_BIASED
+#else
+#ifdef ROUND_BIASED_without_Round_Up
+#undef  ROUND_BIASED
+#define ROUND_BIASED
+#endif
 #endif
 
 #ifdef RND_PRODQUOT
@@ -488,12 +510,12 @@ extern double rnd_prod(double, double), rnd_quot(double, double);
 #define ALL_ON 0xffff
 #endif
 
-#ifndef MULTIPLE_THREADS
+//#ifndef MULTIPLE_THREADS
 #define ACQUIRE_DTOA_LOCK(n)	/*nothing*/
 #define FREE_DTOA_LOCK(n)	/*nothing*/
-#endif
+//#endif
 
-#define Kmax (sizeof(size_t) << 3)
+#define Kmax 9
 
  struct
 Bigint {
@@ -516,12 +538,15 @@ extern void memcpy_D2A ANSI((void*, const void*, size_t));
 
 #define Balloc Balloc_D2A
 #define Bfree Bfree_D2A
+#define InfName InfName_D2A
+#define NanName NanName_D2A
 #define ULtoQ ULtoQ_D2A
 #define ULtof ULtof_D2A
 #define ULtod ULtod_D2A
 #define ULtodd ULtodd_D2A
 #define ULtox ULtox_D2A
 #define ULtoxL ULtoxL_D2A
+#define add_nanbits add_nanbits_D2A
 #define any_on any_on_D2A
 #define b2d b2d_D2A
 #define bigtens bigtens_D2A
@@ -560,9 +585,11 @@ extern void memcpy_D2A ANSI((void*, const void*, size_t));
 #define trailz trailz_D2A
 #define ulp ulp_D2A
 
+ extern char *add_nanbits ANSI((char*, size_t, ULong*, int));
  extern char *dtoa_result;
  extern CONST double bigtens[], tens[], tinytens[];
  extern unsigned char hexdig[];
+ extern const char *InfName[6], *NanName[3];
 
  extern Bigint *Balloc ANSI((int));
  extern void Bfree ANSI((Bigint*));
@@ -577,14 +604,14 @@ 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));
- extern char *g__fmt ANSI((char*, char*, char*, int, ULong));
+ extern char *g__fmt ANSI((char*, char*, char*, int, ULong, size_t));
- extern int gethex ANSI((CONST char**, CONST FPI*, Long*, Bigint**, int));
+ extern int gethex ANSI((CONST char**, FPI*, Long*, Bigint**, int));
  extern void hexdig_init_D2A(Void);
- extern int hexnan ANSI((CONST char**, CONST FPI*, ULong*));
+ extern int hexnan ANSI((CONST char**, FPI*, ULong*));
  extern int hi0bits_D2A ANSI((ULong));
  extern Bigint *i2b ANSI((int));
  extern Bigint *increment ANSI((Bigint*));
@@ -599,14 +626,14 @@ extern void memcpy_D2A ANSI((void*, const void*, size_t));
  extern double ratio ANSI((Bigint*, Bigint*));
  extern void rshift ANSI((Bigint*, int));
  extern char *rv_alloc ANSI((int));
- extern Bigint *s2b ANSI((CONST char*, int, int, ULong));
+ extern Bigint *s2b ANSI((CONST char*, int, int, ULong, int));
  extern Bigint *set_ones ANSI((Bigint*, int));
  extern char *strcp ANSI((char*, const char*));
- extern int strtoIg ANSI((CONST char*, char**, CONST FPI*, Long*, Bigint**, int*));
+ extern int strtoIg ANSI((CONST char*, char**, FPI*, Long*, Bigint**, int*));
  extern double strtod ANSI((const char *s00, char **se));
  extern Bigint *sum ANSI((Bigint*, Bigint*));
  extern int trailz ANSI((Bigint*));
- extern double ulp ANSI((double));
+ extern double ulp ANSI((U*));
 
 #ifdef __cplusplus
 }
@@ -621,6 +648,10 @@ extern void memcpy_D2A ANSI((void*, const void*, size_t));
  * (On HP Series 700/800 machines, -DNAN_WORD0=0x7ff40000 works.)
  */
 #ifdef IEEE_Arith
+#ifndef NO_INFNAN_CHECK
+#undef INFNAN_CHECK
+#define INFNAN_CHECK
+#endif
 #ifdef IEEE_MC68k
 #define _0 0
 #define _1 1

gdtoa/gethex.c

@@ -40,22 +40,34 @@ THIS SOFTWARE.
 gethex(sp, fpi, exp, bp, sign)
 	CONST char **sp; FPI *fpi; Long *exp; Bigint **bp; int sign;
 #else
-gethex( CONST char **sp, CONST FPI *fpi, Long *exp, Bigint **bp, int sign)
+gethex( CONST char **sp, FPI *fpi, Long *exp, Bigint **bp, int sign)
 #endif
 {
 	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
-	unsigned char decimalpoint = *localeconv()->decimal_point;
+	int i;
+#ifdef NO_LOCALE_CACHE
+	const unsigned char *decimalpoint = (unsigned char*)localeconv()->decimal_point;
 #else
-#define decimalpoint '.'
+	const unsigned char *decimalpoint;
+	static unsigned char *decimalpoint_cache;
+	if (!(s0 = decimalpoint_cache)) {
+		s0 = (unsigned char*)localeconv()->decimal_point;
+		if ((decimalpoint_cache = (char*)MALLOC(strlen(s0) + 1))) {
+			strcpy(decimalpoint_cache, s0);
+			s0 = decimalpoint_cache;
+			}
+		}
+	decimalpoint = s0;
+#endif
 #endif
 
-	if (!hexdig['0'])
+	/**** if (!hexdig['0']) hexdig_init_D2A(); ****/
-		hexdig_init_D2A();
+	*bp = 0;
 	havedig = 0;
 	s0 = *(CONST unsigned char **)sp + 2;
 	while(s0[havedig] == '0')
@@ -65,11 +77,21 @@ gethex( CONST char **sp, CONST FPI *fpi, Long *exp, Bigint **bp, int sign)
 	decpt = 0;
 	zret = 0;
 	e = 0;
-	if (!hexdig[*s]) {
+	if (hexdig[*s])
+		havedig++;
+	else {
 		zret = 1;
-		if (*s != decimalpoint)
+#ifdef USE_LOCALE
+		for(i = 0; decimalpoint[i]; ++i) {
+			if (s[i] != decimalpoint[i])
+				goto pcheck;
+			}
+		decpt = s += i;
+#else
+		if (*s != '.')
 			goto pcheck;
 		decpt = ++s;
+#endif
 		if (!hexdig[*s])
 			goto pcheck;
 		while(*s == '0')
@@ -81,19 +103,28 @@ gethex( CONST char **sp, CONST FPI *fpi, Long *exp, Bigint **bp, int sign)
 		}
 	while(hexdig[*s])
 		s++;
-	if (*s == decimalpoint && !decpt) {
+#ifdef USE_LOCALE
+	if (*s == *decimalpoint && !decpt) {
+		for(i = 1; decimalpoint[i]; ++i) {
+			if (s[i] != decimalpoint[i])
+				goto pcheck;
+			}
+		decpt = s += i;
+#else
+	if (*s == '.' && !decpt) {
 		decpt = ++s;
+#endif
 		while(hexdig[*s])
 			s++;
-		}
+		}/*}*/
 	if (decpt)
 		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,29 +137,87 @@ gethex( CONST char **sp, CONST 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;
+		*sp = (char*)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;
 			  }
-	n = (int)(s1 - s0 - 1);
+			goto retz;
-	for(k = 0; n > 7; n >>= 1)
+ 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 > (1 << (kshift-2)) - 1; n >>= 1)
 		k++;
 	b = Balloc(k);
 	x = b->x;
 	n = 0;
 	L = 0;
+#ifdef USE_LOCALE
+	for(i = 0; decimalpoint[i+1]; ++i);
+#endif
 	while(s1 > s0) {
-		if (*--s1 == decimalpoint)
+#ifdef USE_LOCALE
+		if (*--s1 == decimalpoint[i]) {
+			s1 -= i;
 			continue;
-		if (n == 32) {
+			}
+#else
+		if (*--s1 == '.')
+			continue;
+#endif
+		if (n == ULbits) {
 			*x++ = L;
 			L = 0;
 			n = 0;
@@ -137,8 +226,8 @@ gethex( CONST char **sp, CONST FPI *fpi, Long *exp, Bigint **bp, int sign)
 		n += 4;
 		}
 	*x++ = L;
-	b->wds = n = (int)(x - b->x);
+	b->wds = n = x - b->x;
-	n = 32*n - hi0bits(L);
+	n = ULbits*n - hi0bits(L);
 	nbits = fpi->nbits;
 	lostbits = 0;
 	x = b->x;
@@ -149,7 +238,7 @@ gethex( CONST char **sp, CONST 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 +254,10 @@ gethex( CONST char **sp, CONST 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 +277,22 @@ gethex( CONST char **sp, CONST 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;
@@ -214,7 +313,7 @@ gethex( CONST char **sp, CONST FPI *fpi, Long *exp, Bigint **bp, int sign)
 			break;
 		  case FPI_Round_near:
 			if (lostbits & 2
-			 && (lostbits & 1) | x[0] & 1)
+			 && (lostbits | x[0]) & 1)
 				up = 1;
 			break;
 		  case FPI_Round_up:
@@ -233,8 +332,8 @@ gethex( CONST char **sp, CONST FPI *fpi, Long *exp, Bigint **bp, int sign)
 					irv =  STRTOG_Normal;
 				}
 			else if (b->wds > k
-			 || (n = nbits & kmask) !=0
+			 || ((n = nbits & kmask) !=0
-			     && hi0bits(x[k-1]) < 32-n) {
+			      && hi0bits(x[k-1]) < 32-n)) {
 				rshift(b,1);
 				if (++e > fpi->emax)
 					goto ovfl;

gdtoa/gmisc.c

@@ -60,7 +60,7 @@ rshift(Bigint *b, int k)
 			while(x < xe)
 				*x1++ = *x++;
 		}
-	if ((b->wds = (int)(x1 - b->x)) == 0)
+	if ((b->wds = x1 - b->x) == 0)
 		b->x[0] = 0;
 	}
 

gdtoa/hd_init.c

@@ -31,6 +31,7 @@ THIS SOFTWARE.
 
 #include "gdtoaimp.h"
 
+#if 0
  unsigned char hexdig[256];
 
  static void
@@ -46,10 +47,31 @@ htinit(unsigned char *h, unsigned char *s, int inc)
 	}
 
  void
-hexdig_init_D2A(Void)
+hexdig_init_D2A(Void)	/* Use of hexdig_init omitted 20121220 to avoid a */
+			/* race condition when multiple threads are used. */
 {
 #define USC (unsigned char *)
 	htinit(hexdig, USC "0123456789", 0x10);
 	htinit(hexdig, USC "abcdef", 0x10 + 10);
 	htinit(hexdig, USC "ABCDEF", 0x10 + 10);
 	}
+#else
+ unsigned char hexdig[256] = {
+	0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+	0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+	0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+	16,17,18,19,20,21,22,23,24,25,0,0,0,0,0,0,
+	0,26,27,28,29,30,31,0,0,0,0,0,0,0,0,0,
+	0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+	0,26,27,28,29,30,31,0,0,0,0,0,0,0,0,0,
+	0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+	0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+	0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+	0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+	0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+	0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+	0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+	0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+	0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
+	};
+#endif

gdtoa/hexnan.c

@@ -54,15 +54,14 @@ L_shift(ULong *x, ULong *x1, int i)
 hexnan(sp, fpi, x0)
 	CONST char **sp; FPI *fpi; ULong *x0;
 #else
-hexnan( CONST char **sp, CONST FPI *fpi, ULong *x0)
+hexnan( CONST char **sp, FPI *fpi, ULong *x0)
 #endif
 {
 	ULong c, h, *x, *x1, *xe;
 	CONST char *s;
 	int havedig, hd0, i, nbits;
 
-	if (!hexdig['0'])
+	/**** if (!hexdig['0']) hexdig_init_D2A(); ****/
-		hexdig_init_D2A();
 	nbits = fpi->nbits;
 	x = x0 + (nbits >> kshift);
 	if (nbits & kmask)
@@ -72,12 +71,15 @@ hexnan( CONST char **sp, CONST FPI *fpi, ULong *x0)
 	havedig = hd0 = i = 0;
 	s = *sp;
 	/* allow optional initial 0x or 0X */
-	while((c = *(CONST unsigned char*)(s+1)) && c <= ' ')
+	while((c = *(CONST unsigned char*)(s+1)) && c <= ' ') {
+		if (!c)
+			goto retnan;
 		++s;
+		}
 	if (s[1] == '0' && (s[2] == 'x' || s[2] == 'X')
 	 && *(CONST unsigned char*)(s+3) > ' ')
 		s += 2;
-	while(c = *(CONST unsigned char*)++s) {
+	while((c = *(CONST unsigned char*)++s)) {
 		if (!(h = hexdig[c])) {
 			if (c <= ' ') {
 				if (hd0 < havedig) {
@@ -92,8 +94,11 @@ hexnan( CONST char **sp, CONST FPI *fpi, ULong *x0)
 					x1 = x;
 					i = 0;
 					}
-				while(*(CONST unsigned char*)(s+1) <= ' ')
+				while((c = *(CONST unsigned char*)(s+1)) <= ' ') {
+					if (!c)
+						goto retnan;
 					++s;
+					}
 				if (s[1] == '0' && (s[2] == 'x' || s[2] == 'X')
 				 && *(CONST unsigned char*)(s+3) > ' ')
 					s += 2;
@@ -107,10 +112,11 @@ hexnan( CONST char **sp, CONST FPI *fpi, ULong *x0)
 			do {
 				if (/*(*/ c == ')') {
 					*sp = s + 1;
-					break;
+					goto break2;
 					}
-				} while(c = *++s);
+				} while((c = *++s));
 #endif
+ retnan:
 			return STRTOG_NaN;
 			}
 		havedig++;
@@ -120,8 +126,11 @@ hexnan( CONST char **sp, CONST FPI *fpi, ULong *x0)
 			i = 1;
 			*--x = 0;
 			}
-		*x = (*x << 4) | h & 0xf;
+		*x = (*x << 4) | (h & 0xf);
 		}
+#ifndef GDTOA_NON_PEDANTIC_NANCHECK
+ break2:
+#endif
 	if (!havedig)
 		return STRTOG_NaN;
 	if (x < x1 && i < 8)

gdtoa/misc.c

@@ -38,11 +38,6 @@ THIS SOFTWARE.
 #endif
 #define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double))
 static double private_mem[PRIVATE_mem], *pmem_next = private_mem;
-#endif
-
-#ifdef MULTIPLE_THREADS
-static void ACQUIRE_DTOA_LOCK(int n);
-static void FREE_DTOA_LOCK(int n);
 #endif
 
  Bigint *
@@ -60,7 +55,9 @@ Balloc
 #endif
 
 	ACQUIRE_DTOA_LOCK(0);
-	if ( (rv = freelist[k]) !=0) {
+	/* The k > Kmax case does not need ACQUIRE_DTOA_LOCK(0), */
+	/* but this case seems very unlikely. */
+	if (k <= Kmax && (rv = freelist[k]) !=0) {
 		freelist[k] = rv->next;
 		}
 	else {
@@ -70,7 +67,7 @@ Balloc
 #else
 		len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
 			/sizeof(double);
-		if (pmem_next - private_mem + len <= PRIVATE_mem) {
+		if (k <= Kmax && pmem_next - private_mem + len <= PRIVATE_mem) {
 			rv = (Bigint*)pmem_next;
 			pmem_next += len;
 			}
@@ -94,12 +91,20 @@ Bfree
 #endif
 {
 	if (v) {
+		if (v->k > Kmax)
+#ifdef FREE
+			FREE((void*)v);
+#else
+			free((void*)v);
+#endif
+		else {
 			ACQUIRE_DTOA_LOCK(0);
 			v->next = freelist[v->k];
 			freelist[v->k] = v;
 			FREE_DTOA_LOCK(0);
 			}
 		}
+	}
 
  int
 lo0bits
@@ -109,8 +114,8 @@ lo0bits
 	(ULong *y)
 #endif
 {
-	register int k;
+	int k;
-	register ULong x = *y;
+	ULong x = *y;
 
 	if (x & 7) {
 		if (x & 1)
@@ -209,12 +214,12 @@ multadd
  int
 hi0bits_D2A
 #ifdef KR_headers
-	(x) register ULong x;
+	(x) ULong x;
 #else
-	(register ULong x)
+	(ULong x)
 #endif
 {
-	register int k = 0;
+	int k = 0;
 
 	if (!(x & 0xffff0000)) {
 		k = 16;
@@ -621,8 +626,8 @@ b2d
 #ifdef VAX
 	ULong d0, d1;
 #else
-#define d0 word0(d)
+#define d0 word0(&d)
-#define d1 word1(d)
+#define d1 word1(&d)
 #endif
 
 	xa0 = a->x;
@@ -635,16 +640,16 @@ b2d
 	*e = 32 - k;
 #ifdef Pack_32
 	if (k < Ebits) {
-		d0 = Exp_1 | y >> Ebits - k;
+		d0 = Exp_1 | y >> (Ebits - k);
 		w = xa > xa0 ? *--xa : 0;
-		d1 = y << (32-Ebits) + k | w >> Ebits - k;
+		d1 = y << ((32-Ebits) + k) | w >> (Ebits - k);
 		goto ret_d;
 		}
 	z = xa > xa0 ? *--xa : 0;
 	if (k -= Ebits) {
-		d0 = Exp_1 | y << k | z >> 32 - k;
+		d0 = Exp_1 | y << k | z >> (32 - k);
 		y = xa > xa0 ? *--xa : 0;
-		d1 = z << k | y >> 32 - k;
+		d1 = z << k | y >> (32 - k);
 		}
 	else {
 		d0 = Exp_1 | y;
@@ -668,10 +673,10 @@ b2d
 #endif
  ret_d:
 #ifdef VAX
-	word0(d) = d0 >> 16 | d0 << 16;
+	word0(&d) = d0 >> 16 | d0 << 16;
-	word1(d) = d1 >> 16 | d1 << 16;
+	word1(&d) = d1 >> 16 | d1 << 16;
 #endif
-	return dval(d);
+	return dval(&d);
 	}
 #undef d0
 #undef d1
@@ -679,29 +684,28 @@ b2d
  Bigint *
 d2b
 #ifdef KR_headers
-	(d, e, bits) double d; int *e, *bits;
+	(dd, e, bits) double dd; int *e, *bits;
 #else
-	(double _d, int *e, int *bits)
+	(double dd, int *e, int *bits)
 #endif
 {
 	Bigint *b;
+	U d;
 #ifndef Sudden_Underflow
 	int i;
 #endif
 	int de, k;
 	ULong *x, y, z;
-	U d;
 #ifdef VAX
 	ULong d0, d1;
 #else
-#define d0 word0(d)
+#define d0 word0(&d)
-#define d1 word1(d)
+#define d1 word1(&d)
 #endif
-
+	d.d = dd;
-	dval(d) = _d;
 #ifdef VAX
-	d0 = word0(d) >> 16 | word0(d) << 16;
+	d0 = word0(&d) >> 16 | word0(&d) << 16;
-	d1 = word1(d) >> 16 | word1(d) << 16;
+	d1 = word1(&d) >> 16 | word1(&d) << 16;
 #endif
 
 #ifdef Pack_32
@@ -725,7 +729,7 @@ d2b
 #ifdef Pack_32
 	if ( (y = d1) !=0) {
 		if ( (k = lo0bits(&y)) !=0) {
-			x[0] = y | z << 32 - k;
+			x[0] = y | z << (32 - k);
 			z >>= k;
 			}
 		else
@@ -736,10 +740,6 @@ d2b
 		     b->wds = (x[1] = z) !=0 ? 2 : 1;
 		}
 	else {
-#ifdef DEBUG
-		if (!z)
-			Bug("Zero passed to d2b");
-#endif
 		k = lo0bits(&z);
 		x[0] = z;
 #ifndef Sudden_Underflow
@@ -798,7 +798,7 @@ d2b
 #endif
 #ifdef IBM
 		*e = (de - Bias - (P-1) << 2) + k;
-		*bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
+		*bits = 4*P + 8 - k - hi0bits(word0(&d) & Frac_mask);
 #else
 		*e = de - Bias - (P-1) + k;
 		*bits = P - k;
@@ -851,71 +851,25 @@ strcp_D2A(a, b) char *a; char *b;
 strcp_D2A(char *a, CONST char *b)
 #endif
 {
-	while(*a = *b++)
+	while((*a = *b++))
 		a++;
 	return a;
 	}
 
 #ifdef NO_STRING_H
 
- KR_VOID *
+ Char *
 #ifdef KR_headers
-memcpy_D2A(a, b, len) char *a; char *b; size_t len;
+memcpy_D2A(a, b, len) Char *a; Char *b; size_t len;
 #else
 memcpy_D2A(void *a1, void *b1, size_t len)
 #endif
 {
-	register char *a = (char*)a1, *ae = a + len;
+	char *a = (char*)a1, *ae = a + len;
-	register char *b = (char*)b1, *a0 = a;
+	char *b = (char*)b1, *a0 = a;
 	while(a < ae)
 		*a++ = *b++;
 	return a0;
 	}
 
 #endif /* NO_STRING_H */
-
-
-#ifdef MULTIPLE_THREADS
-#ifdef _WIN32
-
-#undef Bias
-#define WIN32_LEAN_AND_MEAN
-#include <windows.h>
-
-static CRITICAL_SECTION dtoa_lock[2];
-static int did_init;
-
-void ACQUIRE_DTOA_LOCK(int n)
-{
-	if (!did_init)
-	{
-		did_init = 1;
-		InitializeCriticalSection(&dtoa_lock[0]);
-		InitializeCriticalSection(&dtoa_lock[1]);
-	}
-	EnterCriticalSection(&dtoa_lock[n]);
-}
-
-void FREE_DTOA_LOCK(int n)
-{
-	LeaveCriticalSection(&dtoa_lock[n]);
-}
-
-#else
-
-#include <pthread.h>
-
-static pthread_mutex_t dtoa_lock[2] = { PTHREAD_MUTEX_INITIALIZER, PTHREAD_MUTEX_INITIALIZER };
-
-void ACQUIRE_DTOA_LOCK(int n)
-{
-	pthread_mutex_lock(&dtoa_lock[n]);
-}
-
-void FREE_DTOA_LOCK(int n)
-{
-	pthread_mutex_unlock(&dtoa_lock[n]);
-}
-
-#endif
-#endif

gdtoa/qnan.c

@@ -55,11 +55,19 @@ typedef unsigned Long Ulong;
 #ifdef IEEE_8087
 #define _0 1
 #define _1 0
+#define _3 3
+#if defined(Gen_ld_QNAN) && !defined(NO_LONG_LONG)
+static int perm[4] = { 0, 1, 2, 3 };
+#endif
 #define HAVE_IEEE
 #endif
 #ifdef IEEE_MC68k
 #define _0 0
 #define _1 1
+#define _3 0
+#if defined(Gen_ld_QNAN) && !defined(NO_LONG_LONG)
+static int perm[4] = { 3, 2, 1, 0 };
+#endif
 #define HAVE_IEEE
 #endif
 
@@ -79,31 +87,32 @@ main(void)
 #endif
 		} U;
 	U a, b, c;
+#if defined(Gen_ld_QNAN) && !defined(NO_LONG_LONG)
 	int i;
+#endif
 
 	a.L[0] = b.L[0] = 0x7f800000;
 	c.f = a.f - b.f;
-	printf("#define f_QNAN 0x%lx\n", UL c.L[0]);
+	printf("#define f_QNAN 0x%lx\n", UL (c.L[0] & 0x7fffffff));
 	a.L[_0] = b.L[_0] = 0x7ff00000;
 	a.L[_1] = b.L[_1] = 0;
 	c.d = a.d - b.d;	/* quiet NaN */
-	printf("#define d_QNAN0 0x%lx\n", UL c.L[0]);
+	c.L[_0] &= 0x7fffffff;
-	printf("#define d_QNAN1 0x%lx\n", UL c.L[1]);
+	printf("#define d_QNAN0 0x%lx\n", UL c.L[_0]);
-#ifdef NO_LONG_LONG
+	printf("#define d_QNAN1 0x%lx\n", UL c.L[_1]);
-	for(i = 0; i < 4; i++)
+#ifndef NO_LONG_LONG
-		printf("#define ld_QNAN%d 0xffffffff\n", i);
+#ifdef Gen_ld_QNAN
-	for(i = 0; i < 5; i++)
+	if (sizeof(a.D) >= 16) {
-		printf("#define ldus_QNAN%d 0xffff\n", i);
-#else
 		b.D = c.D = a.d;
 		if (printf("") < 0)
 			c.D = 37;	/* never executed; just defeat optimization */
-	a.L[2] = a.L[3] = 0;
+		a.L[0] = a.L[1] = a.L[2] = a.L[3] = 0;
 		a.D = b.D - c.D;
+		a.L[_3] &= 0x7fffffff;
 		for(i = 0; i < 4; i++)
-		printf("#define ld_QNAN%d 0x%lx\n", i, UL a.L[i]);
+			printf("#define ld_QNAN%d 0x%lx\n", i, UL a.L[perm[i]]);
-	for(i = 0; i < 5; i++)
+		}
-		printf("#define ldus_QNAN%d 0x%x\n", i, a.u[i]);
+#endif
 #endif
 #endif /* HAVE_IEEE */
 	return 0;

gdtoa/smisc.c

@@ -34,9 +34,9 @@ THIS SOFTWARE.
  Bigint *
 s2b
 #ifdef KR_headers
-	(s, nd0, nd, y9) CONST char *s; int nd0, nd; ULong y9;
+	(s, nd0, nd, y9, dplen) CONST char *s; int dplen, nd0, nd; ULong y9;
 #else
-	(CONST char *s, int nd0, int nd, ULong y9)
+	(CONST char *s, int nd0, int nd, ULong y9, int dplen)
 #endif
 {
 	Bigint *b;
@@ -60,10 +60,10 @@ s2b
 		s += 9;
 		do b = multadd(b, 10, *s++ - '0');
 			while(++i < nd0);
-		s++;
+		s += dplen;
 		}
 	else
-		s += 10;
+		s += dplen + 9;
 	for(; i < nd; i++)
 		b = multadd(b, 10, *s++ - '0');
 	return b;
@@ -80,30 +80,30 @@ ratio
 	U da, db;
 	int k, ka, kb;
 
-	dval(da) = b2d(a, &ka);
+	dval(&da) = b2d(a, &ka);
-	dval(db) = b2d(b, &kb);
+	dval(&db) = b2d(b, &kb);
 	k = ka - kb + ULbits*(a->wds - b->wds);
 #ifdef IBM
 	if (k > 0) {
-		word0(da) += (k >> 2)*Exp_msk1;
+		word0(&da) += (k >> 2)*Exp_msk1;
 		if (k &= 3)
-			dval(da) *= 1 << k;
+			dval(&da) *= 1 << k;
 		}
 	else {
 		k = -k;
-		word0(db) += (k >> 2)*Exp_msk1;
+		word0(&db) += (k >> 2)*Exp_msk1;
 		if (k &= 3)
-			dval(db) *= 1 << k;
+			dval(&db) *= 1 << k;
 		}
 #else
 	if (k > 0)
-		word0(da) += k*Exp_msk1;
+		word0(&da) += k*Exp_msk1;
 	else {
 		k = -k;
-		word0(db) += k*Exp_msk1;
+		word0(&db) += k*Exp_msk1;
 		}
 #endif
-	return dval(da) / dval(db);
+	return dval(&da) / dval(&db);
 	}
 
 #ifdef INFNAN_CHECK

gdtoa/strtoIQ.c

@@ -38,7 +38,7 @@ strtoIQ(s, sp, a, b) CONST char *s; char **sp; void *a; void *b;
 strtoIQ(CONST char *s, char **sp, void *a, void *b)
 #endif
 {
-	static CONST FPI fpi = { 113, 1-16383-113+1, 32766-16383-113+1, 1, SI };
+	static FPI fpi = { 113, 1-16383-113+1, 32766-16383-113+1, 1, SI };
 	Long exp[2];
 	Bigint *B[2];
 	int k, rv[2];

gdtoa/strtoId.c

@@ -38,7 +38,7 @@ strtoId(s, sp, f0, f1) CONST char *s; char **sp; double *f0, *f1;
 strtoId(CONST char *s, char **sp, double *f0, double *f1)
 #endif
 {
-	static CONST 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[2];
 	Bigint *B[2];
 	int k, rv[2];

gdtoa/strtoIdd.c

@@ -39,9 +39,9 @@ strtoIdd(CONST char *s, char **sp, double *f0, double *f1)
 #endif
 {
 #ifdef Sudden_Underflow
-	static CONST FPI fpi = { 106, 1-1023, 2046-1023-106+1, 1, 1 };
+	static FPI fpi = { 106, 1-1023, 2046-1023-106+1, 1, 1 };
 #else
-	static CONST FPI fpi = { 106, 1-1023-53+1, 2046-1023-106+1, 1, 0 };
+	static FPI fpi = { 106, 1-1023-53+1, 2046-1023-106+1, 1, 0 };
 #endif
 	Long exp[2];
 	Bigint *B[2];

gdtoa/strtoIf.c

@@ -38,7 +38,7 @@ strtoIf(s, sp, f0, f1) CONST char *s; char **sp; float *f0, *f1;
 strtoIf(CONST char *s, char **sp, float *f0, float *f1)
 #endif
 {
-	static CONST FPI fpi = { 24, 1-127-24+1,  254-127-24+1, 1, SI };
+	static FPI fpi = { 24, 1-127-24+1,  254-127-24+1, 1, SI };
 	Long exp[2];
 	Bigint *B[2];
 	int k, rv[2];

gdtoa/strtoIg.c

@@ -35,7 +35,7 @@ THIS SOFTWARE.
 #ifdef KR_headers
 strtoIg(s00, se, fpi, exp, B, rvp) CONST char *s00; char **se; FPI *fpi; Long *exp; Bigint **B; int *rvp;
 #else
-strtoIg(CONST char *s00, char **se, CONST FPI *fpi, Long *exp, Bigint **B, int *rvp)
+strtoIg(CONST char *s00, char **se, FPI *fpi, Long *exp, Bigint **B, int *rvp)
 #endif
 {
 	Bigint *b, *b1;
@@ -61,16 +61,20 @@ strtoIg(CONST char *s00, char **se, CONST FPI *fpi, Long *exp, Bigint **B, int *
 	if (rv & STRTOG_Inexlo) {
 		swap = 0;
 		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[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
-		 || nb1 && b1->x[nw1] & 1L << nb1) {
+		 || (nb1 && b1->x[nw1] & 1L << nb1)) {
 			if (++e1 > fpi->emax)
 				rv1 = STRTOG_Infinite | STRTOG_Inexhi;
 			rshift(b1, 1);

gdtoa/strtoIx.c

@@ -38,7 +38,7 @@ strtoIx(s, sp, a, b) CONST char *s; char **sp; void *a; void *b;
 strtoIx(CONST char *s, char **sp, void *a, void *b)
 #endif
 {
-	static CONST FPI fpi = { 64, 1-16383-64+1, 32766 - 16383 - 64 + 1, 1, SI };
+	static FPI fpi = { 64, 1-16383-64+1, 32766 - 16383 - 64 + 1, 1, SI };
 	Long exp[2];
 	Bigint *B[2];
 	int k, rv[2];

gdtoa/strtoIxL.c

@@ -38,7 +38,7 @@ strtoIxL(s, sp, a, b) CONST char *s; char **sp; void *a; void *b;
 strtoIxL(CONST char *s, char **sp, void *a, void *b)
 #endif
 {
-	static CONST FPI fpi = { 64, 1-16383-64+1, 32766 - 16383 - 64 + 1, 1, SI };
+	static FPI fpi = { 64, 1-16383-64+1, 32766 - 16383 - 64 + 1, 1, SI };
 	Long exp[2];
 	Bigint *B[2];
 	int k, rv[2];

gdtoa/strtod.c

@@ -42,24 +42,43 @@ 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
 
-#include <float.h>
-
 #ifdef Honor_FLT_ROUNDS
-#define Rounding rounding
 #undef Check_FLT_ROUNDS
 #define Check_FLT_ROUNDS
 #else
 #define Rounding Flt_Rounds
 #endif
 
+#ifdef Avoid_Underflow /*{*/
+ static double
+sulp
+#ifdef KR_headers
+	(x, scale) U *x; int scale;
+#else
+	(U *x, int scale)
+#endif
+{
+	U u;
+	double rv;
+	int i;
+
+	rv = ulp(x);
+	if (!scale || (i = 2*P + 1 - ((word0(x) & Exp_mask) >> Exp_shift)) <= 0)
+		return rv; /* Is there an example where i <= 0 ? */
+	word0(&u) = Exp_1 + (i << Exp_shift);
+	word1(&u) = 0;
+	return rv * u.d;
+	}
+#endif /*}*/
+
  double
 strtod
 #ifdef KR_headers
@@ -74,21 +93,55 @@ strtod
 	int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, decpt, dsign,
 		 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
 	CONST char *s, *s0, *s1;
-	double aadj, adj;
+	double aadj;
-	U rv, rv0, aadj1;
 	Long L;
+	U adj, aadj1, rv, rv0;
 	ULong y, z;
 	Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
+#ifdef Avoid_Underflow
+	ULong Lsb, Lsb1;
+#endif
 #ifdef SET_INEXACT
 	int inexact, oldinexact;
 #endif
-#ifdef Honor_FLT_ROUNDS
+#ifdef USE_LOCALE /*{{*/
-	int rounding;
+#ifdef NO_LOCALE_CACHE
-#endif
+	char *decimalpoint = localeconv()->decimal_point;
+	int dplen = strlen(decimalpoint);
+#else
+	char *decimalpoint;
+	static char *decimalpoint_cache;
+	static int dplen;
+	if (!(s0 = decimalpoint_cache)) {
+		s0 = localeconv()->decimal_point;
+		if ((decimalpoint_cache = (char*)MALLOC(strlen(s0) + 1))) {
+			strcpy(decimalpoint_cache, s0);
+			s0 = decimalpoint_cache;
+			}
+		dplen = strlen(s0);
+		}
+	decimalpoint = (char*)s0;
+#endif /*NO_LOCALE_CACHE*/
+#else  /*USE_LOCALE}{*/
+#define dplen 1
+#endif /*USE_LOCALE}}*/
+
+#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 /*}*/
 
-	//_control87(_PC_53, _MCW_PC);
 	sign = nz0 = nz = decpt = 0;
-	dval(rv) = 0.;
+	dval(&rv) = 0.;
 	for(s = s00;;s++) switch(*s) {
 		case '-':
 			sign = 1;
@@ -111,22 +164,18 @@ strtod
 		}
  break2:
 	if (*s == '0') {
-#ifndef NO_HEX_FP
+#ifndef NO_HEX_FP /*{*/
 		{
-		static CONST 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;
 		ULong bits[2];
 		switch(s[1]) {
 		  case 'x':
 		  case 'X':
 			{
-#if defined(FE_DOWNWARD) && defined(FE_TONEAREST) && defined(FE_TOWARDZERO) && defined(FE_UPWARD)
+#ifdef Honor_FLT_ROUNDS
 			FPI fpi1 = fpi;
-			switch(fegetround()) {
+			fpi1.rounding = Rounding;
-			  case FE_TOWARDZERO:	fpi1.rounding = 0; break;
-			  case FE_UPWARD:	fpi1.rounding = 2; break;
-			  case FE_DOWNWARD:	fpi1.rounding = 3;
-			  }
 #else
 #define fpi1 fpi
 #endif
@@ -146,7 +195,7 @@ strtod
 			goto ret;
 		  }
 		}
-#endif
+#endif /*}*/
 		nz0 = 1;
 		while(*++s == '0') ;
 		if (!*s)
@@ -161,13 +210,17 @@ strtod
 			z = 10*z + c - '0';
 	nd0 = nd;
 #ifdef USE_LOCALE
-	if (c == *localeconv()->decimal_point)
+	if (c == *decimalpoint) {
+		for(i = 1; decimalpoint[i]; ++i)
+			if (s[i] != decimalpoint[i])
+				goto dig_done;
+		s += i;
+		c = *s;
 #else
-	if (c == '.')
+	if (c == '.') {
-#endif
-		{
-		decpt = 1;
 		c = *++s;
+#endif
+		decpt = 1;
 		if (!nd) {
 			for(; c == '0'; c = *++s)
 				nz++;
@@ -196,7 +249,7 @@ strtod
 				nz = 0;
 				}
 			}
-		}
+		}/*}*/
  dig_done:
 	e = 0;
 	if (c == 'e' || c == 'E') {
@@ -240,7 +293,7 @@ strtod
 #ifdef INFNAN_CHECK
 			/* Check for Nan and Infinity */
 			ULong bits[2];
-			static CONST FPI fpinan =	/* only 52 explicit bits */
+			static FPI fpinan =	/* only 52 explicit bits */
 				{ 52, 1-1023-53+1, 2046-1023-53+1, 1, SI };
 			if (!decpt)
 			 switch(c) {
@@ -250,8 +303,8 @@ strtod
 					--s;
 					if (!match(&s,"inity"))
 						++s;
-					word0(rv) = 0x7ff00000;
+					word0(&rv) = 0x7ff00000;
-					word1(rv) = 0;
+					word1(&rv) = 0;
 					goto ret;
 					}
 				break;
@@ -262,13 +315,13 @@ strtod
 					if (*s == '(' /*)*/
 					 && hexnan(&s, &fpinan, bits)
 							== STRTOG_NaNbits) {
-						word0(rv) = 0x7ff00000 | bits[1];
+						word0(&rv) = 0x7ff00000 | bits[1];
-						word1(rv) = bits[0];
+						word1(&rv) = bits[0];
 						}
 					else {
 #endif
-						word0(rv) = NAN_WORD0;
+						word0(&rv) = NAN_WORD0;
-						word1(rv) = NAN_WORD1;
+						word1(&rv) = NAN_WORD1;
 #ifndef No_Hex_NaN
 						}
 #endif
@@ -292,13 +345,13 @@ strtod
 	if (!nd0)
 		nd0 = nd;
 	k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
-	dval(rv) = y;
+	dval(&rv) = y;
 	if (k > 9) {
 #ifdef SET_INEXACT
 		if (k > DBL_DIG)
 			oldinexact = get_inexact();
 #endif
-		dval(rv) = tens[k - 9] * dval(rv) + z;
+		dval(&rv) = tens[k - 9] * dval(&rv) + z;
 		}
 	bd0 = 0;
 	if (nd <= DBL_DIG
@@ -310,6 +363,7 @@ strtod
 			) {
 		if (!e)
 			goto ret;
+#ifndef ROUND_BIASED_without_Round_Up
 		if (e > 0) {
 			if (e <= Ten_pmax) {
 #ifdef VAX
@@ -318,11 +372,11 @@ strtod
 #ifdef Honor_FLT_ROUNDS
 				/* round correctly FLT_ROUNDS = 2 or 3 */
 				if (sign) {
-					rv = -rv;
+					rv.d = -rv.d;
 					sign = 0;
 					}
 #endif
-				/* rv = */ rounded_product(dval(rv), tens[e]);
+				/* rv = */ rounded_product(dval(&rv), tens[e]);
 				goto ret;
 #endif
 				}
@@ -334,25 +388,25 @@ strtod
 #ifdef Honor_FLT_ROUNDS
 				/* round correctly FLT_ROUNDS = 2 or 3 */
 				if (sign) {
-					rv = -rv;
+					rv.d = -rv.d;
 					sign = 0;
 					}
 #endif
 				e -= i;
-				dval(rv) *= tens[i];
+				dval(&rv) *= tens[i];
 #ifdef VAX
 				/* VAX exponent range is so narrow we must
 				 * worry about overflow here...
 				 */
  vax_ovfl_check:
-				word0(rv) -= P*Exp_msk1;
+				word0(&rv) -= P*Exp_msk1;
-				/* rv = */ rounded_product(dval(rv), tens[e]);
+				/* rv = */ rounded_product(dval(&rv), tens[e]);
-				if ((word0(rv) & Exp_mask)
+				if ((word0(&rv) & Exp_mask)
 				 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
 					goto ovfl;
-				word0(rv) += P*Exp_msk1;
+				word0(&rv) += P*Exp_msk1;
 #else
-				/* rv = */ rounded_product(dval(rv), tens[e]);
+				/* rv = */ rounded_product(dval(&rv), tens[e]);
 #endif
 				goto ret;
 				}
@@ -362,14 +416,15 @@ strtod
 #ifdef Honor_FLT_ROUNDS
 			/* round correctly FLT_ROUNDS = 2 or 3 */
 			if (sign) {
-				rv = -rv;
+				rv.d = -rv.d;
 				sign = 0;
 				}
 #endif
-			/* rv = */ rounded_quotient(dval(rv), tens[-e]);
+			/* rv = */ rounded_quotient(dval(&rv), tens[-e]);
 			goto ret;
 			}
 #endif
+#endif /* ROUND_BIASED_without_Round_Up */
 		}
 	e1 += nd - k;
 
@@ -383,12 +438,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)
+			if (Rounding != 2)
-				rounding = 0;
+				Rounding = 0;
 		}
 #endif
 #endif /*IEEE_Arith*/
@@ -397,67 +452,73 @@ strtod
 
 	if (e1 > 0) {
 		if ( (i = e1 & 15) !=0)
-			dval(rv) *= tens[i];
+			dval(&rv) *= tens[i];
 		if (e1 &= ~15) {
 			if (e1 > DBL_MAX_10_EXP) {
  ovfl:
-#ifndef NO_ERRNO
-				errno = ERANGE;
-#endif
 				/* 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;
+					word0(&rv) = Big0;
-					word1(rv) = Big1;
+					word1(&rv) = Big1;
 					break;
 				  default:
-					word0(rv) = Exp_mask;
+					word0(&rv) = Exp_mask;
-					word1(rv) = 0;
+					word1(&rv) = 0;
 				  }
 #else /*Honor_FLT_ROUNDS*/
-				word0(rv) = Exp_mask;
+				word0(&rv) = Exp_mask;
-				word1(rv) = 0;
+				word1(&rv) = 0;
 #endif /*Honor_FLT_ROUNDS*/
 #ifdef SET_INEXACT
 				/* set overflow bit */
-				dval(rv0) = 1e300;
+				dval(&rv0) = 1e300;
-				dval(rv0) *= dval(rv0);
+				dval(&rv0) *= dval(&rv0);
 #endif
 #else /*IEEE_Arith*/
-				word0(rv) = Big0;
+				word0(&rv) = Big0;
-				word1(rv) = Big1;
+				word1(&rv) = Big1;
 #endif /*IEEE_Arith*/
-				if (bd0)
+ range_err:
-					goto retfree;
+				if (bd0) {
+					Bfree(bb);
+					Bfree(bd);
+					Bfree(bs);
+					Bfree(bd0);
+					Bfree(delta);
+					}
+#ifndef NO_ERRNO
+				errno = ERANGE;
+#endif
 				goto ret;
 				}
 			e1 >>= 4;
 			for(j = 0; e1 > 1; j++, e1 >>= 1)
 				if (e1 & 1)
-					dval(rv) *= bigtens[j];
+					dval(&rv) *= bigtens[j];
 		/* The last multiplication could overflow. */
-			word0(rv) -= P*Exp_msk1;
+			word0(&rv) -= P*Exp_msk1;
-			dval(rv) *= bigtens[j];
+			dval(&rv) *= bigtens[j];
-			if ((z = word0(rv) & Exp_mask)
+			if ((z = word0(&rv) & Exp_mask)
 			 > Exp_msk1*(DBL_MAX_EXP+Bias-P))
 				goto ovfl;
 			if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
 				/* set to largest number */
 				/* (Can't trust DBL_MAX) */
-				word0(rv) = Big0;
+				word0(&rv) = Big0;
-				word1(rv) = Big1;
+				word1(&rv) = Big1;
 				}
 			else
-				word0(rv) += P*Exp_msk1;
+				word0(&rv) += P*Exp_msk1;
 			}
 		}
 	else if (e1 < 0) {
 		e1 = -e1;
 		if ( (i = e1 & 15) !=0)
-			dval(rv) /= tens[i];
+			dval(&rv) /= tens[i];
 		if (e1 >>= 4) {
 			if (e1 >= 1 << n_bigtens)
 				goto undfl;
@@ -466,44 +527,39 @@ strtod
 				scale = 2*P;
 			for(j = 0; e1 > 0; j++, e1 >>= 1)
 				if (e1 & 1)
-					dval(rv) *= tinytens[j];
+					dval(&rv) *= tinytens[j];
-			if (scale && (j = 2*P + 1 - ((word0(rv) & Exp_mask)
+			if (scale && (j = 2*P + 1 - ((word0(&rv) & Exp_mask)
 						>> Exp_shift)) > 0) {
 				/* scaled rv is denormal; zap j low bits */
 				if (j >= 32) {
-					word1(rv) = 0;
+					word1(&rv) = 0;
 					if (j >= 53)
-					 word0(rv) = (P+2)*Exp_msk1;
+					 word0(&rv) = (P+2)*Exp_msk1;
 					else
-					 word0(rv) &= 0xffffffff << j-32;
+					 word0(&rv) &= 0xffffffff << (j-32);
 					}
 				else
-					word1(rv) &= 0xffffffff << j;
+					word1(&rv) &= 0xffffffff << j;
 				}
 #else
 			for(j = 0; e1 > 1; j++, e1 >>= 1)
 				if (e1 & 1)
-					dval(rv) *= tinytens[j];
+					dval(&rv) *= tinytens[j];
 			/* The last multiplication could underflow. */
-			dval(rv0) = dval(rv);
+			dval(&rv0) = dval(&rv);
-			dval(rv) *= tinytens[j];
+			dval(&rv) *= tinytens[j];
-			if (!dval(rv)) {
+			if (!dval(&rv)) {
-				dval(rv) = 2.*dval(rv0);
+				dval(&rv) = 2.*dval(&rv0);
-				dval(rv) *= tinytens[j];
+				dval(&rv) *= tinytens[j];
 #endif
-				if (!dval(rv)) {
+				if (!dval(&rv)) {
  undfl:
-					dval(rv) = 0.;
+					dval(&rv) = 0.;
-#ifndef NO_ERRNO
+					goto range_err;
-					errno = ERANGE;
-#endif
-					if (bd0)
-						goto retfree;
-					goto ret;
 					}
 #ifndef Avoid_Underflow
-				word0(rv) = Tiny0;
+				word0(&rv) = Tiny0;
-				word1(rv) = Tiny1;
+				word1(&rv) = Tiny1;
 				/* The refinement below will clean
 				 * this approximation up.
 				 */
@@ -516,12 +572,12 @@ strtod
 
 	/* Put digits into bd: true value = bd * 10^e */
 
-	bd0 = s2b(s0, nd0, nd, y);
+	bd0 = s2b(s0, nd0, nd, y, dplen);
 
 	for(;;) {
 		bd = Balloc(bd0->k);
 		Bcopy(bd, bd0);
-		bb = d2b(dval(rv), &bbe, &bbbits);	/* rv = bb * 2^bbe */
+		bb = d2b(dval(&rv), &bbe, &bbbits);	/* rv = bb * 2^bbe */
 		bs = i2b(1);
 
 		if (e >= 0) {
@@ -538,16 +594,23 @@ strtod
 			bd2 -= bbe;
 		bs2 = bb2;
 #ifdef Honor_FLT_ROUNDS
-		if (rounding != 1)
+		if (Rounding != 1)
 			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;
+		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
@@ -557,7 +620,7 @@ strtod
 #endif
 #else /*Sudden_Underflow*/
 		j = bbe;
-		i = j + bbbits - 1;	/* logb(rv) */
+		i = j + bbbits - 1;	/* logb(&rv) */
 		if (i < Emin)	/* denormal */
 			j += P - Emin;
 		else
@@ -596,7 +659,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) {
@@ -606,17 +669,17 @@ strtod
 #endif
 					break;
 					}
-				if (rounding) {
+				if (Rounding) {
 					if (dsign) {
-						adj = 1.;
+						dval(&adj) = 1.;
 						goto apply_adj;
 						}
 					}
 				else if (!dsign) {
-					adj = -1.;
+					dval(&adj) = -1.;
-					if (!word1(rv)
+					if (!word1(&rv)
-					 && !(word0(rv) & Frac_mask)) {
+					 && !(word0(&rv) & Frac_mask)) {
-						y = word0(rv) & Exp_mask;
+						y = word0(&rv) & Exp_mask;
 #ifdef Avoid_Underflow
 						if (!scale || y > 2*P*Exp_msk1)
 #else
@@ -625,63 +688,66 @@ strtod
 						  {
 						  delta = lshift(delta,Log2P);
 						  if (cmp(delta, bs) <= 0)
-							adj = -0.5;
+							dval(&adj) = -0.5;
 						  }
 						}
  apply_adj:
 #ifdef Avoid_Underflow
-					if (scale && (y = word0(rv) & Exp_mask)
+					if (scale && (y = word0(&rv) & Exp_mask)
 						<= 2*P*Exp_msk1)
-					  word0(adj) += (2*P+1)*Exp_msk1 - y;
+					  word0(&adj) += (2*P+1)*Exp_msk1 - y;
 #else
 #ifdef Sudden_Underflow
-					if ((word0(rv) & Exp_mask) <=
+					if ((word0(&rv) & Exp_mask) <=
 							P*Exp_msk1) {
-						word0(rv) += P*Exp_msk1;
+						word0(&rv) += P*Exp_msk1;
-						dval(rv) += adj*ulp(dval(rv));
+						dval(&rv) += adj*ulp(&rv);
-						word0(rv) -= P*Exp_msk1;
+						word0(&rv) -= P*Exp_msk1;
 						}
 					else
 #endif /*Sudden_Underflow*/
 #endif /*Avoid_Underflow*/
-					dval(rv) += adj*ulp(dval(rv));
+					dval(&rv) += adj.d*ulp(&rv);
 					}
 				break;
 				}
-			adj = ratio(delta, bs);
+			dval(&adj) = ratio(delta, bs);
-			if (adj < 1.)
+			if (adj.d < 1.)
-				adj = 1.;
+				dval(&adj) = 1.;
-			if (adj <= 0x7ffffffe) {
+			if (adj.d <= 0x7ffffffe) {
-				/* adj = rounding ? ceil(adj) : floor(adj); */
+				/* dval(&adj) = Rounding ? ceil(&adj) : floor(&adj); */
-				y = adj;
+				y = adj.d;
-				if (y != adj) {
+				if (y != adj.d) {
-					if (!((rounding>>1) ^ dsign))
+					if (!((Rounding>>1) ^ dsign))
 						y++;
-					adj = y;
+					dval(&adj) = y;
 					}
 				}
 #ifdef Avoid_Underflow
-			if (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
+			if (scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1)
-				word0(adj) += (2*P+1)*Exp_msk1 - y;
+				word0(&adj) += (2*P+1)*Exp_msk1 - y;
 #else
 #ifdef Sudden_Underflow
-			if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
+			if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) {
-				word0(rv) += P*Exp_msk1;
+				word0(&rv) += P*Exp_msk1;
-				adj *= ulp(dval(rv));
+				dval(&adj) *= ulp(&rv);
 				if (dsign)
-					dval(rv) += adj;
+					dval(&rv) += adj;
 				else
-					dval(rv) -= adj;
+					dval(&rv) -= adj;
-				word0(rv) -= P*Exp_msk1;
+				word0(&rv) -= P*Exp_msk1;
 				goto cont;
 				}
 #endif /*Sudden_Underflow*/
 #endif /*Avoid_Underflow*/
-			adj *= ulp(dval(rv));
+			dval(&adj) *= ulp(&rv);
-			if (dsign)
+			if (dsign) {
-				dval(rv) += adj;
+				if (word0(&rv) == Big0 && word1(&rv) == Big1)
+					goto ovfl;
+				dval(&rv) += adj.d;
+				}
 			else
-				dval(rv) -= adj;
+				dval(&rv) -= adj.d;
 			goto cont;
 			}
 #endif /*Honor_FLT_ROUNDS*/
@@ -690,12 +756,12 @@ strtod
 			/* Error is less than half an ulp -- check for
 			 * special case of mantissa a power of two.
 			 */
-			if (dsign || word1(rv) || word0(rv) & Bndry_mask
+			if (dsign || word1(&rv) || word0(&rv) & Bndry_mask
 #ifdef IEEE_Arith
 #ifdef Avoid_Underflow
-			 || (word0(rv) & Exp_mask) <= (2*P+1)*Exp_msk1
+			 || (word0(&rv) & Exp_mask) <= (2*P+1)*Exp_msk1
 #else
-			 || (word0(rv) & Exp_mask) <= Exp_msk1
+			 || (word0(&rv) & Exp_mask) <= Exp_msk1
 #endif
 #endif
 				) {
@@ -720,32 +786,34 @@ strtod
 		if (i == 0) {
 			/* exactly half-way between */
 			if (dsign) {
-				if ((word0(rv) & Bndry_mask1) == Bndry_mask1
+				if ((word0(&rv) & Bndry_mask1) == Bndry_mask1
-				 &&  word1(rv) == (
+				 &&  word1(&rv) == (
 #ifdef Avoid_Underflow
-			(scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
+			(scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1)
 		? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
 #endif
 						   0xffffffff)) {
 					/*boundary case -- increment exponent*/
-					word0(rv) = (word0(rv) & Exp_mask)
+					if (word0(&rv) == Big0 && word1(&rv) == Big1)
+						goto ovfl;
+					word0(&rv) = (word0(&rv) & Exp_mask)
 						+ Exp_msk1
 #ifdef IBM
 						| Exp_msk1 >> 4
 #endif
 						;
-					word1(rv) = 0;
+					word1(&rv) = 0;
 #ifdef Avoid_Underflow
 					dsign = 0;
 #endif
 					break;
 					}
 				}
-			else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
+			else if (!(word0(&rv) & Bndry_mask) && !word1(&rv)) {
  drop_down:
 				/* boundary case -- decrement exponent */
 #ifdef Sudden_Underflow /*{{*/
-				L = word0(rv) & Exp_mask;
+				L = word0(&rv) & Exp_mask;
 #ifdef IBM
 				if (L <  Exp_msk1)
 #else
@@ -760,7 +828,7 @@ strtod
 #else /*Sudden_Underflow}{*/
 #ifdef Avoid_Underflow
 				if (scale) {
-					L = word0(rv) & Exp_mask;
+					L = word0(&rv) & Exp_mask;
 					if (L <= (2*P+1)*Exp_msk1) {
 						if (L > (P+2)*Exp_msk1)
 							/* round even ==> */
@@ -771,10 +839,10 @@ strtod
 						}
 					}
 #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
 				goto cont;
 #else
@@ -782,16 +850,33 @@ strtod
 #endif
 				}
 #ifndef ROUND_BIASED
-			if (!(word1(rv) & LSB))
+#ifdef Avoid_Underflow
+			if (Lsb1) {
+				if (!(word0(&rv) & Lsb1))
+					break;
+				}
+			else if (!(word1(&rv) & Lsb))
+				break;
+#else
+			if (!(word1(&rv) & LSB))
 				break;
+#endif
 #endif
 			if (dsign)
-				dval(rv) += ulp(dval(rv));
+#ifdef Avoid_Underflow
+				dval(&rv) += sulp(&rv, scale);
+#else
+				dval(&rv) += ulp(&rv);
+#endif
 #ifndef ROUND_BIASED
 			else {
-				dval(rv) -= ulp(dval(rv));
+#ifdef Avoid_Underflow
+				dval(&rv) -= sulp(&rv, scale);
+#else
+				dval(&rv) -= ulp(&rv);
+#endif
 #ifndef Sudden_Underflow
-				if (!dval(rv))
+				if (!dval(&rv))
 					goto undfl;
 #endif
 				}
@@ -803,14 +888,14 @@ strtod
 			}
 		if ((aadj = ratio(delta, bs)) <= 2.) {
 			if (dsign)
-				aadj = dval(aadj1) = 1.;
+				aadj = dval(&aadj1) = 1.;
-			else if (word1(rv) || word0(rv) & Bndry_mask) {
+			else if (word1(&rv) || word0(&rv) & Bndry_mask) {
 #ifndef Sudden_Underflow
-				if (word1(rv) == Tiny1 && !word0(rv))
+				if (word1(&rv) == Tiny1 && !word0(&rv))
 					goto undfl;
 #endif
 				aadj = 1.;
-				dval(aadj1) = -1.;
+				dval(&aadj1) = -1.;
 				}
 			else {
 				/* special case -- power of FLT_RADIX to be */
@@ -820,45 +905,45 @@ strtod
 					aadj = 1./FLT_RADIX;
 				else
 					aadj *= 0.5;
-				dval(aadj1) = -aadj;
+				dval(&aadj1) = -aadj;
 				}
 			}
 		else {
 			aadj *= 0.5;
-			dval(aadj1) = dsign ? aadj : -aadj;
+			dval(&aadj1) = dsign ? aadj : -aadj;
 #ifdef Check_FLT_ROUNDS
 			switch(Rounding) {
 				case 2: /* towards +infinity */
-					aadj1 -= 0.5;
+					dval(&aadj1) -= 0.5;
 					break;
 				case 0: /* towards 0 */
 				case 3: /* towards -infinity */
-					aadj1 += 0.5;
+					dval(&aadj1) += 0.5;
 				}
 #else
 			if (Flt_Rounds == 0)
-				dval(aadj1) += 0.5;
+				dval(&aadj1) += 0.5;
 #endif /*Check_FLT_ROUNDS*/
 			}
-		y = word0(rv) & Exp_mask;
+		y = word0(&rv) & Exp_mask;
 
 		/* Check for overflow */
 
 		if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
-			dval(rv0) = dval(rv);
+			dval(&rv0) = dval(&rv);
-			word0(rv) -= P*Exp_msk1;
+			word0(&rv) -= P*Exp_msk1;
-			adj = dval(aadj1) * ulp(dval(rv));
+			dval(&adj) = dval(&aadj1) * ulp(&rv);
-			dval(rv) += adj;
+			dval(&rv) += dval(&adj);
-			if ((word0(rv) & Exp_mask) >=
+			if ((word0(&rv) & Exp_mask) >=
 					Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
-				if (word0(rv0) == Big0 && word1(rv0) == Big1)
+				if (word0(&rv0) == Big0 && word1(&rv0) == Big1)
 					goto ovfl;
-				word0(rv) = Big0;
+				word0(&rv) = Big0;
-				word1(rv) = Big1;
+				word1(&rv) = Big1;
 				goto cont;
 				}
 			else
-				word0(rv) += P*Exp_msk1;
+				word0(&rv) += P*Exp_msk1;
 			}
 		else {
 #ifdef Avoid_Underflow
@@ -867,58 +952,58 @@ strtod
 					if ((z = (ULong)aadj) <= 0)
 						z = 1;
 					aadj = z;
-					dval(aadj1) = dsign ? aadj : -aadj;
+					dval(&aadj1) = dsign ? aadj : -aadj;
 					}
-				word0(aadj1) += (2*P+1)*Exp_msk1 - y;
+				word0(&aadj1) += (2*P+1)*Exp_msk1 - y;
 				}
-			adj = dval(aadj1) * ulp(dval(rv));
+			dval(&adj) = dval(&aadj1) * ulp(&rv);
-			dval(rv) += adj;
+			dval(&rv) += dval(&adj);
 #else
 #ifdef Sudden_Underflow
-			if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
+			if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) {
-				dval(rv0) = dval(rv);
+				dval(&rv0) = dval(&rv);
-				word0(rv) += P*Exp_msk1;
+				word0(&rv) += P*Exp_msk1;
-				adj = aadj1 * ulp(dval(rv));
+				dval(&adj) = dval(&aadj1) * ulp(&rv);
-				dval(rv) += adj;
+				dval(&rv) += adj;
 #ifdef IBM
-				if ((word0(rv) & Exp_mask) <  P*Exp_msk1)
+				if ((word0(&rv) & Exp_mask) <  P*Exp_msk1)
 #else
-				if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
+				if ((word0(&rv) & Exp_mask) <= P*Exp_msk1)
 #endif
 					{
-					if (word0(rv0) == Tiny0
+					if (word0(&rv0) == Tiny0
-					 && word1(rv0) == Tiny1)
+					 && word1(&rv0) == Tiny1)
 						goto undfl;
-					word0(rv) = Tiny0;
+					word0(&rv) = Tiny0;
-					word1(rv) = Tiny1;
+					word1(&rv) = Tiny1;
 					goto cont;
 					}
 				else
-					word0(rv) -= P*Exp_msk1;
+					word0(&rv) -= P*Exp_msk1;
 				}
 			else {
-				adj = aadj1 * ulp(dval(rv));
+				dval(&adj) = dval(&aadj1) * ulp(&rv);
-				dval(rv) += adj;
+				dval(&rv) += adj;
 				}
 #else /*Sudden_Underflow*/
-			/* Compute adj so that the IEEE rounding rules will
+			/* Compute dval(&adj) so that the IEEE rounding rules will
-			 * correctly round rv + adj in some half-way cases.
+			 * correctly round rv + dval(&adj) in some half-way cases.
-			 * If rv * ulp(rv) is denormalized (i.e.,
+			 * If rv * ulp(&rv) is denormalized (i.e.,
 			 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
 			 * trouble from bits lost to denormalization;
 			 * example: 1.2e-307 .
 			 */
 			if (y <= (P-1)*Exp_msk1 && aadj > 1.) {
-				aadj1 = (double)(int)(aadj + 0.5);
+				dval(&aadj1) = (double)(int)(aadj + 0.5);
 				if (!dsign)
-					aadj1 = -aadj1;
+					dval(&aadj1) = -dval(&aadj1);
 				}
-			adj = aadj1 * ulp(dval(rv));
+			dval(&adj) = dval(&aadj1) * ulp(&rv);
-			dval(rv) += adj;
+			dval(&rv) += adj;
 #endif /*Sudden_Underflow*/
 #endif /*Avoid_Underflow*/
 			}
-		z = word0(rv) & Exp_mask;
+		z = word0(&rv) & Exp_mask;
 #ifndef SET_INEXACT
 #ifdef Avoid_Underflow
 		if (!scale)
@@ -928,7 +1013,7 @@ strtod
 			L = (Long)aadj;
 			aadj -= L;
 			/* The tolerances below are conservative. */
-			if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
+			if (dsign || word1(&rv) || word0(&rv) & Bndry_mask) {
 				if (aadj < .4999999 || aadj > .5000001)
 					break;
 				}
@@ -942,12 +1027,17 @@ strtod
 		Bfree(bs);
 		Bfree(delta);
 		}
+	Bfree(bb);
+	Bfree(bd);
+	Bfree(bs);
+	Bfree(bd0);
+	Bfree(delta);
 #ifdef SET_INEXACT
 	if (inexact) {
 		if (!oldinexact) {
-			word0(rv0) = Exp_1 + (70 << Exp_shift);
+			word0(&rv0) = Exp_1 + (70 << Exp_shift);
-			word1(rv0) = 0;
+			word1(&rv0) = 0;
-			dval(rv0) += 1.;
+			dval(&rv0) += 1.;
 			}
 		}
 	else if (!oldinexact)
@@ -955,32 +1045,30 @@ strtod
 #endif
 #ifdef Avoid_Underflow
 	if (scale) {
-		word0(rv0) = Exp_1 - 2*P*Exp_msk1;
+		word0(&rv0) = Exp_1 - 2*P*Exp_msk1;
-		word1(rv0) = 0;
+		word1(&rv0) = 0;
-		dval(rv) *= dval(rv0);
+		dval(&rv) *= dval(&rv0);
 #ifndef NO_ERRNO
 		/* try to avoid the bug of testing an 8087 register value */
-		if (word0(rv) == 0 && word1(rv) == 0)
+#ifdef IEEE_Arith
+		if (!(word0(&rv) & Exp_mask))
+#else
+		if (word0(&rv) == 0 && word1(&rv) == 0)
+#endif
 			errno = ERANGE;
 #endif
 		}
 #endif /* Avoid_Underflow */
 #ifdef SET_INEXACT
-	if (inexact && !(word0(rv) & Exp_mask)) {
+	if (inexact && !(word0(&rv) & Exp_mask)) {
 		/* set underflow bit */
-		dval(rv0) = 1e-300;
+		dval(&rv0) = 1e-300;
-		dval(rv0) *= dval(rv0);
+		dval(&rv0) *= dval(&rv0);
 		}
 #endif
- retfree:
-	Bfree(bb);
-	Bfree(bd);
-	Bfree(bs);
-	Bfree(bd0);
-	Bfree(delta);
  ret:
 	if (se)
 		*se = (char *)s;
-	return sign ? -dval(rv) : dval(rv);
+	return sign ? -dval(&rv) : dval(&rv);
 	}
 

gdtoa/strtodI.c

@@ -33,16 +33,16 @@ THIS SOFTWARE.
 
  static double
 #ifdef KR_headers
-ulpdown(d) double *d;
+ulpdown(d) U *d;
 #else
-ulpdown(double *d)
+ulpdown(U *d)
 #endif
 {
 	double u;
-	ULong *L = (ULong*)d;
+	ULong *L = d->L;
 
-	u = ulp(*d);
+	u = ulp(d);
-	if (!(L[_1] | L[_0] & 0xfffff)
+	if (!(L[_1] | (L[_0] & 0xfffff))
 	 && (L[_0] & 0x7ff00000) > 0x00100000)
 		u *= 0.5;
 	return u;
@@ -55,14 +55,10 @@ strtodI(s, sp, dd) CONST char *s; char **sp; double *dd;
 strtodI(CONST char *s, char **sp, double *dd)
 #endif
 {
-	static CONST 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 };
 	ULong bits[2], sign;
 	Long exp;
 	int j, k;
-	typedef union {
-		double d[2];
-		ULong L[4];
-		} U;
 	U *u;
 
 	k = strtodg(s, sp, &fpi, &exp, bits);
@@ -70,17 +66,17 @@ strtodI(CONST char *s, char **sp, double *dd)
 	sign = k & STRTOG_Neg ? 0x80000000L : 0;
 	switch(k & STRTOG_Retmask) {
 	  case STRTOG_NoNumber:
-		u->d[0] = u->d[1] = 0.;
+		dval(&u[0]) = dval(&u[1]) = 0.;
 		break;
 
 	  case STRTOG_Zero:
-		u->d[0] = u->d[1] = 0.;
+		dval(&u[0]) = dval(&u[1]) = 0.;
 #ifdef Sudden_Underflow
 		if (k & STRTOG_Inexact) {
 			if (sign)
-				u->L[_0] = 0x80100000L;
+				word0(&u[0]) = 0x80100000L;
 			else
-				u->L[2+_0] = 0x100000L;
+				word0(&u[1]) = 0x100000L;
 			}
 		break;
 #else
@@ -88,80 +84,80 @@ strtodI(CONST char *s, char **sp, double *dd)
 #endif
 
 	  case STRTOG_Denormal:
-		u->L[_1] = bits[0];
+		word1(&u[0]) = bits[0];
-		u->L[_0] = bits[1];
+		word0(&u[0]) = bits[1];
 		goto contain;
 
 	  case STRTOG_Normal:
-		u->L[_1] = bits[0];
+		word1(&u[0]) = bits[0];
-		u->L[_0] = (bits[1] & ~0x100000) | ((exp + 0x3ff + 52) << 20);
+		word0(&u[0]) = (bits[1] & ~0x100000) | ((exp + 0x3ff + 52) << 20);
 	  contain:
 		j = k & STRTOG_Inexact;
 		if (sign) {
-			u->L[_0] |= sign;
+			word0(&u[0]) |= sign;
 			j = STRTOG_Inexact - j;
 			}
 		switch(j) {
 		  case STRTOG_Inexlo:
 #ifdef Sudden_Underflow
 			if ((u->L[_0] & 0x7ff00000) < 0x3500000) {
-				u->L[2+_0] = u->L[_0] + 0x3500000;
+				word0(&u[1]) = word0(&u[0]) + 0x3500000;
-				u->L[2+_1] = u->L[_1];
+				word1(&u[1]) = word1(&u[0]);
-				u->d[1] += ulp(u->d[1]);
+				dval(&u[1]) += ulp(&u[1]);
-				u->L[2+_0] -= 0x3500000;
+				word0(&u[1]) -= 0x3500000;
-				if (!(u->L[2+_0] & 0x7ff00000)) {
+				if (!(word0(&u[1]) & 0x7ff00000)) {
-					u->L[2+_0] = sign;
+					word0(&u[1]) = sign;
-					u->L[2+_1] = 0;
+					word1(&u[1]) = 0;
 					}
 				}
 			else
 #endif
-			u->d[1] = u->d[0] + ulp(u->d[0]);
+			dval(&u[1]) = dval(&u[0]) + ulp(&u[0]);
 			break;
 		  case STRTOG_Inexhi:
-			u->d[1] = u->d[0];
+			dval(&u[1]) = dval(&u[0]);
 #ifdef Sudden_Underflow
-			if ((u->L[_0] & 0x7ff00000) < 0x3500000) {
+			if ((word0(&u[0]) & 0x7ff00000) < 0x3500000) {
-				u->L[_0] += 0x3500000;
+				word0(&u[0]) += 0x3500000;
-				u->d[0] -= ulpdown(u->d);
+				dval(&u[0]) -= ulpdown(u);
-				u->L[_0] -= 0x3500000;
+				word0(&u[0]) -= 0x3500000;
-				if (!(u->L[_0] & 0x7ff00000)) {
+				if (!(word0(&u[0]) & 0x7ff00000)) {
-					u->L[_0] = sign;
+					word0(&u[0]) = sign;
-					u->L[_1] = 0;
+					word1(&u[0]) = 0;
 					}
 				}
 			else
 #endif
-			u->d[0] -= ulpdown(u->d);
+			dval(&u[0]) -= ulpdown(u);
 			break;
 		  default:
-			u->d[1] = u->d[0];
+			dval(&u[1]) = dval(&u[0]);
 		  }
 		break;
 
 	  case STRTOG_Infinite:
-		u->L[_0] = u->L[2+_0] = sign | 0x7ff00000;
+		word0(&u[0]) = word0(&u[1]) = sign | 0x7ff00000;
-		u->L[_1] = u->L[2+_1] = 0;
+		word1(&u[0]) = word1(&u[1]) = 0;
 		if (k & STRTOG_Inexact) {
 			if (sign) {
-				u->L[2+_0] = 0xffefffffL;
+				word0(&u[1]) = 0xffefffffL;
-				u->L[2+_1] = 0xffffffffL;
+				word1(&u[1]) = 0xffffffffL;
 				}
 			else {
-				u->L[_0] = 0x7fefffffL;
+				word0(&u[0]) = 0x7fefffffL;
-				u->L[_1] = 0xffffffffL;
+				word1(&u[0]) = 0xffffffffL;
 				}
 			}
 		break;
 
 	  case STRTOG_NaN:
-		u->L[0] = u->L[2] = d_QNAN0;
+		u->L[0] = (u+1)->L[0] = d_QNAN0;
-		u->L[1] = u->L[3] = d_QNAN1;
+		u->L[1] = (u+1)->L[1] = d_QNAN1;
 		break;
 
 	  case STRTOG_NaNbits:
-		u->L[_0] = u->L[2+_0] = 0x7ff00000 | sign | bits[1];
+		word0(&u[0]) = word0(&u[1]) = 0x7ff00000 | sign | bits[1];
-		u->L[_1] = u->L[2+_1] = bits[0];
+		word1(&u[0]) = word1(&u[1]) = bits[0];
 	  }
 	return k;
 	}

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
@@ -173,9 +172,9 @@ set_ones(Bigint *b, int n)
 rvOK
 #ifdef KR_headers
  (d, fpi, exp, bits, exact, rd, irv)
- double d; FPI *fpi; Long *exp; ULong *bits; int exact, rd, *irv;
+ U *d; FPI *fpi; Long *exp; ULong *bits; int exact, rd, *irv;
 #else
- (double d, CONST FPI *fpi, Long *exp, ULong *bits, int exact, int rd, int *irv)
+ (U *d, FPI *fpi, Long *exp, ULong *bits, int exact, int rd, int *irv)
 #endif
 {
 	Bigint *b;
@@ -183,7 +182,7 @@ rvOK
 	int bdif, e, j, k, k1, nb, rv;
 
 	carry = rv = 0;
-	b = d2b(d, &e, &bdif);
+	b = d2b(dval(d), &e, &bdif);
 	bdif -= nb = fpi->nbits;
 	e += bdif;
 	if (bdif <= 0) {
@@ -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;
@@ -292,15 +291,12 @@ rvOK
 
  static int
 #ifdef KR_headers
-mantbits(d) double d;
+mantbits(d) U *d;
 #else
-mantbits(double _d)
+mantbits(U *d)
 #endif
 {
 	ULong L;
-	U d;
-
-	dval(d) = _d;
 #ifdef VAX
 	L = word1(d) << 16 | word1(d) >> 16;
 	if (L)
@@ -322,7 +318,7 @@ strtodg
 	(s00, se, fpi, exp, bits)
 	CONST char *s00; char **se; FPI *fpi; Long *exp; ULong *bits;
 #else
-	(CONST char *s00, char **se, CONST FPI *fpi, Long *exp, ULong *bits)
+	(CONST char *s00, char **se, FPI *fpi, Long *exp, ULong *bits)
 #endif
 {
 	int abe, abits, asub;
@@ -332,14 +328,35 @@ strtodg
 	int sudden_underflow;
 	CONST char *s, *s0, *s1;
 	double adj0, tol;
-	U adj, rv;
 	Long L;
-	ULong y, z;
+	U adj, rv;
+	ULong *b, *be, y, z;
 	Bigint *ab, *bb, *bb1, *bd, *bd0, *bs, *delta, *rvb, *rvb0;
+#ifdef USE_LOCALE /*{{*/
+#ifdef NO_LOCALE_CACHE
+	char *decimalpoint = localeconv()->decimal_point;
+	int dplen = strlen(decimalpoint);
+#else
+	char *decimalpoint;
+	static char *decimalpoint_cache;
+	static int dplen;
+	if (!(s0 = decimalpoint_cache)) {
+		s0 = localeconv()->decimal_point;
+		if ((decimalpoint_cache = (char*)MALLOC(strlen(s0) + 1))) {
+			strcpy(decimalpoint_cache, s0);
+			s0 = decimalpoint_cache;
+			}
+		dplen = strlen(s0);
+		}
+	decimalpoint = (char*)s0;
+#endif /*NO_LOCALE_CACHE*/
+#else  /*USE_LOCALE}{*/
+#define dplen 1
+#endif /*USE_LOCALE}}*/
 
 	irv = STRTOG_Zero;
 	denorm = sign = nz0 = nz = 0;
-	dval(rv) = 0.;
+	dval(&rv) = 0.;
 	rvb = 0;
 	nbits = fpi->nbits;
 	for(s = s00;;s++) switch(*s) {
@@ -394,13 +411,17 @@ strtodg
 			z = 10*z + c - '0';
 	nd0 = nd;
 #ifdef USE_LOCALE
-	if (c == *localeconv()->decimal_point)
+	if (c == *decimalpoint) {
+		for(i = 1; decimalpoint[i]; ++i)
+			if (s[i] != decimalpoint[i])
+				goto dig_done;
+		s += i;
+		c = *s;
 #else
-	if (c == '.')
+	if (c == '.') {
-#endif
-		{
-		decpt = 1;
 		c = *++s;
+#endif
+		decpt = 1;
 		if (!nd) {
 			for(; c == '0'; c = *++s)
 				nz++;
@@ -429,7 +450,7 @@ strtodg
 				nz = 0;
 				}
 			}
-		}
+		}/*}*/
  dig_done:
 	e = 0;
 	if (c == 'e' || c == 'E') {
@@ -527,13 +548,13 @@ strtodg
 	if (!nd0)
 		nd0 = nd;
 	k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
-	dval(rv) = y;
+	dval(&rv) = y;
 	if (k > 9)
-		dval(rv) = tens[k - 9] * dval(rv) + z;
+		dval(&rv) = tens[k - 9] * dval(&rv) + z;
 	bd0 = 0;
 	if (nbits <= P && nd <= DBL_DIG) {
 		if (!e) {
-			if (rvOK(dval(rv), fpi, exp, bits, 1, rd, &irv))
+			if (rvOK(&rv, fpi, exp, bits, 1, rd, &irv))
 				goto ret;
 			}
 		else if (e > 0) {
@@ -541,9 +562,9 @@ strtodg
 #ifdef VAX
 				goto vax_ovfl_check;
 #else
-				i = fivesbits[e] + mantbits(dval(rv)) <= P;
+				i = fivesbits[e] + mantbits(&rv) <= P;
-				/* rv = */ rounded_product(dval(rv), tens[e]);
+				/* rv = */ rounded_product(dval(&rv), tens[e]);
-				if (rvOK(dval(rv), fpi, exp, bits, i, rd, &irv))
+				if (rvOK(&rv, fpi, exp, bits, i, rd, &irv))
 					goto ret;
 				e1 -= e;
 				goto rv_notOK;
@@ -556,32 +577,32 @@ strtodg
 				 */
 				e2 = e - i;
 				e1 -= i;
-				dval(rv) *= tens[i];
+				dval(&rv) *= tens[i];
 #ifdef VAX
 				/* VAX exponent range is so narrow we must
 				 * worry about overflow here...
 				 */
  vax_ovfl_check:
-				dval(adj) = dval(rv);
+				dval(&adj) = dval(&rv);
-				word0(adj) -= P*Exp_msk1;
+				word0(&adj) -= P*Exp_msk1;
-				/* adj = */ rounded_product(dval(adj), tens[e2]);
+				/* adj = */ rounded_product(dval(&adj), tens[e2]);
-				if ((word0(adj) & Exp_mask)
+				if ((word0(&adj) & Exp_mask)
 				 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
 					goto rv_notOK;
-				word0(adj) += P*Exp_msk1;
+				word0(&adj) += P*Exp_msk1;
-				dval(rv) = dval(adj);
+				dval(&rv) = dval(&adj);
 #else
-				/* rv = */ rounded_product(dval(rv), tens[e2]);
+				/* rv = */ rounded_product(dval(&rv), tens[e2]);
 #endif
-				if (rvOK(dval(rv), fpi, exp, bits, 0, rd, &irv))
+				if (rvOK(&rv, fpi, exp, bits, 0, rd, &irv))
 					goto ret;
 				e1 -= e2;
 				}
 			}
 #ifndef Inaccurate_Divide
 		else if (e >= -Ten_pmax) {
-			/* rv = */ rounded_quotient(dval(rv), tens[-e]);
+			/* rv = */ rounded_quotient(dval(&rv), tens[-e]);
-			if (rvOK(dval(rv), fpi, exp, bits, 0, rd, &irv))
+			if (rvOK(&rv, fpi, exp, bits, 0, rd, &irv))
 				goto ret;
 			e1 -= e;
 			}
@@ -595,45 +616,45 @@ strtodg
 	e2 = 0;
 	if (e1 > 0) {
 		if ( (i = e1 & 15) !=0)
-			dval(rv) *= tens[i];
+			dval(&rv) *= tens[i];
 		if (e1 &= ~15) {
 			e1 >>= 4;
-			while(e1 >= (1 << n_bigtens-1)) {
+			while(e1 >= (1 << (n_bigtens-1))) {
-				e2 += ((word0(rv) & Exp_mask)
+				e2 += ((word0(&rv) & Exp_mask)
 					>> Exp_shift1) - Bias;
-				word0(rv) &= ~Exp_mask;
+				word0(&rv) &= ~Exp_mask;
-				word0(rv) |= Bias << Exp_shift1;
+				word0(&rv) |= Bias << Exp_shift1;
-				dval(rv) *= bigtens[n_bigtens-1];
+				dval(&rv) *= bigtens[n_bigtens-1];
-				e1 -= 1 << n_bigtens-1;
+				e1 -= 1 << (n_bigtens-1);
 				}
-			e2 += ((word0(rv) & Exp_mask) >> Exp_shift1) - Bias;
+			e2 += ((word0(&rv) & Exp_mask) >> Exp_shift1) - Bias;
-			word0(rv) &= ~Exp_mask;
+			word0(&rv) &= ~Exp_mask;
-			word0(rv) |= Bias << Exp_shift1;
+			word0(&rv) |= Bias << Exp_shift1;
 			for(j = 0; e1 > 0; j++, e1 >>= 1)
 				if (e1 & 1)
-					dval(rv) *= bigtens[j];
+					dval(&rv) *= bigtens[j];
 			}
 		}
 	else if (e1 < 0) {
 		e1 = -e1;
 		if ( (i = e1 & 15) !=0)
-			dval(rv) /= tens[i];
+			dval(&rv) /= tens[i];
 		if (e1 &= ~15) {
 			e1 >>= 4;
-			while(e1 >= (1 << n_bigtens-1)) {
+			while(e1 >= (1 << (n_bigtens-1))) {
-				e2 += ((word0(rv) & Exp_mask)
+				e2 += ((word0(&rv) & Exp_mask)
 					>> Exp_shift1) - Bias;
-				word0(rv) &= ~Exp_mask;
+				word0(&rv) &= ~Exp_mask;
-				word0(rv) |= Bias << Exp_shift1;
+				word0(&rv) |= Bias << Exp_shift1;
-				dval(rv) *= tinytens[n_bigtens-1];
+				dval(&rv) *= tinytens[n_bigtens-1];
-				e1 -= 1 << n_bigtens-1;
+				e1 -= 1 << (n_bigtens-1);
 				}
-			e2 += ((word0(rv) & Exp_mask) >> Exp_shift1) - Bias;
+			e2 += ((word0(&rv) & Exp_mask) >> Exp_shift1) - Bias;
-			word0(rv) &= ~Exp_mask;
+			word0(&rv) &= ~Exp_mask;
-			word0(rv) |= Bias << Exp_shift1;
+			word0(&rv) |= Bias << Exp_shift1;
 			for(j = 0; e1 > 0; j++, e1 >>= 1)
 				if (e1 & 1)
-					dval(rv) *= tinytens[j];
+					dval(&rv) *= tinytens[j];
 			}
 		}
 #ifdef IBM
@@ -644,7 +665,7 @@ strtodg
 	 */
 	e2 <<= 2;
 #endif
-	rvb = d2b(dval(rv), &rve, &rvbits);	/* rv = rvb * 2^rve */
+	rvb = d2b(dval(&rv), &rve, &rvbits);	/* rv = rvb * 2^rve */
 	rve += e2;
 	if ((j = rvbits - nbits) > 0) {
 		rshift(rvb, j);
@@ -688,7 +709,7 @@ strtodg
 
 	/* Put digits into bd: true value = bd * 10^e */
 
-	bd0 = s2b(s0, nd0, nd, y);
+	bd0 = s2b(s0, nd0, nd, y, dplen);
 
 	for(;;) {
 		bd = Balloc(bd0->k);
@@ -822,7 +843,7 @@ strtodg
 				}
 			else
 				irv = STRTOG_Normal | STRTOG_Inexhi;
-			if (bbbits < nbits && !denorm || !(rvb->x[0] & 1))
+			if ((bbbits < nbits && !denorm) || !(rvb->x[0] & 1))
 				break;
 			if (dsign) {
 				rvb = increment(rvb);
@@ -839,7 +860,7 @@ strtodg
 				}
 			break;
 			}
-		if ((dval(adj) = ratio(delta, bs)) <= 2.) {
+		if ((dval(&adj) = ratio(delta, bs)) <= 2.) {
  adj1:
 			inex = STRTOG_Inexlo;
 			if (dsign) {
@@ -853,15 +874,15 @@ strtodg
 				irv = STRTOG_Underflow | STRTOG_Inexlo;
 				break;
 				}
-			adj0 = dval(adj) = 1.;
+			adj0 = dval(&adj) = 1.;
 			}
 		else {
-			adj0 = dval(adj) *= 0.5;
+			adj0 = dval(&adj) *= 0.5;
 			if (dsign) {
 				asub = 0;
 				inex = STRTOG_Inexlo;
 				}
-			if (dval(adj) < 2147483647.) {
+			if (dval(&adj) < 2147483647.) {
 				L = (Long)adj0;
 				adj0 -= L;
 				switch(rd) {
@@ -880,12 +901,12 @@ strtodg
 						inex = STRTOG_Inexact - inex;
 						}
 				  }
-				dval(adj) = L;
+				dval(&adj) = L;
 				}
 			}
 		y = rve + rvbits;
 
-		/* adj *= ulp(dval(rv)); */
+		/* adj *= ulp(dval(&rv)); */
 		/* if (asub) rv -= adj; else rv += adj; */
 
 		if (!denorm && rvbits < nbits) {
@@ -893,7 +914,7 @@ strtodg
 			rve -= j;
 			rvbits = nbits;
 			}
-		ab = d2b(dval(adj), &abe, &abits);
+		ab = d2b(dval(&adj), &abe, &abits);
 		if (abe < 0)
 			rshift(ab, -abe);
 		else if (abe > 0)
@@ -947,15 +968,15 @@ strtodg
 		z = rve + rvbits;
 		if (y == z && L) {
 			/* Can we stop now? */
-			tol = dval(adj) * 5e-16; /* > max rel error */
+			tol = dval(&adj) * 5e-16; /* > max rel error */
-			dval(adj) = adj0 - .5;
+			dval(&adj) = adj0 - .5;
-			if (dval(adj) < -tol) {
+			if (dval(&adj) < -tol) {
 				if (adj0 > tol) {
 					irv |= inex;
 					break;
 					}
 				}
-			else if (dval(adj) > tol && adj0 < 1. - tol) {
+			else if (dval(&adj) > tol && adj0 < 1. - tol) {
 				irv |= inex;
 				break;
 				}
@@ -980,6 +1001,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 + 31) >> 5);
+		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;
@@ -994,12 +1038,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)

gdtoa/strtodnrp.c

@@ -44,7 +44,7 @@ strtod(s, sp) CONST char *s; char **sp;
 strtod(CONST char *s, char **sp)
 #endif
 {
-	static CONST 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 };
 	ULong bits[2];
 	Long exp;
 	int k;

gdtoa/strtof.c

@@ -38,13 +38,18 @@ strtof(s, sp) CONST char *s; char **sp;
 strtof(CONST char *s, char **sp)
 #endif
 {
-	static CONST FPI fpi = { 24, 1-127-24+1,  254-127-24+1, 1, SI };
+	static FPI fpi0 = { 24, 1-127-24+1,  254-127-24+1, 1, SI };
 	ULong bits[1];
 	Long exp;
 	int k;
 	union { ULong L[1]; float f; } u;
+#ifdef Honor_FLT_ROUNDS
+#include "gdtoa_fltrnds.h"
+#else
+#define fpi &fpi0
+#endif
 
-	k = strtodg(s, sp, &fpi, &exp, bits);
+	k = strtodg(s, sp, fpi, &exp, bits);
 	switch(k & STRTOG_Retmask) {
 	  case STRTOG_NoNumber:
 	  case STRTOG_Zero:
@@ -53,7 +58,7 @@ strtof(CONST char *s, char **sp)
 
 	  case STRTOG_Normal:
 	  case STRTOG_NaNbits:
-		u.L[0] = bits[0] & 0x7fffff | exp + 0x7f + 23 << 23;
+		u.L[0] = (bits[0] & 0x7fffff) | ((exp + 0x7f + 23) << 23);
 		break;
 
 	  case STRTOG_Denormal:

gdtoa/strtopQ.c

@@ -49,6 +49,9 @@ THIS SOFTWARE.
 #define _3 0
 #endif
 
+ extern ULong NanDflt_Q_D2A[4];
+
+
  int
 #ifdef KR_headers
 strtopQ(s, sp, V) CONST char *s; char **sp; void *V;
@@ -56,13 +59,18 @@ strtopQ(s, sp, V) CONST char *s; char **sp; void *V;
 strtopQ(CONST char *s, char **sp, void *V)
 #endif
 {
-	static CONST FPI fpi = { 113, 1-16383-113+1, 32766 - 16383 - 113 + 1, 1, SI };
+	static FPI fpi0 = { 113, 1-16383-113+1, 32766 - 16383 - 113 + 1, 1, SI };
 	ULong bits[4];
 	Long exp;
 	int k;
 	ULong *L = (ULong*)V;
+#ifdef Honor_FLT_ROUNDS
+#include "gdtoa_fltrnds.h"
+#else
+#define fpi &fpi0
+#endif
 
-	k = strtodg(s, sp, &fpi, &exp, bits);
+	k = strtodg(s, sp, fpi, &exp, bits);
 	switch(k & STRTOG_Retmask) {
 	  case STRTOG_NoNumber:
 	  case STRTOG_Zero:
@@ -90,10 +98,10 @@ strtopQ(CONST char *s, char **sp, void *V)
 		break;
 
 	  case STRTOG_NaN:
-		L[0] = ld_QNAN0;
+		L[_0] = NanDflt_Q_D2A[3];
-		L[1] = ld_QNAN1;
+		L[_1] = NanDflt_Q_D2A[2];
-		L[2] = ld_QNAN2;
+		L[_2] = NanDflt_Q_D2A[1];
-		L[3] = ld_QNAN3;
+		L[_3] = NanDflt_Q_D2A[0];
 	  }
 	if (k & STRTOG_Neg)
 		L[_0] |= 0x80000000L;

gdtoa/strtopd.c

@@ -38,12 +38,17 @@ strtopd(s, sp, d) char *s; char **sp; double *d;
 strtopd(CONST char *s, char **sp, double *d)
 #endif
 {
-	static CONST FPI fpi0 = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI };
+	static FPI fpi0 = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI };
 	ULong bits[2];
 	Long exp;
 	int k;
+#ifdef Honor_FLT_ROUNDS
+#include "gdtoa_fltrnds.h"
+#else
+#define fpi &fpi0
+#endif
 
-	k = strtodg(s, sp, &fpi0, &exp, bits);
+	k = strtodg(s, sp, fpi, &exp, bits);
 	ULtod((ULong*)d, bits, exp, k);
 	return k;
 	}

gdtoa/strtopdd.c

@@ -39,9 +39,9 @@ strtopdd(CONST char *s, char **sp, double *dd)
 #endif
 {
 #ifdef Sudden_Underflow
-	static CONST FPI fpi = { 106, 1-1023, 2046-1023-106+1, 1, 1 };
+	static FPI fpi0 = { 106, 1-1023, 2046-1023-106+1, 1, 1 };
 #else
-	static CONST FPI fpi = { 106, 1-1023-53+1, 2046-1023-106+1, 1, 0 };
+	static FPI fpi0 = { 106, 1-1023-53+1, 2046-1023-106+1, 1, 0 };
 #endif
 	ULong bits[4];
 	Long exp;
@@ -51,8 +51,13 @@ strtopdd(CONST char *s, char **sp, double *dd)
 		ULong L[4];
 		} U;
 	U *u;
+#ifdef Honor_FLT_ROUNDS
+#include "gdtoa_fltrnds.h"
+#else
+#define fpi &fpi0
+#endif
 
-	rv = strtodg(s, sp, &fpi, &exp, bits);
+	rv = strtodg(s, sp, fpi, &exp, bits);
 	u = (U*)dd;
 	switch(rv & STRTOG_Retmask) {
 	  case STRTOG_NoNumber:
@@ -62,8 +67,8 @@ strtopdd(CONST char *s, char **sp, double *dd)
 
 	  case STRTOG_Normal:
 		u->L[_1] = (bits[1] >> 21 | bits[2] << 11) & 0xffffffffL;
-		u->L[_0] = bits[2] >> 21 | bits[3] << 11 & 0xfffff
+		u->L[_0] = (bits[2] >> 21) | ((bits[3] << 11) & 0xfffff)
-			  | exp + 0x3ff + 105 << 20;
+			  | ((exp + 0x3ff + 105) << 20);
 		exp += 0x3ff + 52;
 		if (bits[1] &= 0x1fffff) {
 			i = hi0bits(bits[1]) - 11;
@@ -74,7 +79,7 @@ strtopdd(CONST char *s, char **sp, double *dd)
 			else
 				exp -= i;
 			if (i > 0) {
-				bits[1] = bits[1] << i | bits[0] >> 32-i;
+				bits[1] = bits[1] << i | bits[0] >> (32-i);
 				bits[0] = bits[0] << i & 0xffffffffL;
 				}
 			}
@@ -87,11 +92,11 @@ strtopdd(CONST char *s, char **sp, double *dd)
 			else
 				exp -= i;
 			if (i < 32) {
-				bits[1] = bits[0] >> 32 - i;
+				bits[1] = bits[0] >> (32 - i);
 				bits[0] = bits[0] << i & 0xffffffffL;
 				}
 			else {
-				bits[1] = bits[0] << i - 32;
+				bits[1] = bits[0] << (i - 32);
 				bits[0] = 0;
 				}
 			}
@@ -100,7 +105,7 @@ strtopdd(CONST char *s, char **sp, double *dd)
 			break;
 			}
 		u->L[2+_1] = bits[0];
-		u->L[2+_0] = bits[1] & 0xfffff | exp << 20;
+		u->L[2+_0] = (bits[1] & 0xfffff) | (exp << 20);
 		break;
 
 	  case STRTOG_Denormal:
@@ -119,10 +124,10 @@ strtopdd(CONST char *s, char **sp, double *dd)
 	  nearly_normal:
 		i = hi0bits(bits[3]) - 11;	/* i >= 12 */
 		j = 32 - i;
-		u->L[_0] = (bits[3] << i | bits[2] >> j) & 0xfffff
+		u->L[_0] = ((bits[3] << i | bits[2] >> j) & 0xfffff)
-			| 65 - i << 20;
+			| ((65 - i) << 20);
 		u->L[_1] = (bits[2] << i | bits[1] >> j) & 0xffffffffL;
-		u->L[2+_0] = bits[1] & (1L << j) - 1;
+		u->L[2+_0] = bits[1] & ((1L << j) - 1);
 		u->L[2+_1] = bits[0];
 		break;
 
@@ -131,34 +136,34 @@ strtopdd(CONST char *s, char **sp, double *dd)
 		if (i < 0) {
 			j = -i;
 			i += 32;
-			u->L[_0] = bits[2] >> j & 0xfffff | (33 + j) << 20;
+			u->L[_0] = (bits[2] >> j & 0xfffff) | (33 + j) << 20;
-			u->L[_1] = (bits[2] << i | bits[1] >> j) & 0xffffffffL;
+			u->L[_1] = ((bits[2] << i) | (bits[1] >> j)) & 0xffffffffL;
-			u->L[2+_0] = bits[1] & (1L << j) - 1;
+			u->L[2+_0] = bits[1] & ((1L << j) - 1);
 			u->L[2+_1] = bits[0];
 			break;
 			}
 		if (i == 0) {
-			u->L[_0] = bits[2] & 0xfffff | 33 << 20;
+			u->L[_0] = (bits[2] & 0xfffff) | (33 << 20);
 			u->L[_1] = bits[1];
 			u->L[2+_0] = 0;
 			u->L[2+_1] = bits[0];
 			break;
 			}
 		j = 32 - i;
-		u->L[_0] = (bits[2] << i | bits[1] >> j) & 0xfffff
+		u->L[_0] = (((bits[2] << i) | (bits[1] >> j)) & 0xfffff)
-				| j + 1 << 20;
+				| ((j + 1) << 20);
 		u->L[_1] = (bits[1] << i | bits[0] >> j) & 0xffffffffL;
 		u->L[2+_0] = 0;
-		u->L[2+_1] = bits[0] & (1L << j) - 1;
+		u->L[2+_1] = bits[0] & ((1L << j) - 1);
 		break;
 
 	  hardly_normal:
 		j = 11 - hi0bits(bits[1]);
 		i = 32 - j;
-		u->L[_0] = bits[1] >> j & 0xfffff | j + 1 << 20;
+		u->L[_0] = (bits[1] >> j & 0xfffff) | ((j + 1) << 20);
 		u->L[_1] = (bits[1] << i | bits[0] >> j) & 0xffffffffL;
 		u->L[2+_0] = 0;
-		u->L[2+_1] = bits[0] & (1L << j) - 1;
+		u->L[2+_1] = bits[0] & ((1L << j) - 1);
 		break;
 
 	  case STRTOG_Infinite:

gdtoa/strtopf.c

@@ -38,12 +38,17 @@ strtopf(s, sp, f) CONST char *s; char **sp; float *f;
 strtopf(CONST char *s, char **sp, float *f)
 #endif
 {
-	static CONST FPI fpi = { 24, 1-127-24+1,  254-127-24+1, 1, SI };
+	static FPI fpi0 = { 24, 1-127-24+1,  254-127-24+1, 1, SI };
 	ULong bits[1], *L;
 	Long exp;
 	int k;
+#ifdef Honor_FLT_ROUNDS
+#include "gdtoa_fltrnds.h"
+#else
+#define fpi &fpi0
+#endif
 
-	k = strtodg(s, sp, &fpi, &exp, bits);
+	k = strtodg(s, sp, fpi, &exp, bits);
 	L = (ULong*)f;
 	switch(k & STRTOG_Retmask) {
 	  case STRTOG_NoNumber:
@@ -53,7 +58,7 @@ strtopf(CONST char *s, char **sp, float *f)
 
 	  case STRTOG_Normal:
 	  case STRTOG_NaNbits:
-		L[0] = bits[0] & 0x7fffff | exp + 0x7f + 23 << 23;
+		L[0] = (bits[0] & 0x7fffff) | ((exp + 0x7f + 23) << 23);
 		break;
 
 	  case STRTOG_Denormal:

gdtoa/strtopx.c

@@ -31,6 +31,8 @@ THIS SOFTWARE.
 
 #include "gdtoaimp.h"
 
+ extern UShort NanDflt_ldus_D2A[5];
+
 #undef _0
 #undef _1
 
@@ -58,13 +60,18 @@ strtopx(s, sp, V) CONST char *s; char **sp; void *V;
 strtopx(CONST char *s, char **sp, void *V)
 #endif
 {
-	static CONST FPI fpi = { 64, 1-16383-64+1, 32766 - 16383 - 64 + 1, 1, SI };
+	static FPI fpi0 = { 64, 1-16383-64+1, 32766 - 16383 - 64 + 1, 1, SI };
 	ULong bits[2];
 	Long exp;
 	int k;
 	UShort *L = (UShort*)V;
+#ifdef Honor_FLT_ROUNDS
+#include "gdtoa_fltrnds.h"
+#else
+#define fpi &fpi0
+#endif
 
-	k = strtodg(s, sp, &fpi, &exp, bits);
+	k = strtodg(s, sp, fpi, &exp, bits);
 	switch(k & STRTOG_Retmask) {
 	  case STRTOG_NoNumber:
 	  case STRTOG_Zero:
@@ -87,15 +94,16 @@ strtopx(CONST char *s, char **sp, void *V)
 
 	  case STRTOG_Infinite:
 		L[_0] = 0x7fff;
-		L[_1] = L[_2] = L[_3] = L[_4] = 0;
+		L[_1] = 0x8000;
+		L[_2] = L[_3] = L[_4] = 0;
 		break;
 
 	  case STRTOG_NaN:
-		L[0] = ldus_QNAN0;
+		L[_4] = NanDflt_ldus_D2A[0];
-		L[1] = ldus_QNAN1;
+		L[_3] = NanDflt_ldus_D2A[1];
-		L[2] = ldus_QNAN2;
+		L[_2] = NanDflt_ldus_D2A[2];
-		L[3] = ldus_QNAN3;
+		L[_1] = NanDflt_ldus_D2A[3];
-		L[4] = ldus_QNAN4;
+		L[_0] = NanDflt_ldus_D2A[4];
 	  }
 	if (k & STRTOG_Neg)
 		L[_0] |= 0x8000;

gdtoa/strtopxL.c

@@ -31,6 +31,8 @@ THIS SOFTWARE.
 
 #include "gdtoaimp.h"
 
+ extern ULong NanDflt_xL_D2A[3];
+
 #undef _0
 #undef _1
 
@@ -54,13 +56,18 @@ strtopxL(s, sp, V) CONST char *s; char **sp; void *V;
 strtopxL(CONST char *s, char **sp, void *V)
 #endif
 {
-	static CONST FPI fpi = { 64, 1-16383-64+1, 32766 - 16383 - 64 + 1, 1, SI };
+	static FPI fpi0 = { 64, 1-16383-64+1, 32766 - 16383 - 64 + 1, 1, SI };
 	ULong bits[2];
 	Long exp;
 	int k;
 	ULong *L = (ULong*)V;
+#ifdef Honor_FLT_ROUNDS
+#include "gdtoa_fltrnds.h"
+#else
+#define fpi &fpi0
+#endif
 
-	k = strtodg(s, sp, &fpi, &exp, bits);
+	k = strtodg(s, sp, fpi, &exp, bits);
 	switch(k & STRTOG_Retmask) {
 	  case STRTOG_NoNumber:
 	  case STRTOG_Zero:
@@ -77,13 +84,14 @@ strtopxL(CONST char *s, char **sp, void *V)
 
 	  case STRTOG_Infinite:
 		L[_0] = 0x7fff << 16;
-		L[_1] = L[_2] = 0;
+		L[_1] = 0x80000000;
+		L[_2] = 0;
 		break;
 
 	  case STRTOG_NaN:
-		L[0] = ld_QNAN0;
+		L[_0] = NanDflt_xL_D2A[2];
-		L[1] = ld_QNAN1;
+		L[_1] = NanDflt_xL_D2A[1];
-		L[2] = ld_QNAN2;
+		L[_2] = NanDflt_xL_D2A[0];
 	  }
 	if (k & STRTOG_Neg)
 		L[_0] |= 0x80000000L;

gdtoa/strtorQ.c

@@ -49,6 +49,8 @@ THIS SOFTWARE.
 #define _3 0
 #endif
 
+ extern ULong NanDflt_Q_D2A[4];
+
  void
 #ifdef KR_headers
 ULtoQ(L, bits, exp, k) ULong *L; ULong *bits; Long exp; int k;
@@ -83,10 +85,10 @@ ULtoQ(ULong *L, ULong *bits, Long exp, int k)
 		break;
 
 	  case STRTOG_NaN:
-		L[0] = ld_QNAN0;
+		L[_0] = NanDflt_Q_D2A[3];
-		L[1] = ld_QNAN1;
+		L[_1] = NanDflt_Q_D2A[2];
-		L[2] = ld_QNAN2;
+		L[_2] = NanDflt_Q_D2A[1];
-		L[3] = ld_QNAN3;
+		L[_3] = NanDflt_Q_D2A[0];
 	  }
 	if (k & STRTOG_Neg)
 		L[_0] |= 0x80000000L;
@@ -99,9 +101,8 @@ strtorQ(s, sp, rounding, L) CONST char *s; char **sp; int rounding; void *L;
 strtorQ(CONST char *s, char **sp, int rounding, void *L)
 #endif
 {
-	static CONST FPI fpi0 = { 113, 1-16383-113+1, 32766-16383-113+1, 1, SI };
+	static FPI fpi0 = { 113, 1-16383-113+1, 32766-16383-113+1, 1, SI };
-	CONST FPI *fpi;
+	FPI *fpi, fpi1;
-	FPI fpi1;
 	ULong bits[4];
 	Long exp;
 	int k;

gdtoa/strtord.c

@@ -31,6 +31,8 @@ THIS SOFTWARE.
 
 #include "gdtoaimp.h"
 
+ extern ULong NanDflt_d_D2A[2];
+
  void
 #ifdef KR_headers
 ULtod(L, bits, exp, k) ULong *L; ULong *bits; Long exp; int k;
@@ -61,8 +63,8 @@ ULtod(ULong *L, ULong *bits, Long exp, int k)
 		break;
 
 	  case STRTOG_NaN:
-		L[0] = d_QNAN0;
+		L[_0] = NanDflt_d_D2A[1];
-		L[1] = d_QNAN1;
+		L[_1] = NanDflt_d_D2A[0];
 	  }
 	if (k & STRTOG_Neg)
 		L[_0] |= 0x80000000L;
@@ -75,9 +77,8 @@ strtord(s, sp, rounding, d) CONST char *s; char **sp; int rounding; double *d;
 strtord(CONST char *s, char **sp, int rounding, double *d)
 #endif
 {
-	static CONST FPI fpi0 = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI };
+	static FPI fpi0 = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI };
-	CONST FPI *fpi;
+	FPI *fpi, fpi1;
-	FPI fpi1;
 	ULong bits[2];
 	Long exp;
 	int k;

gdtoa/strtordd.c

@@ -31,6 +31,8 @@ THIS SOFTWARE.
 
 #include "gdtoaimp.h"
 
+ extern ULong NanDflt_d_D2A[2];
+
  void
 #ifdef KR_headers
 ULtodd(L, bits, exp, k) ULong *L; ULong *bits; Long exp; int k;
@@ -48,8 +50,8 @@ ULtodd(ULong *L, ULong *bits, Long exp, int k)
 
 	  case STRTOG_Normal:
 		L[_1] = (bits[1] >> 21 | bits[2] << 11) & (ULong)0xffffffffL;
-		L[_0] = bits[2] >> 21 | bits[3] << 11 & 0xfffff
+		L[_0] = (bits[2] >> 21) | (bits[3] << 11 & 0xfffff)
-			  | exp + 0x3ff + 105 << 20;
+			  | ((exp + 0x3ff + 105) << 20);
 		exp += 0x3ff + 52;
 		if (bits[1] &= 0x1fffff) {
 			i = hi0bits(bits[1]) - 11;
@@ -60,7 +62,7 @@ ULtodd(ULong *L, ULong *bits, Long exp, int k)
 			else
 				exp -= i;
 			if (i > 0) {
-				bits[1] = bits[1] << i | bits[0] >> 32-i;
+				bits[1] = bits[1] << i | bits[0] >> (32-i);
 				bits[0] = bits[0] << i & (ULong)0xffffffffL;
 				}
 			}
@@ -73,11 +75,11 @@ ULtodd(ULong *L, ULong *bits, Long exp, int k)
 			else
 				exp -= i;
 			if (i < 32) {
-				bits[1] = bits[0] >> 32 - i;
+				bits[1] = bits[0] >> (32 - i);
 				bits[0] = bits[0] << i & (ULong)0xffffffffL;
 				}
 			else {
-				bits[1] = bits[0] << i - 32;
+				bits[1] = bits[0] << (i - 32);
 				bits[0] = 0;
 				}
 			}
@@ -86,7 +88,7 @@ ULtodd(ULong *L, ULong *bits, Long exp, int k)
 			break;
 			}
 		L[2+_1] = bits[0];
-		L[2+_0] = bits[1] & 0xfffff | exp << 20;
+		L[2+_0] = (bits[1] & 0xfffff) | (exp << 20);
 		break;
 
 	  case STRTOG_Denormal:
@@ -105,10 +107,10 @@ ULtodd(ULong *L, ULong *bits, Long exp, int k)
 	  nearly_normal:
 		i = hi0bits(bits[3]) - 11;	/* i >= 12 */
 		j = 32 - i;
-		L[_0] = (bits[3] << i | bits[2] >> j) & 0xfffff
+		L[_0] = ((bits[3] << i | bits[2] >> j) & 0xfffff)
-			| 65 - i << 20;
+			| ((65 - i) << 20);
 		L[_1] = (bits[2] << i | bits[1] >> j) & 0xffffffffL;
-		L[2+_0] = bits[1] & ((ULong)1L << j) - 1;
+		L[2+_0] = bits[1] & (((ULong)1L << j) - 1);
 		L[2+_1] = bits[0];
 		break;
 
@@ -117,34 +119,34 @@ ULtodd(ULong *L, ULong *bits, Long exp, int k)
 		if (i < 0) {
 			j = -i;
 			i += 32;
-			L[_0] = bits[2] >> j & 0xfffff | (33 + j) << 20;
+			L[_0] = (bits[2] >> j & 0xfffff) | ((33 + j) << 20);
 			L[_1] = (bits[2] << i | bits[1] >> j) & 0xffffffffL;
-			L[2+_0] = bits[1] & ((ULong)1L << j) - 1;
+			L[2+_0] = bits[1] & (((ULong)1L << j) - 1);
 			L[2+_1] = bits[0];
 			break;
 			}
 		if (i == 0) {
-			L[_0] = bits[2] & 0xfffff | 33 << 20;
+			L[_0] = (bits[2] & 0xfffff) | (33 << 20);
 			L[_1] = bits[1];
 			L[2+_0] = 0;
 			L[2+_1] = bits[0];
 			break;
 			}
 		j = 32 - i;
-		L[_0] = (bits[2] << i | bits[1] >> j) & 0xfffff
+		L[_0] = (((bits[2] << i) | (bits[1] >> j)) & 0xfffff)
-				| j + 1 << 20;
+				| ((j + 1) << 20);
 		L[_1] = (bits[1] << i | bits[0] >> j) & 0xffffffffL;
 		L[2+_0] = 0;
-		L[2+_1] = bits[0] & (1L << j) - 1;
+		L[2+_1] = bits[0] & ((1L << j) - 1);
 		break;
 
 	  hardly_normal:
 		j = 11 - hi0bits(bits[1]);
 		i = 32 - j;
-		L[_0] = bits[1] >> j & 0xfffff | j + 1 << 20;
+		L[_0] = (bits[1] >> j & 0xfffff) | ((j + 1) << 20);
 		L[_1] = (bits[1] << i | bits[0] >> j) & 0xffffffffL;
 		L[2+_0] = 0;
-		L[2+_1] = bits[0] & ((ULong)1L << j) - 1;
+		L[2+_1] = bits[0] & (((ULong)1L << j) - 1);
 		break;
 
 	  case STRTOG_Infinite:
@@ -153,16 +155,17 @@ ULtodd(ULong *L, ULong *bits, Long exp, int k)
 		break;
 
 	  case STRTOG_NaN:
-		L[0] = L[2] = d_QNAN0;
+		L[_0] = L[_0+2] = NanDflt_d_D2A[1];
-		L[1] = L[3] = d_QNAN1;
+		L[_1] = L[_1+2] = NanDflt_d_D2A[0];
 		break;
 
 	  case STRTOG_NaNbits:
-		L[_1] = (bits[1] >> 21 | bits[2] << 11) & (ULong)0xffffffffL;
+		L[_1] = (bits[1] >> 20 | bits[2] << 12) & (ULong)0xffffffffL;
-		L[_0] = bits[2] >> 21 | bits[3] << 11
+		L[_0] = bits[2] >> 20 | bits[3] << 12;
-			  | (ULong)0x7ff00000L;
+		L[_0] |= (L[_1] | L[_0]) ? (ULong)0x7ff00000L : (ULong)0x7ff80000L;
-		L[2+_1] = bits[0];
+		L[2+_1] = bits[0] & (ULong)0xffffffffL;
-		L[2+_0] = bits[1] | (ULong)0x7ff00000L;
+		L[2+_0] = bits[1] & 0xfffffL;
+		L[2+_0] |= (L[2+_1] | L[2+_0]) ? (ULong)0x7ff00000L : (ULong)0x7ff80000L;
 	  }
 	if (k & STRTOG_Neg) {
 		L[_0] |= 0x80000000L;
@@ -178,12 +181,11 @@ strtordd(CONST char *s, char **sp, int rounding, double *dd)
 #endif
 {
 #ifdef Sudden_Underflow
-	static CONST FPI fpi0 = { 106, 1-1023, 2046-1023-106+1, 1, 1 };
+	static FPI fpi0 = { 106, 1-1023, 2046-1023-106+1, 1, 1 };
 #else
-	static CONST FPI fpi0 = { 106, 1-1023-53+1, 2046-1023-106+1, 1, 0 };
+	static FPI fpi0 = { 106, 1-1023-53+1, 2046-1023-106+1, 1, 0 };
 #endif
-	CONST FPI *fpi;
+	FPI *fpi, fpi1;
-	FPI fpi1;
 	ULong bits[4];
 	Long exp;
 	int k;

gdtoa/strtorf.c

@@ -31,6 +31,8 @@ THIS SOFTWARE.
 
 #include "gdtoaimp.h"
 
+ extern ULong NanDflt_f_D2A[1];
+
  void
 #ifdef KR_headers
 ULtof(L, bits, exp, k) ULong *L; ULong *bits; Long exp; int k;
@@ -46,7 +48,7 @@ ULtof(ULong *L, ULong *bits, Long exp, int k)
 
 	  case STRTOG_Normal:
 	  case STRTOG_NaNbits:
-		L[0] = bits[0] & 0x7fffff | exp + 0x7f + 23 << 23;
+		L[0] = (bits[0] & 0x7fffff) | ((exp + 0x7f + 23) << 23);
 		break;
 
 	  case STRTOG_Denormal:
@@ -58,7 +60,7 @@ ULtof(ULong *L, ULong *bits, Long exp, int k)
 		break;
 
 	  case STRTOG_NaN:
-		L[0] = f_QNAN;
+		L[0] = NanDflt_f_D2A[0];
 	  }
 	if (k & STRTOG_Neg)
 		L[0] |= 0x80000000L;
@@ -71,9 +73,8 @@ strtorf(s, sp, rounding, f) CONST char *s; char **sp; int rounding; float *f;
 strtorf(CONST char *s, char **sp, int rounding, float *f)
 #endif
 {
-	static CONST FPI fpi0 = { 24, 1-127-24+1,  254-127-24+1, 1, SI };
+	static FPI fpi0 = { 24, 1-127-24+1,  254-127-24+1, 1, SI };
-	CONST FPI *fpi;
+	FPI *fpi, fpi1;
-	FPI fpi1;
 	ULong bits[1];
 	Long exp;
 	int k;

gdtoa/strtorx.c

@@ -51,6 +51,8 @@ THIS SOFTWARE.
 #define _4 0
 #endif
 
+ extern UShort NanDflt_ldus_D2A[5];
+
  void
 #ifdef KR_headers
 ULtox(L, bits, exp, k) UShort *L; ULong *bits; Long exp; int k;
@@ -80,15 +82,16 @@ ULtox(UShort *L, ULong *bits, Long exp, int k)
 
 	  case STRTOG_Infinite:
 		L[_0] = 0x7fff;
-		L[_1] = L[_2] = L[_3] = L[_4] = 0;
+		L[_1] = 0x8000;
+		L[_2] = L[_3] = L[_4] = 0;
 		break;
 
 	  case STRTOG_NaN:
-		L[0] = ldus_QNAN0;
+		L[_4] = NanDflt_ldus_D2A[0];
-		L[1] = ldus_QNAN1;
+		L[_3] = NanDflt_ldus_D2A[1];
-		L[2] = ldus_QNAN2;
+		L[_2] = NanDflt_ldus_D2A[2];
-		L[3] = ldus_QNAN3;
+		L[_1] = NanDflt_ldus_D2A[3];
-		L[4] = ldus_QNAN4;
+		L[_0] = NanDflt_ldus_D2A[4];
 	  }
 	if (k & STRTOG_Neg)
 		L[_0] |= 0x8000;
@@ -101,9 +104,8 @@ strtorx(s, sp, rounding, L) CONST char *s; char **sp; int rounding; void *L;
 strtorx(CONST char *s, char **sp, int rounding, void *L)
 #endif
 {
-	static CONST FPI fpi0 = { 64, 1-16383-64+1, 32766 - 16383 - 64 + 1, 1, SI };
+	static FPI fpi0 = { 64, 1-16383-64+1, 32766 - 16383 - 64 + 1, 1, SI };
-	CONST FPI *fpi;
+	FPI *fpi, fpi1;
-	FPI fpi1;
 	ULong bits[2];
 	Long exp;
 	int k;

gdtoa/strtorxL.c

@@ -47,6 +47,8 @@ THIS SOFTWARE.
 #define _2 0
 #endif
 
+ extern ULong NanDflt_xL_D2A[3];
+
  void
 #ifdef KR_headers
 ULtoxL(L, bits, exp, k) ULong *L; ULong *bits; Long exp; int k;
@@ -69,14 +71,15 @@ ULtoxL(ULong *L, ULong *bits, Long exp, int k)
 		break;
 
 	  case STRTOG_Infinite:
-		L[_0] = 0x7fff << 16;
+		L[_0] = 0x7fff0000;
-		L[_1] = L[_2] = 0;
+		L[_1] = 0x80000000;
+		L[_2] = 0;
 		break;
 
 	  case STRTOG_NaN:
-		L[0] = ld_QNAN0;
+		L[_0] = NanDflt_xL_D2A[2];
-		L[1] = ld_QNAN1;
+		L[_1] = NanDflt_xL_D2A[1];
-		L[2] = ld_QNAN2;
+		L[_2] = NanDflt_xL_D2A[0];
 	  }
 	if (k & STRTOG_Neg)
 		L[_0] |= 0x80000000L;
@@ -89,9 +92,8 @@ strtorxL(s, sp, rounding, L) CONST char *s; char **sp; int rounding; void *L;
 strtorxL(CONST char *s, char **sp, int rounding, void *L)
 #endif
 {
-	static CONST FPI fpi0 = { 64, 1-16383-64+1, 32766 - 16383 - 64 + 1, 1, SI };
+	static FPI fpi0 = { 64, 1-16383-64+1, 32766 - 16383 - 64 + 1, 1, SI };
-	CONST FPI *fpi;
+	FPI *fpi, fpi1;
-	FPI fpi1;
 	ULong bits[2];
 	Long exp;
 	int k;

gdtoa/ulp.c

@@ -34,15 +34,14 @@ THIS SOFTWARE.
  double
 ulp
 #ifdef KR_headers
-	(x) double x;
+	(x) U *x;
 #else
-	(double _x)
+	(U *x)
 #endif
 {
 	Long L;
-	U x, a;
+	U a;
 
-	dval(x) = _x;
 	L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
 #ifndef Sudden_Underflow
 	if (L > 0) {
@@ -50,22 +49,22 @@ ulp
 #ifdef IBM
 		L |= Exp_msk1 >> 4;
 #endif
-		word0(a) = L;
+		word0(&a) = L;
-		word1(a) = 0;
+		word1(&a) = 0;
 #ifndef Sudden_Underflow
 		}
 	else {
 		L = -L >> Exp_shift;
 		if (L < Exp_shift) {
-			word0(a) = 0x80000 >> L;
+			word0(&a) = 0x80000 >> L;
-			word1(a) = 0;
+			word1(&a) = 0;
 			}
 		else {
-			word0(a) = 0;
+			word0(&a) = 0;
 			L -= Exp_shift;
-			word1(a) = L >= 31 ? 1 : 1 << 31 - L;
+			word1(&a) = L >= 31 ? 1 : 1 << (31 - L);
 			}
 		}
 #endif
-	return dval(a);
+	return dval(&a);
 	}