2006-06-22 Jeff Johnston <jjohnstn@redhat.com>

        * libc/stdlib/Makefile.am: Add new gdtoa routines.
        * libc/stdlib/Makefile.in: Regenerated.
        * libc/stdlib/gd_qnan.h: New file.
        * libc/stdlib/gdtoa-gethex.c: Ditto.
        * libc/stdlib/gdtoa-hexnan.c: Ditto.
        * libc/stdlib/gdtoa.h: Ditto.
        * libc/stdlib/mprec.c: Add new helper routines needed by
        the new gdtoa code.
        * libc/stdlib/mprec.h: Integrate some defines and prototypes
        used by gdtoa routines here.
        * libc/stdlib/strtod.c: Rebased on David M. Gay's gdtoa-strtod.c
        which adds C99 support such as nan, inf, and hexadecimal input
        format.

9 files changed, 1816 insertions, 635 deletions

diff --git a/newlib/libc/stdlib/Makefile.am b/newlib/libc/stdlib/Makefile.am
index fdc52e04e..0874729c5 100644
--- a/newlib/libc/stdlib/Makefile.am
+++ b/newlib/libc/stdlib/Makefile.am
@@ -27,6 +27,8 @@ GENERAL_SOURCES = \
 	envlock.c	\
 	eprintf.c	\
 	exit.c 		\
+	gdtoa-gethex.c	\
+	gdtoa-hexnan.c	\
 	getenv.c  	\
 	getenv_r.c	\
 	labs.c 		\
@@ -256,6 +258,8 @@ $(lpfx)mbtowc_r.$(oext): mbtowc_r.c mbctype.h
 
 $(lpfx)mprec.$(oext): mprec.c mprec.h
 $(lpfx)strtod.$(oext): strtod.c mprec.h
+$(lpfx)gdtoa-gethex.$(oext): gdtoa-gethex.c mprec.h
+$(lpfx)gdtoa-hexnan.$(oext): gdtoa-hexnan.c mprec.h
 $(lpfx)wctomb_r.$(oext): wctomb_r.c mbctype.h
 $(lpfx)drand48.$(oext): drand48.c rand48.h
 $(lpfx)erand48.$(oext): erand48.c rand48.h
diff --git a/newlib/libc/stdlib/Makefile.in b/newlib/libc/stdlib/Makefile.in
index 56f4a365f..0c6c5e75a 100644
--- a/newlib/libc/stdlib/Makefile.in
+++ b/newlib/libc/stdlib/Makefile.in
@@ -74,6 +74,7 @@ am__objects_1 = lib_a-__adjust.$(OBJEXT) lib_a-__atexit.$(OBJEXT) \
 	lib_a-dtoa.$(OBJEXT) lib_a-dtoastub.$(OBJEXT) \
 	lib_a-environ.$(OBJEXT) lib_a-envlock.$(OBJEXT) \
 	lib_a-eprintf.$(OBJEXT) lib_a-exit.$(OBJEXT) \
+	lib_a-gdtoa-gethex.$(OBJEXT) lib_a-gdtoa-hexnan.$(OBJEXT) \
 	lib_a-getenv.$(OBJEXT) lib_a-getenv_r.$(OBJEXT) \
 	lib_a-labs.$(OBJEXT) lib_a-ldiv.$(OBJEXT) \
 	lib_a-ldtoa.$(OBJEXT) lib_a-malloc.$(OBJEXT) \
@@ -124,12 +125,12 @@ LTLIBRARIES = $(noinst_LTLIBRARIES)
 am__objects_7 = __adjust.lo __atexit.lo __call_atexit.lo __exp10.lo \
 	__ten_mu.lo _Exit.lo abort.lo abs.lo assert.lo atexit.lo \
 	atof.lo atoff.lo atoi.lo atol.lo calloc.lo div.lo dtoa.lo \
-	dtoastub.lo environ.lo envlock.lo eprintf.lo exit.lo getenv.lo \
-	getenv_r.lo labs.lo ldiv.lo ldtoa.lo malloc.lo mblen.lo \
-	mblen_r.lo mbstowcs.lo mbstowcs_r.lo mbtowc.lo mbtowc_r.lo \
-	mlock.lo mprec.lo mstats.lo rand.lo rand_r.lo realloc.lo \
-	strtod.lo strtol.lo strtoul.lo wcstombs.lo wcstombs_r.lo \
-	wctomb.lo wctomb_r.lo
+	dtoastub.lo environ.lo envlock.lo eprintf.lo exit.lo \
+	gdtoa-gethex.lo gdtoa-hexnan.lo getenv.lo getenv_r.lo labs.lo \
+	ldiv.lo ldtoa.lo malloc.lo mblen.lo mblen_r.lo mbstowcs.lo \
+	mbstowcs_r.lo mbtowc.lo mbtowc_r.lo mlock.lo mprec.lo \
+	mstats.lo rand.lo rand_r.lo realloc.lo strtod.lo strtol.lo \
+	strtoul.lo wcstombs.lo wcstombs_r.lo wctomb.lo wctomb_r.lo
 am__objects_8 = cxa_atexit.lo cxa_finalize.lo drand48.lo ecvtbuf.lo \
 	efgcvt.lo erand48.lo jrand48.lo lcong48.lo lrand48.lo \
 	mrand48.lo msize.lo mtrim.lo nrand48.lo rand48.lo seed48.lo \
@@ -249,6 +250,7 @@ PACKAGE_TARNAME = @PACKAGE_TARNAME@
 PACKAGE_VERSION = @PACKAGE_VERSION@
 PATH_SEPARATOR = @PATH_SEPARATOR@
 RANLIB = @RANLIB@
+READELF = @READELF@
 SET_MAKE = @SET_MAKE@
 SHELL = @SHELL@
 STRIP = @STRIP@
@@ -259,6 +261,7 @@ ac_ct_AR = @ac_ct_AR@
 ac_ct_AS = @ac_ct_AS@
 ac_ct_CC = @ac_ct_CC@
 ac_ct_RANLIB = @ac_ct_RANLIB@
+ac_ct_READELF = @ac_ct_READELF@
 ac_ct_STRIP = @ac_ct_STRIP@
 aext = @aext@
 am__fastdepCC_FALSE = @am__fastdepCC_FALSE@
@@ -329,6 +332,8 @@ GENERAL_SOURCES = \
 	envlock.c	\
 	eprintf.c	\
 	exit.c 		\
+	gdtoa-gethex.c	\
+	gdtoa-hexnan.c	\
 	getenv.c  	\
 	getenv_r.c	\
 	labs.c 		\
@@ -685,6 +690,18 @@ lib_a-exit.o: exit.c
 lib_a-exit.obj: exit.c
 	$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-exit.obj `if test -f 'exit.c'; then $(CYGPATH_W) 'exit.c'; else $(CYGPATH_W) '$(srcdir)/exit.c'; fi`
 
+lib_a-gdtoa-gethex.o: gdtoa-gethex.c
+	$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-gdtoa-gethex.o `test -f 'gdtoa-gethex.c' || echo '$(srcdir)/'`gdtoa-gethex.c
+
+lib_a-gdtoa-gethex.obj: gdtoa-gethex.c
+	$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-gdtoa-gethex.obj `if test -f 'gdtoa-gethex.c'; then $(CYGPATH_W) 'gdtoa-gethex.c'; else $(CYGPATH_W) '$(srcdir)/gdtoa-gethex.c'; fi`
+
+lib_a-gdtoa-hexnan.o: gdtoa-hexnan.c
+	$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-gdtoa-hexnan.o `test -f 'gdtoa-hexnan.c' || echo '$(srcdir)/'`gdtoa-hexnan.c
+
+lib_a-gdtoa-hexnan.obj: gdtoa-hexnan.c
+	$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-gdtoa-hexnan.obj `if test -f 'gdtoa-hexnan.c'; then $(CYGPATH_W) 'gdtoa-hexnan.c'; else $(CYGPATH_W) '$(srcdir)/gdtoa-hexnan.c'; fi`
+
 lib_a-getenv.o: getenv.c
 	$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-getenv.o `test -f 'getenv.c' || echo '$(srcdir)/'`getenv.c
 
@@ -1298,6 +1315,8 @@ $(lpfx)mbtowc_r.$(oext): mbtowc_r.c mbctype.h
 
 $(lpfx)mprec.$(oext): mprec.c mprec.h
 $(lpfx)strtod.$(oext): strtod.c mprec.h
+$(lpfx)gdtoa-gethex.$(oext): gdtoa-gethex.c mprec.h
+$(lpfx)gdtoa-hexnan.$(oext): gdtoa-hexnan.c mprec.h
 $(lpfx)wctomb_r.$(oext): wctomb_r.c mbctype.h
 $(lpfx)drand48.$(oext): drand48.c rand48.h
 $(lpfx)erand48.$(oext): erand48.c rand48.h
diff --git a/newlib/libc/stdlib/gd_qnan.h b/newlib/libc/stdlib/gd_qnan.h
new file mode 100644
index 000000000..68f16cd12
--- /dev/null
+++ b/newlib/libc/stdlib/gd_qnan.h
@@ -0,0 +1,33 @@
+#ifdef __IEEE_BIG_ENDIAN
+
+#define f_QNAN 0x7fc00000
+#define d_QNAN0 0x7ff80000
+#define d_QNAN1 0x0
+#define ld_QNAN0 0x7ff80000
+#define ld_QNAN1 0x0
+#define ld_QNAN2 0x0
+#define ld_QNAN3 0x0
+#define ldus_QNAN0 0x7ff8
+#define ldus_QNAN1 0x0
+#define ldus_QNAN2 0x0
+#define ldus_QNAN3 0x0
+#define ldus_QNAN4 0x0
+
+#elif defined(__IEEE_LITTLE_ENDIAN)
+
+#define f_QNAN 0xffc00000
+#define d_QNAN0 0x0
+#define d_QNAN1 0xfff80000
+#define ld_QNAN0 0x0
+#define ld_QNAN1 0xc0000000
+#define ld_QNAN2 0xffff
+#define ld_QNAN3 0x0
+#define ldus_QNAN0 0x0
+#define ldus_QNAN1 0x0
+#define ldus_QNAN2 0x0
+#define ldus_QNAN3 0xc000
+#define ldus_QNAN4 0xffff
+
+#else
+#error IEEE endian not defined
+#endif
diff --git a/newlib/libc/stdlib/gdtoa-gethex.c b/newlib/libc/stdlib/gdtoa-gethex.c
new file mode 100644
index 000000000..92f30fca0
--- /dev/null
+++ b/newlib/libc/stdlib/gdtoa-gethex.c
@@ -0,0 +1,351 @@
+/****************************************************************
+
+The author of this software is David M. Gay.
+
+Copyright (C) 1998 by Lucent Technologies
+All Rights Reserved
+
+Permission to use, copy, modify, and distribute this software and
+its documentation for any purpose and without fee is hereby
+granted, provided that the above copyright notice appear in all
+copies and that both that the copyright notice and this
+permission notice and warranty disclaimer appear in supporting
+documentation, and that the name of Lucent or any of its entities
+not be used in advertising or publicity pertaining to
+distribution of the software without specific, written prior
+permission.
+
+LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
+INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
+IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
+SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
+IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
+ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
+THIS SOFTWARE.
+
+****************************************************************/
+
+/* Please send bug reports to David M. Gay (dmg at acm dot org,
+ * with " at " changed at "@" and " dot " changed to ".").	*/
+
+#include <_ansi.h>
+#include <reent.h>
+#include <string.h>
+#include "mprec.h"
+#include "gdtoa.h"
+#include "gd_qnan.h"
+
+#ifdef USE_LOCALE
+#include "locale.h"
+#endif
+
+unsigned char hexdig[256];
+
+static void
+_DEFUN (htinit, (h, s, inc),
+	unsigned char *h _AND
+	unsigned char *s _AND
+	int inc)
+{
+	int i, j;
+	for(i = 0; (j = s[i]) !=0; i++)
+		h[j] = i + inc;
+}
+
+void
+_DEFUN_VOID (hexdig_init)
+{
+#define USC (unsigned char *)
+	htinit(hexdig, USC "0123456789", 0x10);
+	htinit(hexdig, USC "abcdef", 0x10 + 10);
+	htinit(hexdig, USC "ABCDEF", 0x10 + 10);
+}
+
+static void
+_DEFUN(rshift, (b, k),
+	_Bigint *b _AND
+	int k)
+{
+	__ULong *x, *x1, *xe, y;
+	int n;
+
+	x = x1 = b->_x;
+	n = k >> kshift;
+	if (n < b->_wds) {
+		xe = x + b->_wds;
+		x += n;
+		if (k &= kmask) {
+			n = ULbits - k;
+			y = *x++ >> k;
+			while(x < xe) {
+				*x1++ = (y | (*x << n)) & ALL_ON;
+				y = *x++ >> k;
+				}
+			if ((*x1 = y) !=0)
+				x1++;
+			}
+		else
+			while(x < xe)
+				*x1++ = *x++;
+		}
+	if ((b->_wds = x1 - b->_x) == 0)
+		b->_x[0] = 0;
+}
+
+static _Bigint *
+_DEFUN (increment, (ptr, b),
+	struct _reent *ptr _AND
+	_Bigint *b)
+{
+	__ULong *x, *xe;
+	_Bigint *b1;
+#ifdef Pack_16
+	__ULong carry = 1, y;
+#endif
+
+	x = b->_x;
+	xe = x + b->_wds;
+#ifdef Pack_32
+	do {
+		if (*x < (__ULong)0xffffffffL) {
+			++*x;
+			return b;
+			}
+		*x++ = 0;
+		} while(x < xe);
+#else
+	do {
+		y = *x + carry;
+		carry = y >> 16;
+		*x++ = y & 0xffff;
+		if (!carry)
+			return b;
+		} while(x < xe);
+	if (carry)
+#endif
+	{
+		if (b->_wds >= b->_maxwds) {
+			b1 = Balloc(ptr, b->_k+1);
+			Bcopy(b1, b);
+			Bfree(ptr, b);
+			b = b1;
+			}
+		b->_x[b->_wds++] = 1;
+		}
+	return b;
+}
+
+
+int
+_DEFUN(gethex, (ptr, sp, fpi, exp, bp, sign),
+	struct _reent *ptr _AND
+	_CONST char **sp _AND
+	FPI *fpi _AND
+	Long *exp _AND
+	_Bigint **bp _AND
+	int sign)
+{
+	_Bigint *b;
+	_CONST unsigned char *decpt, *s0, *s, *s1;
+	int esign, havedig, irv, k, n, nbits, up, zret;
+	__ULong L, lostbits, *x;
+	Long e, e1;
+#ifdef USE_LOCALE
+	unsigned char decimalpoint = *localeconv()->decimal_point;
+#else
+#define decimalpoint '.'
+#endif
+
+	if (!hexdig['0'])
+		hexdig_init();
+	havedig = 0;
+	s0 = *(_CONST unsigned char **)sp + 2;
+	while(s0[havedig] == '0')
+		havedig++;
+	s0 += havedig;
+	s = s0;
+	decpt = 0;
+	zret = 0;
+	e = 0;
+	if (!hexdig[*s]) {
+		zret = 1;
+		if (*s != decimalpoint)
+			goto pcheck;
+		decpt = ++s;
+		if (!hexdig[*s])
+			goto pcheck;
+		while(*s == '0')
+			s++;
+		if (hexdig[*s])
+			zret = 0;
+		havedig = 1;
+		s0 = s;
+		}
+	while(hexdig[*s])
+		s++;
+	if (*s == decimalpoint && !decpt) {
+		decpt = ++s;
+		while(hexdig[*s])
+			s++;
+		}
+	if (decpt)
+		e = -(((Long)(s-decpt)) << 2);
+ pcheck:
+	s1 = s;
+	switch(*s) {
+	  case 'p':
+	  case 'P':
+		esign = 0;
+		switch(*++s) {
+		  case '-':
+			esign = 1;
+			/* no break */
+		  case '+':
+			s++;
+		  }
+		if ((n = hexdig[*s]) == 0 || n > 0x19) {
+			s = s1;
+			break;
+			}
+		e1 = n - 0x10;
+		while((n = hexdig[*++s]) !=0 && n <= 0x19)
+			e1 = 10*e1 + n - 0x10;
+		if (esign)
+			e1 = -e1;
+		e += e1;
+	  }
+	*sp = (char*)s;
+	if (zret)
+		return havedig ? STRTOG_Zero : STRTOG_NoNumber;
+	n = s1 - s0 - 1;
+	for(k = 0; n > 7; n >>= 1)
+		k++;
+	b = Balloc(ptr, k);
+	x = b->_x;
+	n = 0;
+	L = 0;
+	while(s1 > s0) {
+		if (*--s1 == decimalpoint)
+			continue;
+		if (n == 32) {
+			*x++ = L;
+			L = 0;
+			n = 0;
+			}
+		L |= (hexdig[*s1] & 0x0f) << n;
+		n += 4;
+		}
+	*x++ = L;
+	b->_wds = n = x - b->_x;
+	n = 32*n - hi0bits(L);
+	nbits = fpi->nbits;
+	lostbits = 0;
+	x = b->_x;
+	if (n > nbits) {
+		n -= nbits;
+		if (any_on(b,n)) {
+			lostbits = 1;
+			k = n - 1;
+			if (x[k>>kshift] & 1 << (k & kmask)) {
+				lostbits = 2;
+				if (k > 1 && any_on(b,k-1))
+					lostbits = 3;
+				}
+			}
+		rshift(b, n);
+		e += n;
+		}
+	else if (n < nbits) {
+		n = nbits - n;
+		b = lshift(ptr, b, n);
+		e -= n;
+		x = b->_x;
+		}
+	if (e > fpi->emax) {
+ ovfl:
+		Bfree(ptr, b);
+		*bp = 0;
+		return STRTOG_Infinite | STRTOG_Overflow | STRTOG_Inexhi;
+		}
+	irv = STRTOG_Normal;
+	if (e < fpi->emin) {
+		irv = STRTOG_Denormal;
+		n = fpi->emin - e;
+		if (n >= nbits) {
+			switch (fpi->rounding) {
+			  case FPI_Round_near:
+				if (n == nbits && (n < 2 || any_on(b,n-1)))
+					goto one_bit;
+				break;
+			  case FPI_Round_up:
+				if (!sign)
+					goto one_bit;
+				break;
+			  case FPI_Round_down:
+				if (sign) {
+ one_bit:
+					*exp = fpi->emin;
+					x[0] = b->_wds = 1;
+					*bp = b;
+					return STRTOG_Denormal | STRTOG_Inexhi
+						| STRTOG_Underflow;
+					}
+			  }
+			Bfree(ptr, b);
+			*bp = 0;
+			return STRTOG_Zero | STRTOG_Inexlo | STRTOG_Underflow;
+			}
+		k = n - 1;
+		if (lostbits)
+			lostbits = 1;
+		else if (k > 0)
+			lostbits = any_on(b,k);
+		if (x[k>>kshift] & 1 << (k & kmask))
+			lostbits |= 2;
+		nbits -= n;
+		rshift(b,n);
+		e = fpi->emin;
+		}
+	if (lostbits) {
+		up = 0;
+		switch(fpi->rounding) {
+		  case FPI_Round_zero:
+			break;
+		  case FPI_Round_near:
+		    if ((lostbits & 2)
+			    && ((lostbits & 1) | (x[0] & 1)))
+				up = 1;
+			break;
+		  case FPI_Round_up:
+			up = 1 - sign;
+			break;
+		  case FPI_Round_down:
+			up = sign;
+		  }
+		if (up) {
+			k = b->_wds;
+			b = increment(ptr, b);
+			x = b->_x;
+			if (irv == STRTOG_Denormal) {
+				if (nbits == fpi->nbits - 1
+				 && x[nbits >> kshift] & 1 << (nbits & kmask))
+					irv =  STRTOG_Normal;
+				}
+			else if ((b->_wds > k)
+			 || ((n = nbits & kmask) !=0
+			     && (hi0bits(x[k-1]) < 32-n))) {
+				rshift(b,1);
+				if (++e > fpi->emax)
+					goto ovfl;
+				}
+			irv |= STRTOG_Inexhi;
+			}
+		else
+			irv |= STRTOG_Inexlo;
+		}
+	*bp = b;
+	*exp = e;
+	return irv;
+}
+
diff --git a/newlib/libc/stdlib/gdtoa-hexnan.c b/newlib/libc/stdlib/gdtoa-hexnan.c
new file mode 100644
index 000000000..058bbd17d
--- /dev/null
+++ b/newlib/libc/stdlib/gdtoa-hexnan.c
@@ -0,0 +1,142 @@
+/****************************************************************
+
+The author of this software is David M. Gay.
+
+Copyright (C) 2000 by Lucent Technologies
+All Rights Reserved
+
+Permission to use, copy, modify, and distribute this software and
+its documentation for any purpose and without fee is hereby
+granted, provided that the above copyright notice appear in all
+copies and that both that the copyright notice and this
+permission notice and warranty disclaimer appear in supporting
+documentation, and that the name of Lucent or any of its entities
+not be used in advertising or publicity pertaining to
+distribution of the software without specific, written prior
+permission.
+
+LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
+INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
+IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
+SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
+IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
+ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
+THIS SOFTWARE.
+
+****************************************************************/
+
+/* Please send bug reports to
+	David M. Gay
+	Bell Laboratories, Room 2C-463
+	600 Mountain Avenue
+	Murray Hill, NJ 07974-0636
+	U.S.A.
+	dmg@bell-labs.com
+ */
+
+/* Modified 06-21-2006 by Jeff Johnston to work with newlib.  */
+
+#include <_ansi.h>
+#include <reent.h>
+#include <string.h>
+#include "mprec.h"
+#include "gdtoa.h"
+
+#ifdef INFNAN_CHECK
+static void
+_DEFUN (L_shift, (x, x1, i),
+	__ULong *x _AND
+	__ULong *x1 _AND
+	int i)
+{
+	int j;
+
+	i = 8 - i;
+	i <<= 2;
+	j = ULbits - i;
+	do {
+		*x |= x[1] << j;
+		x[1] >>= i;
+		} while(++x < x1);
+}
+
+int
+_DEFUN (hexnan, (sp, fpi, x0),
+	_CONST char **sp _AND
+	FPI *fpi _AND
+	__ULong *x0)
+{
+	__ULong c, h, *x, *x1, *xe;
+	_CONST char *s;
+	int havedig, hd0, i, nbits;
+
+	if (!hexdig['0'])
+		hexdig_init();
+	nbits = fpi->nbits;
+	x = x0 + (nbits >> kshift);
+	if (nbits & kmask)
+		x++;
+	*--x = 0;
+	x1 = xe = x;
+	havedig = hd0 = i = 0;
+	s = *sp;
+	while((c = *(_CONST unsigned char*)++s)) {
+		if (!(h = hexdig[c])) {
+			if (c <= ' ') {
+				if (hd0 < havedig) {
+					if (x < x1 && i < 8)
+						L_shift(x, x1, i);
+					if (x <= x0) {
+						i = 8;
+						continue;
+						}
+					hd0 = havedig;
+					*--x = 0;
+					x1 = x;
+					i = 0;
+					}
+				continue;
+				}
+			if (/*(*/ c == ')' && havedig) {
+				*sp = s + 1;
+				break;
+				}
+			return STRTOG_NaN;
+			}
+		havedig++;
+		if (++i > 8) {
+			if (x <= x0)
+				continue;
+			i = 1;
+			*--x = 0;
+			}
+		*x = ((*x << 4) | (h & 0xf));
+		}
+	if (!havedig)
+		return STRTOG_NaN;
+	if (x < x1 && i < 8)
+		L_shift(x, x1, i);
+	if (x > x0) {
+		x1 = x0;
+		do *x1++ = *x++;
+			while(x <= xe);
+		do *x1++ = 0;
+			while(x1 <= xe);
+		}
+	else {
+		/* truncate high-order word if necessary */
+		if ( (i = nbits & (ULbits-1)) !=0)
+			*xe &= ((__ULong)0xffffffff) >> (ULbits - i);
+		}
+	for(x1 = xe;; --x1) {
+		if (*x1 != 0)
+			break;
+		if (x1 == x0) {
+			*x1 = 1;
+			break;
+			}
+		}
+	return STRTOG_NaNbits;
+}
+#endif /* INFNAN_CHECK */
diff --git a/newlib/libc/stdlib/gdtoa.h b/newlib/libc/stdlib/gdtoa.h
new file mode 100644
index 000000000..5029e58de
--- /dev/null
+++ b/newlib/libc/stdlib/gdtoa.h
@@ -0,0 +1,72 @@
+/****************************************************************
+
+The author of this software is David M. Gay.
+
+Copyright (C) 1998 by Lucent Technologies
+All Rights Reserved
+
+Permission to use, copy, modify, and distribute this software and
+its documentation for any purpose and without fee is hereby
+granted, provided that the above copyright notice appear in all
+copies and that both that the copyright notice and this
+permission notice and warranty disclaimer appear in supporting
+documentation, and that the name of Lucent or any of its entities
+not be used in advertising or publicity pertaining to
+distribution of the software without specific, written prior
+permission.
+
+LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
+INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
+IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
+SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
+IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
+ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
+THIS SOFTWARE.
+
+****************************************************************/
+
+/* Please send bug reports to David M. Gay (dmg at acm dot org,
+ * with " at " changed at "@" and " dot " changed to ".").	*/
+
+#ifndef GDTOA_H_INCLUDED
+#define GDTOA_H_INCLUDED
+
+
+ enum {	/* return values from strtodg */
+	STRTOG_Zero	= 0,
+	STRTOG_Normal	= 1,
+	STRTOG_Denormal	= 2,
+	STRTOG_Infinite	= 3,
+	STRTOG_NaN	= 4,
+	STRTOG_NaNbits	= 5,
+	STRTOG_NoNumber	= 6,
+	STRTOG_Retmask	= 7,
+
+	/* The following may be or-ed into one of the above values. */
+
+	STRTOG_Neg	= 0x08,
+	STRTOG_Inexlo	= 0x10,
+	STRTOG_Inexhi	= 0x20,
+	STRTOG_Inexact	= 0x30,
+	STRTOG_Underflow= 0x40,
+	STRTOG_Overflow	= 0x80
+	};
+
+ typedef struct
+FPI {
+	int nbits;
+	int emin;
+	int emax;
+	int rounding;
+	int sudden_underflow;
+	} FPI;
+
+enum {	/* FPI.rounding values: same as FLT_ROUNDS */
+	FPI_Round_zero = 0,
+	FPI_Round_near = 1,
+	FPI_Round_up = 2,
+	FPI_Round_down = 3
+	};
+
+#endif /* GDTOA_H_INCLUDED */
diff --git a/newlib/libc/stdlib/mprec.c b/newlib/libc/stdlib/mprec.c
index 0ef28c745..6e84ece5b 100644
--- a/newlib/libc/stdlib/mprec.c
+++ b/newlib/libc/stdlib/mprec.c
@@ -985,3 +985,61 @@ _DEFUN (_mprec_log10, (dig),
     }
   return v;
 }
+
+void
+_DEFUN (copybits, (c, n, b),
+	__ULong *c _AND
+	int n _AND
+	_Bigint *b)
+{
+	__ULong *ce, *x, *xe;
+#ifdef Pack_16
+	int nw, nw1;
+#endif
+
+	ce = c + ((n-1) >> kshift) + 1;
+	x = b->_x;
+#ifdef Pack_32
+	xe = x + b->_wds;
+	while(x < xe)
+		*c++ = *x++;
+#else
+	nw = b->_wds;
+	nw1 = nw & 1;
+	for(xe = x + (nw - nw1); x < xe; x += 2)
+		Storeinc(c, x[1], x[0]);
+	if (nw1)
+		*c++ = *x;
+#endif
+	while(c < ce)
+		*c++ = 0;
+}
+
+__ULong
+_DEFUN (any_on, (b, k),
+	_Bigint *b _AND
+	int k)
+{
+	int n, nwds;
+	__ULong *x, *x0, x1, x2;
+
+	x = b->_x;
+	nwds = b->_wds;
+	n = k >> kshift;
+	if (n > nwds)
+		n = nwds;
+	else if (n < nwds && (k &= kmask)) {
+		x1 = x2 = x[n];
+		x1 >>= k;
+		x1 <<= k;
+		if (x1 != x2)
+			return 1;
+		}
+	x0 = x;
+	x += n;
+	while(x > x0)
+		if (*--x)
+			return 1;
+	return 0;
+}
+
diff --git a/newlib/libc/stdlib/mprec.h b/newlib/libc/stdlib/mprec.h
index 4ca48f22f..4fd496831 100644
--- a/newlib/libc/stdlib/mprec.h
+++ b/newlib/libc/stdlib/mprec.h
@@ -78,6 +78,39 @@ union double_union
 #define word1(x) (x.i[1])
 #endif
 
+
+/* The following is taken from gdtoaimp.h for use with new strtod.  */
+typedef __int32_t Long;
+typedef union { double d; __ULong L[2]; } U;
+
+#ifdef YES_ALIAS
+#define dval(x) x
+#ifdef IEEE_8087
+#define dword0(x) ((__ULong *)&x)[1]
+#define dword1(x) ((__ULong *)&x)[0]
+#else
+#define dword0(x) ((__ULong *)&x)[0]
+#define dword1(x) ((__ULong *)&x)[1]
+#endif
+#else /* !YES_ALIAS */
+#ifdef IEEE_8087
+#define dword0(x) ((U*)&x)->L[1]
+#define dword1(x) ((U*)&x)->L[0]
+#else
+#define dword0(x) ((U*)&x)->L[0]
+#define dword1(x) ((U*)&x)->L[1]
+#endif
+#define dval(x) ((U*)&x)->d
+#endif /* YES_ALIAS */
+
+
+#undef SI
+#ifdef Sudden_Underflow
+#define SI 1
+#else
+#define SI 0
+#endif
+
 /* The following definition of Storeinc is appropriate for MIPS processors.
  * An alternative that might be better on some machines is
  * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
@@ -161,12 +194,22 @@ union double_union
 #define Quick_max 14
 #define Int_max 14
 #define Infinite(x) (word0(x) == ((__uint32_t)0x7ff00000L)) /* sufficient test for here */
+
+#ifndef Flt_Rounds
+#ifdef FLT_ROUNDS
+#define Flt_Rounds FLT_ROUNDS
+#else
+#define Flt_Rounds 1
+#endif
+#endif /*Flt_Rounds*/
+
 #endif
 
 #else
 #undef  Sudden_Underflow
 #define Sudden_Underflow
 #ifdef IBM
+#define Flt_Rounds 0
 #define Exp_shift  24
 #define Exp_shift1 24
 #define Exp_msk1   ((__uint32_t)0x1000000L)
@@ -191,6 +234,7 @@ union double_union
 #define Quick_max 14
 #define Int_max 15
 #else /* VAX */
+#define Flt_Rounds 1
 #define Exp_shift  23
 #define Exp_shift1 7
 #define Exp_msk1    0x80
@@ -219,6 +263,58 @@ union double_union
 
 #ifndef IEEE_Arith
 #define ROUND_BIASED
+#else
+#define Scale_Bit 0x10
+#if defined(_DOUBLE_IS_32BITS) && defined(__v800)
+#define n_bigtens 2
+#else
+#define n_bigtens 5
+#endif
+#endif
+
+#ifdef IBM
+#define n_bigtens 3
+#endif
+
+#ifdef VAX
+#define n_bigtens 2
+#endif
+
+#ifndef __NO_INFNAN_CHECK
+#define INFNAN_CHECK
+#endif
+
+/*
+ * NAN_WORD0 and NAN_WORD1 are only referenced in strtod.c.  Prior to
+ * 20050115, they used to be hard-wired here (to 0x7ff80000 and 0,
+ * respectively), but now are determined by compiling and running
+ * qnan.c to generate gd_qnan.h, which specifies d_QNAN0 and d_QNAN1.
+ * Formerly gdtoaimp.h recommended supplying suitable -DNAN_WORD0=...
+ * and -DNAN_WORD1=...  values if necessary.  This should still work.
+ * (On HP Series 700/800 machines, -DNAN_WORD0=0x7ff40000 works.)
+ */
+#ifdef IEEE_Arith
+#ifdef IEEE_MC68k
+#define _0 0
+#define _1 1
+#ifndef NAN_WORD0
+#define NAN_WORD0 d_QNAN0
+#endif
+#ifndef NAN_WORD1
+#define NAN_WORD1 d_QNAN1
+#endif
+#else
+#define _0 1
+#define _1 0
+#ifndef NAN_WORD0
+#define NAN_WORD0 d_QNAN1
+#endif
+#ifndef NAN_WORD1
+#define NAN_WORD1 d_QNAN0
+#endif
+#endif
+#else
+#undef INFNAN_CHECK
 #endif
 
 #ifdef RND_PRODQUOT
@@ -249,6 +345,17 @@ extern double rnd_prod(double, double), rnd_quot(double, double);
 #endif
 #endif
 
+#ifdef Pack_32
+#define ULbits 32
+#define kshift 5
+#define kmask 31
+#define ALL_ON 0xffffffff
+#else
+#define ULbits 16
+#define kshift 4
+#define kmask 15
+#define ALL_ON 0xffff
+#endif
 
 #ifdef __cplusplus
 extern "C" double strtod(const char *s00, char **se);
@@ -261,26 +368,34 @@ typedef struct _Bigint _Bigint;
 
 #define Balloc	_Balloc
 #define Bfree	_Bfree
-#define multadd _multadd
-#define s2b	_s2b
-#define lo0bits _lo0bits
-#define hi0bits _hi0bits
-#define i2b	_i2b
-#define mult	_multiply
-#define pow5mult	_pow5mult
-#define lshift	_lshift
+#define multadd __multadd
+#define s2b	__s2b
+#define lo0bits __lo0bits
+#define hi0bits __hi0bits
+#define i2b	__i2b
+#define mult	__multiply
+#define pow5mult	__pow5mult
+#define lshift	__lshift
 #define cmp	__mcmp
 #define diff	__mdiff
-#define ulp 	_ulp
-#define b2d	_b2d
-#define d2b	_d2b
-#define ratio	_ratio
+#define ulp 	__ulp
+#define b2d	__b2d
+#define d2b	__d2b
+#define ratio	__ratio
+#define any_on	__any_on
+#define gethex  __gethex
+#define copybits 	__copybits
+#define hexnan	__hexnan
+#define hexdig_init 	__hexdig_init
+
+#define hexdig  __hexdig
 
 #define tens __mprec_tens
 #define bigtens __mprec_bigtens
 #define tinytens __mprec_tinytens
 
 struct _reent ;
+struct FPI;
 double 		_EXFUN(ulp,(double x));
 double		_EXFUN(b2d,(_Bigint *a , int *e));
 _Bigint *	_EXFUN(Balloc,(struct _reent *p, int k));
@@ -292,23 +407,25 @@ _Bigint *	_EXFUN(mult, (struct _reent *, _Bigint *, _Bigint *));
 _Bigint *	_EXFUN(pow5mult, (struct _reent *, _Bigint *, int k));
 int 		_EXFUN(hi0bits,(__ULong));
 int 		_EXFUN(lo0bits,(__ULong *));
-_Bigint *        _EXFUN(d2b,(struct _reent *p, double d, int *e, int *bits));
-_Bigint *        _EXFUN(lshift,(struct _reent *p, _Bigint *b, int k));
-_Bigint *        _EXFUN(diff,(struct _reent *p, _Bigint *a, _Bigint *b));
-int             _EXFUN(cmp,(_Bigint *a, _Bigint *b));
-
+_Bigint *	_EXFUN(d2b,(struct _reent *p, double d, int *e, int *bits));
+_Bigint *	_EXFUN(lshift,(struct _reent *p, _Bigint *b, int k));
+_Bigint *	_EXFUN(diff,(struct _reent *p, _Bigint *a, _Bigint *b));
+int		_EXFUN(cmp,(_Bigint *a, _Bigint *b));
+int		_EXFUN(gethex,(struct _reent *p, _CONST char **sp, struct FPI *fpi, Long *exp, _Bigint **bp, int sign));     
 double		_EXFUN(ratio,(_Bigint *a, _Bigint *b));
-#define Bcopy(x,y) memcpy((char *)&x->_sign, (char *)&y->_sign, y->_wds*sizeof(__Long) + 2*sizeof(int))
-
-#if defined(_DOUBLE_IS_32BITS) && defined(__v800)
-#define n_bigtens 2
-#else
-#define n_bigtens 5
+__ULong		_EXFUN(any_on,(_Bigint *b, int k));
+void		_EXFUN(copybits,(__ULong *c, int n, _Bigint *b));
+void		_EXFUN(hexdig_init,(void));
+#ifdef INFNAN_CHECK
+int		_EXFUN(hexnan,(_CONST char **sp, struct FPI *fpi, __ULong *x0));
 #endif
 
+#define Bcopy(x,y) memcpy((char *)&x->_sign, (char *)&y->_sign, y->_wds*sizeof(__Long) + 2*sizeof(int))
+
 extern _CONST double tinytens[];
 extern _CONST double bigtens[];
 extern _CONST double tens[];
+extern unsigned char hexdig[];
 
 
 double _EXFUN(_mprec_log10,(int));
diff --git a/newlib/libc/stdlib/strtod.c b/newlib/libc/stdlib/strtod.c
index 9b70dfc3c..57297e05b 100644
--- a/newlib/libc/stdlib/strtod.c
+++ b/newlib/libc/stdlib/strtod.c
@@ -70,37 +70,130 @@ Supporting OS subroutines required: <<close>>, <<fstat>>, <<isatty>>,
 */
 
 /****************************************************************
- *
- * The author of this software is David M. Gay.
- *
- * Copyright (c) 1991 by AT&T.
- *
- * Permission to use, copy, modify, and distribute this software for any
- * purpose without fee is hereby granted, provided that this entire notice
- * is included in all copies of any software which is or includes a copy
- * or modification of this software and in all copies of the supporting
- * documentation for such software.
- *
- * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
- * WARRANTY.  IN PARTICULAR, NEITHER THE AUTHOR NOR AT&T MAKES ANY
- * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
- * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
- *
- ***************************************************************/
-
-/* Please send bug reports to
-	David M. Gay
-	AT&T Bell Laboratories, Room 2C-463
-	600 Mountain Avenue
-	Murray Hill, NJ 07974-2070
-	U.S.A.
-	dmg@research.att.com or research!dmg
- */
+
+The author of this software is David M. Gay.
+
+Copyright (C) 1998-2001 by Lucent Technologies
+All Rights Reserved
+
+Permission to use, copy, modify, and distribute this software and
+its documentation for any purpose and without fee is hereby
+granted, provided that the above copyright notice appear in all
+copies and that both that the copyright notice and this
+permission notice and warranty disclaimer appear in supporting
+documentation, and that the name of Lucent or any of its entities
+not be used in advertising or publicity pertaining to
+distribution of the software without specific, written prior
+permission.
+
+LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
+INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
+IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
+SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
+IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
+ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
+THIS SOFTWARE.
+
+****************************************************************/
+
+/* Please send bug reports to David M. Gay (dmg at acm dot org,
+ * with " at " changed at "@" and " dot " changed to ".").	*/
+
+/* Original file gdtoa-strtod.c Modified 06-21-2006 by Jeff Johnston to work within newlib.  */
 
 #include <_ansi.h>
-#include <reent.h>
+#include <errno.h>
 #include <string.h>
 #include "mprec.h"
+#include "gdtoa.h"
+#include "gd_qnan.h"
+
+/* #ifndef NO_FENV_H */
+/* #include <fenv.h> */
+/* #endif */
+
+#ifdef USE_LOCALE
+#include "locale.h"
+#endif
+
+#ifdef IEEE_Arith
+#ifndef NO_IEEE_Scale
+#define Avoid_Underflow
+#undef tinytens
+/* The factor of 2^53 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
+		};
+#endif
+#endif
+
+#ifdef Honor_FLT_ROUNDS
+#define Rounding rounding
+#undef Check_FLT_ROUNDS
+#define Check_FLT_ROUNDS
+#else
+#define Rounding Flt_Rounds
+#endif
+
+static void
+_DEFUN (ULtod, (L, bits, exp, k),
+	__ULong *L _AND
+	__ULong *bits _AND
+	Long exp _AND
+	int k)
+{
+	switch(k & STRTOG_Retmask) {
+	  case STRTOG_NoNumber:
+	  case STRTOG_Zero:
+		L[0] = L[1] = 0;
+		break;
+
+	  case STRTOG_Denormal:
+		L[_1] = bits[0];
+		L[_0] = bits[1];
+		break;
+
+	  case STRTOG_Normal:
+	  case STRTOG_NaNbits:
+		L[_1] = bits[0];
+		L[_0] = (bits[1] & ~0x100000) | ((exp + 0x3ff + 52) << 20);
+		break;
+
+	  case STRTOG_Infinite:
+		L[_0] = 0x7ff00000;
+		L[_1] = 0;
+		break;
+
+	  case STRTOG_NaN:
+		L[_0] = 0x7fffffff;
+		L[_1] = (__ULong)-1;
+	  }
+	if (k & STRTOG_Neg)
+		L[_0] |= 0x80000000L;
+}
+
+#ifdef INFNAN_CHECK
+static int
+_DEFUN (match, (sp, t),
+	_CONST char **sp _AND
+	char *t)
+{
+	int c, d;
+	_CONST char *s = *sp;
+
+	while( (d = *t++) !=0) {
+		if ((c = *++s) >= 'A' && c <= 'Z')
+			c += 'a' - 'A';
+		if (c != d)
+			return 0;
+		}
+	*sp = s + 1;
+	return 1;
+}
+#endif /* INFNAN_CHECK */
+
 
 double
 _DEFUN (_strtod_r, (ptr, s00, se),
@@ -108,627 +201,919 @@ _DEFUN (_strtod_r, (ptr, s00, se),
 	_CONST char *s00 _AND
 	char **se)
 {
-  int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign, e1, esign, i, j,
-    k, nd, nd0, nf, nz, nz0, sign;
-  long e;
-  _CONST char *s, *s0, *s1, *s2;
-  double aadj, aadj1, adj;
-  long L;
-  unsigned long z;
-  __ULong y;
-  union double_union rv, rv0;
-  int nanflag;
-
-  _Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
-  sign = nz0 = nz = nanflag = 0;
-  rv.d = 0.;
-  for (s = s00;; s++)
-    switch (*s)
-      {
-      case '-':
-	sign = 1;
-	/* no break */
-      case '+':
-	if (*++s)
-	  goto break2;
-	/* no break */
-      case 0:
-	s = s00;
-	goto ret;
-      case '\t':
-      case '\n':
-      case '\v':
-      case '\f':
-      case '\r':
-      case ' ':
-	continue;
-      default:
-	goto break2;
-      }
-break2:
-  if (*s == 'n' || *s == 'N')
-    {
-      ++s; 
-      if (*s == 'a' || *s == 'A')
-        {
-          ++s;
-          if (*s == 'n' || *s == 'N')
-            {
-              nanflag = 1;
-              ++s;
-              goto ret;
-            }
-        }
-      s = s00;
-      goto ret;
-    }
-  else if (*s == '0')
-    {
-      nz0 = 1;
-      while (*++s == '0');
-      if (!*s)
-	goto ret;
-    }
-  s0 = s;
-  y = z = 0;
-  for (nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
-    if (nd < 9)
-      y = 10 * y + c - '0';
-    else if (nd < 16)
-      z = 10 * z + c - '0';
-  nd0 = nd;
-  if (c == '.')
-    {
-      c = *++s;
-      if (!nd)
-	{
-	  for (; c == '0'; c = *++s)
-	    nz++;
-	  if (c > '0' && c <= '9')
-	    {
-	      s0 = s;
-	      nf += nz;
-	      nz = 0;
-	      goto have_dig;
-	    }
-	  goto dig_done;
-	}
-      for (; c >= '0' && c <= '9'; c = *++s)
-	{
-	have_dig:
-	  nz++;
-	  if (c -= '0')
-	    {
-	      nf += nz;
-	      for (i = 1; i < nz; i++)
-		if (nd++ < 9)
-		  y *= 10;
-		else if (nd <= DBL_DIG + 1)
-		  z *= 10;
-	      if (nd++ < 9)
-		y = 10 * y + c;
-	      else if (nd <= DBL_DIG + 1)
-		z = 10 * z + c;
-	      nz = 0;
-	    }
-	}
-    }
-dig_done:
-  e = 0;
-  if (c == 'e' || c == 'E')
-    {
-      if (!nd && !nz && !nz0)
-	{
-	  s = s00;
-	  goto ret;
-	}
-      s2 = s;
-      esign = 0;
-      switch (c = *++s)
-	{
-	case '-':
-	  esign = 1;
-	case '+':
-	  c = *++s;
-	}
-      if (c >= '0' && c <= '9')
-	{
-	  while (c == '0')
-	    c = *++s;
-	  if (c > '0' && c <= '9')
-	    {
-	      e = c - '0';
-	      s1 = s;
-	      while ((c = *++s) >= '0' && c <= '9')
-		e = 10 * e + c - '0';
-	      if (s - s1 > 8)
-		/* Avoid confusion from exponents
-		 * so large that e might overflow.
-		 */
-		e = 9999999L;
-	      if (esign)
-		e = -e;
-	    }
-	  else
-	    e = 0;
-	}
-      else
-	s = s2;
-    }
-  if (!nd)
-    {
-      if (!nz && !nz0)
-	s = s00;
-      goto ret;
-    }
-  e1 = e -= nf;
-
-  /* Now we have nd0 digits, starting at s0, followed by a
-   * decimal point, followed by nd-nd0 digits.  The number we're
-   * after is the integer represented by those digits times
-   * 10**e */
-
-  if (!nd0)
-    nd0 = nd;
-  k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
-  rv.d = y;
-  if (k > 9)
-    rv.d = tens[k - 9] * rv.d + z;
-  bd0 = 0;
-  if (nd <= DBL_DIG
+#ifdef Avoid_Underflow
+	int scale;
+#endif
+	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, aadj1, adj, rv, rv0;
+	Long L;
+	__ULong y, z;
+	_Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
+#ifdef SET_INEXACT
+	int inexact, oldinexact;
+#endif
+#ifdef Honor_FLT_ROUNDS
+	int rounding;
+#endif
+
+	delta = bs = bd = NULL;
+	sign = nz0 = nz = decpt = 0;
+	dval(rv) = 0.;
+	for(s = s00;;s++) switch(*s) {
+		case '-':
+			sign = 1;
+			/* no break */
+		case '+':
+			if (*++s)
+				goto break2;
+			/* no break */
+		case 0:
+			goto ret0;
+		case '\t':
+		case '\n':
+		case '\v':
+		case '\f':
+		case '\r':
+		case ' ':
+			continue;
+		default:
+			goto break2;
+		}
+ break2:
+	if (*s == '0') {
+#ifndef NO_HEX_FP
+		{
+		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)
+			FPI fpi1 = fpi;
+			switch(fegetround()) {
+			  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
+			switch((i = gethex(ptr, &s, &fpi1, &exp, &bb, sign)) & STRTOG_Retmask) {
+			  case STRTOG_NoNumber:
+				s = s00;
+				sign = 0;
+			  case STRTOG_Zero:
+				break;
+			  default:
+				if (bb) {
+					copybits(bits, fpi.nbits, bb);
+					Bfree(ptr,bb);
+					}
+				ULtod(((U*)&rv)->L, bits, exp, i);
+			  }}
+			goto ret;
+		  }
+		}
+#endif
+		nz0 = 1;
+		while(*++s == '0') ;
+		if (!*s)
+			goto ret;
+		}
+	s0 = s;
+	y = z = 0;
+	for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
+		if (nd < 9)
+			y = 10*y + c - '0';
+		else if (nd < 16)
+			z = 10*z + c - '0';
+	nd0 = nd;
+#ifdef USE_LOCALE
+	if (c == *localeconv()->decimal_point)
+#else
+	if (c == '.')
+#endif
+		{
+		decpt = 1;
+		c = *++s;
+		if (!nd) {
+			for(; c == '0'; c = *++s)
+				nz++;
+			if (c > '0' && c <= '9') {
+				s0 = s;
+				nf += nz;
+				nz = 0;
+				goto have_dig;
+				}
+			goto dig_done;
+			}
+		for(; c >= '0' && c <= '9'; c = *++s) {
+ have_dig:
+			nz++;
+			if (c -= '0') {
+				nf += nz;
+				for(i = 1; i < nz; i++)
+					if (nd++ < 9)
+						y *= 10;
+					else if (nd <= DBL_DIG + 1)
+						z *= 10;
+				if (nd++ < 9)
+					y = 10*y + c;
+				else if (nd <= DBL_DIG + 1)
+					z = 10*z + c;
+				nz = 0;
+				}
+			}
+		}
+ dig_done:
+	e = 0;
+	if (c == 'e' || c == 'E') {
+		if (!nd && !nz && !nz0) {
+			goto ret0;
+			}
+		s00 = s;
+		esign = 0;
+		switch(c = *++s) {
+			case '-':
+				esign = 1;
+			case '+':
+				c = *++s;
+			}
+		if (c >= '0' && c <= '9') {
+			while(c == '0')
+				c = *++s;
+			if (c > '0' && c <= '9') {
+				L = c - '0';
+				s1 = s;
+				while((c = *++s) >= '0' && c <= '9')
+					L = 10*L + c - '0';
+				if (s - s1 > 8 || L > 19999)
+					/* Avoid confusion from exponents
+					 * so large that e might overflow.
+					 */
+					e = 19999; /* safe for 16 bit ints */
+				else
+					e = (int)L;
+				if (esign)
+					e = -e;
+				}
+			else
+				e = 0;
+			}
+		else
+			s = s00;
+		}
+	if (!nd) {
+		if (!nz && !nz0) {
+#ifdef INFNAN_CHECK
+			/* Check for Nan and Infinity */
+			__ULong bits[2];
+			static FPI fpinan =	/* only 52 explicit bits */
+				{ 52, 1-1023-53+1, 2046-1023-53+1, 1, SI };
+			if (!decpt)
+			 switch(c) {
+			  case 'i':
+			  case 'I':
+				if (match(&s,"nf")) {
+					--s;
+					if (!match(&s,"inity"))
+						++s;
+					dword0(rv) = 0x7ff00000;
+					dword1(rv) = 0;
+					goto ret;
+					}
+				break;
+			  case 'n':
+			  case 'N':
+				if (match(&s, "an")) {
+#ifndef No_Hex_NaN
+					if (*s == '(' /*)*/
+					 && hexnan(&s, &fpinan, bits)
+							== STRTOG_NaNbits) {
+						dword0(rv) = 0x7ff00000 | bits[1];
+						dword1(rv) = bits[0];
+						}
+					else {
+#endif
+						dword0(rv) = NAN_WORD0;
+						dword1(rv) = NAN_WORD1;
+#ifndef No_Hex_NaN
+						}
+#endif
+					goto ret;
+					}
+			  }
+#endif /* INFNAN_CHECK */
+ ret0:
+			s = s00;
+			sign = 0;
+			}
+		goto ret;
+		}
+	e1 = e -= nf;
+
+	/* Now we have nd0 digits, starting at s0, followed by a
+	 * decimal point, followed by nd-nd0 digits.  The number we're
+	 * after is the integer represented by those digits times
+	 * 10**e */
+
+	if (!nd0)
+		nd0 = nd;
+	k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
+	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;
+		}
+	bd0 = 0;
+	if (nd <= DBL_DIG
 #ifndef RND_PRODQUOT
-      && FLT_ROUNDS == 1
-#endif
-    )
-    {
-      if (!e)
-	goto ret;
-      if (e > 0)
-	{
-	  if (e <= Ten_pmax)
-	    {
+#ifndef Honor_FLT_ROUNDS
+		&& Flt_Rounds == 1
+#endif
+#endif
+			) {
+		if (!e)
+			goto ret;
+		if (e > 0) {
+			if (e <= Ten_pmax) {
 #ifdef VAX
-	      goto vax_ovfl_check;
+				goto vax_ovfl_check;
 #else
-	      /* rv.d = */ rounded_product (rv.d, tens[e]);
-	      goto ret;
-#endif
-	    }
-	  i = DBL_DIG - nd;
-	  if (e <= Ten_pmax + i)
-	    {
-	      /* A fancier test would sometimes let us do
+#ifdef Honor_FLT_ROUNDS
+				/* round correctly FLT_ROUNDS = 2 or 3 */
+				if (sign) {
+					rv = -rv;
+					sign = 0;
+					}
+#endif
+				/* rv = */ rounded_product(dval(rv), tens[e]);
+				goto ret;
+#endif
+				}
+			i = DBL_DIG - nd;
+			if (e <= Ten_pmax + i) {
+				/* A fancier test would sometimes let us do
 				 * this for larger i values.
 				 */
-	      e -= i;
-	      rv.d *= tens[i];
+#ifdef Honor_FLT_ROUNDS
+				/* round correctly FLT_ROUNDS = 2 or 3 */
+				if (sign) {
+					rv = -rv;
+					sign = 0;
+					}
+#endif
+				e -= 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;
-	      /* rv.d = */ rounded_product (rv.d, tens[e]);
-	      if ((word0 (rv) & Exp_mask)
-		  > Exp_msk1 * (DBL_MAX_EXP + Bias - 1 - P))
-		goto ovfl;
-	      word0 (rv) += P * Exp_msk1;
+				/* VAX exponent range is so narrow we must
+				 * worry about overflow here...
+				 */
+ vax_ovfl_check:
+				dword0(rv) -= P*Exp_msk1;
+				/* rv = */ rounded_product(dval(rv), tens[e]);
+				if ((dword0(rv) & Exp_mask)
+				 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
+					goto ovfl;
+				dword0(rv) += P*Exp_msk1;
 #else
-	      /* rv.d = */ rounded_product (rv.d, tens[e]);
+				/* rv = */ rounded_product(dval(rv), tens[e]);
 #endif
-	      goto ret;
-	    }
-	}
+				goto ret;
+				}
+			}
 #ifndef Inaccurate_Divide
-      else if (e >= -Ten_pmax)
-	{
-	  /* rv.d = */ rounded_quotient (rv.d, tens[-e]);
-	  goto ret;
-	}
-#endif
-    }
-  e1 += nd - k;
-
-  /* Get starting approximation = rv.d * 10**e1 */
-
-  if (e1 > 0)
-    {
-      if ((i = e1 & 15) != 0)
-	rv.d *= tens[i];
-      if (e1 &= ~15)
-	{
-	  if (e1 > DBL_MAX_10_EXP)
-	    {
-	    ovfl:
-	      ptr->_errno = ERANGE;
-#ifdef _HAVE_STDC
-	      rv.d = HUGE_VAL;
-#else
-	      /* Can't trust HUGE_VAL */
+		else if (e >= -Ten_pmax) {
+#ifdef Honor_FLT_ROUNDS
+			/* round correctly FLT_ROUNDS = 2 or 3 */
+			if (sign) {
+				rv = -rv;
+				sign = 0;
+				}
+#endif
+			/* rv = */ rounded_quotient(dval(rv), tens[-e]);
+			goto ret;
+			}
+#endif
+		}
+	e1 += nd - k;
+
 #ifdef IEEE_Arith
-	      word0 (rv) = Exp_mask;
-#ifndef _DOUBLE_IS_32BITS
-	      word1 (rv) = 0;
+#ifdef SET_INEXACT
+	inexact = 1;
+	if (k <= DBL_DIG)
+		oldinexact = get_inexact();
 #endif
-#else
-	      word0 (rv) = Big0;
-#ifndef _DOUBLE_IS_32BITS
-	      word1 (rv) = Big1;
-#endif
-#endif
-#endif
-	      if (bd0)
-		goto retfree;
-	      goto ret;
-	    }
-	  if (e1 >>= 4)
-	    {
-	      for (j = 0; e1 > 1; j++, e1 >>= 1)
-		if (e1 & 1)
-		  rv.d *= bigtens[j];
-	      /* The last multiplication could overflow. */
-	      word0 (rv) -= P * Exp_msk1;
-	      rv.d *= bigtens[j];
-	      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;
-#ifndef _DOUBLE_IS_32BITS
-		  word1 (rv) = Big1;
+#ifdef Avoid_Underflow
+	scale = 0;
 #endif
+#ifdef Honor_FLT_ROUNDS
+	if ((rounding = Flt_Rounds) >= 2) {
+		if (sign)
+			rounding = rounding == 2 ? 0 : 2;
+		else
+			if (rounding != 2)
+				rounding = 0;
 		}
-	      else
-		word0 (rv) += P * Exp_msk1;
-	    }
-
-	}
-    }
-  else if (e1 < 0)
-    {
-      e1 = -e1;
-      if ((i = e1 & 15) != 0)
-	rv.d /= tens[i];
-      if (e1 &= ~15)
-	{
-	  e1 >>= 4;
-	  if (e1 >= 1 << n_bigtens)
-            goto undfl;
-	  for (j = 0; e1 > 1; j++, e1 >>= 1)
-	    if (e1 & 1)
-	      rv.d *= tinytens[j];
-	  /* The last multiplication could underflow. */
-	  rv0.d = rv.d;
-	  rv.d *= tinytens[j];
-	  if (!rv.d)
-	    {
-	      rv.d = 2. * rv0.d;
-	      rv.d *= tinytens[j];
-	      if (!rv.d)
-		{
-		undfl:
-		  rv.d = 0.;
-		  ptr->_errno = ERANGE;
-		  if (bd0)
-		    goto retfree;
-		  goto ret;
+#endif
+#endif /*IEEE_Arith*/
+
+	/* Get starting approximation = rv * 10**e1 */
+
+	if (e1 > 0) {
+		if ( (i = e1 & 15) !=0)
+			dval(rv) *= tens[i];
+		if (e1 &= ~15) {
+			if (e1 > DBL_MAX_10_EXP) {
+ ovfl:
+#ifndef NO_ERRNO
+				ptr->_errno = ERANGE;
+#endif
+				/* Can't trust HUGE_VAL */
+#ifdef IEEE_Arith
+#ifdef Honor_FLT_ROUNDS
+				switch(rounding) {
+				  case 0: /* toward 0 */
+				  case 3: /* toward -infinity */
+					dword0(rv) = Big0;
+					dword1(rv) = Big1;
+					break;
+				  default:
+					dword0(rv) = Exp_mask;
+					dword1(rv) = 0;
+				  }
+#else /*Honor_FLT_ROUNDS*/
+				dword0(rv) = Exp_mask;
+				dword1(rv) = 0;
+#endif /*Honor_FLT_ROUNDS*/
+#ifdef SET_INEXACT
+				/* set overflow bit */
+				dval(rv0) = 1e300;
+				dval(rv0) *= dval(rv0);
+#endif
+#else /*IEEE_Arith*/
+				dword0(rv) = Big0;
+				dword1(rv) = Big1;
+#endif /*IEEE_Arith*/
+				if (bd0)
+					goto retfree;
+				goto ret;
+				}
+			e1 >>= 4;
+			for(j = 0; e1 > 1; j++, e1 >>= 1)
+				if (e1 & 1)
+					dval(rv) *= bigtens[j];
+		/* The last multiplication could overflow. */
+			dword0(rv) -= P*Exp_msk1;
+			dval(rv) *= bigtens[j];
+			if ((z = dword0(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) */
+				dword0(rv) = Big0;
+				dword1(rv) = Big1;
+				}
+			else
+				dword0(rv) += P*Exp_msk1;
+			}
 		}
-#ifndef _DOUBLE_IS_32BITS
-	      word0 (rv) = Tiny0;
-	      word1 (rv) = Tiny1;
-#else
-	      word0 (rv) = Tiny1;
-#endif
-	      /* The refinement below will clean
-	       * this approximation up.
-	       */
-	    }
-	}
-    }
-
-  /* Now the hard part -- adjusting rv to the correct value.*/
-
-  /* Put digits into bd: true value = bd * 10^e */
-
-  bd0 = s2b (ptr, s0, nd0, nd, y);
-
-  for (;;)
-    {
-      bd = Balloc (ptr, bd0->_k);
-      Bcopy (bd, bd0);
-      bb = d2b (ptr, rv.d, &bbe, &bbbits);	/* rv.d = bb * 2^bbe */
-      bs = i2b (ptr, 1);
-
-      if (e >= 0)
-	{
-	  bb2 = bb5 = 0;
-	  bd2 = bd5 = e;
-	}
-      else
-	{
-	  bb2 = bb5 = -e;
-	  bd2 = bd5 = 0;
-	}
-      if (bbe >= 0)
-	bb2 += bbe;
-      else
-	bd2 -= bbe;
-      bs2 = bb2;
-#ifdef Sudden_Underflow
-#ifdef IBM
-      j = 1 + 4 * P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
+	else if (e1 < 0) {
+		e1 = -e1;
+		if ( (i = e1 & 15) !=0)
+			dval(rv) /= tens[i];
+		if (e1 >>= 4) {
+			if (e1 >= 1 << n_bigtens)
+				goto undfl;
+#ifdef Avoid_Underflow
+			if (e1 & Scale_Bit)
+				scale = 2*P;
+			for(j = 0; e1 > 0; j++, e1 >>= 1)
+				if (e1 & 1)
+					dval(rv) *= tinytens[j];
+			if (scale && (j = 2*P + 1 - ((dword0(rv) & Exp_mask)
+						>> Exp_shift)) > 0) {
+				/* scaled rv is denormal; zap j low bits */
+				if (j >= 32) {
+					dword1(rv) = 0;
+					if (j >= 53)
+					 dword0(rv) = (P+2)*Exp_msk1;
+					else
+					 dword0(rv) &= 0xffffffff << (j-32);
+					}
+				else
+					dword1(rv) &= 0xffffffff << j;
+				}
 #else
-      j = P + 1 - bbbits;
+			for(j = 0; e1 > 1; j++, e1 >>= 1)
+				if (e1 & 1)
+					dval(rv) *= tinytens[j];
+			/* The last multiplication could underflow. */
+			dval(rv0) = dval(rv);
+			dval(rv) *= tinytens[j];
+			if (!dval(rv)) {
+				dval(rv) = 2.*dval(rv0);
+				dval(rv) *= tinytens[j];
 #endif
-#else
-      i = bbe + bbbits - 1;	/* logb(rv.d) */
-      if (i < Emin)		/* denormal */
-	j = bbe + (P - Emin);
-      else
-	j = P + 1 - bbbits;
-#endif
-      bb2 += j;
-      bd2 += j;
-      i = bb2 < bd2 ? bb2 : bd2;
-      if (i > bs2)
-	i = bs2;
-      if (i > 0)
-	{
-	  bb2 -= i;
-	  bd2 -= i;
-	  bs2 -= i;
-	}
-      if (bb5 > 0)
-	{
-	  bs = pow5mult (ptr, bs, bb5);
-	  bb1 = mult (ptr, bs, bb);
-	  Bfree (ptr, bb);
-	  bb = bb1;
-	}
-      if (bb2 > 0)
-	bb = lshift (ptr, bb, bb2);
-      if (bd5 > 0)
-	bd = pow5mult (ptr, bd, bd5);
-      if (bd2 > 0)
-	bd = lshift (ptr, bd, bd2);
-      if (bs2 > 0)
-	bs = lshift (ptr, bs, bs2);
-      delta = diff (ptr, bb, bd);
-      dsign = delta->_sign;
-      delta->_sign = 0;
-      i = cmp (delta, bs);
-      if (i < 0)
-	{
-	  /* 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)
-	    break;
-	  delta = lshift (ptr, delta, Log2P);
-	  if (cmp (delta, bs) > 0)
-	    goto drop_down;
-	  break;
-	}
-      if (i == 0)
-	{
-	  /* exactly half-way between */
-	  if (dsign)
-	    {
-	      if ((word0 (rv) & Bndry_mask1) == Bndry_mask1
-		  && word1 (rv) == 0xffffffff)
-		{
-		  /*boundary case -- increment exponent*/
-		  word0 (rv) = (word0 (rv) & Exp_mask)
-		    + Exp_msk1
-#ifdef IBM
-		    | Exp_msk1 >> 4
+				if (!dval(rv)) {
+ undfl:
+					dval(rv) = 0.;
+#ifndef NO_ERRNO
+					ptr->_errno = ERANGE;
 #endif
-		    ;
-#ifndef _DOUBLE_IS_32BITS
-		  word1 (rv) = 0;
+					if (bd0)
+						goto retfree;
+					goto ret;
+					}
+#ifndef Avoid_Underflow
+				dword0(rv) = Tiny0;
+				dword1(rv) = Tiny1;
+				/* The refinement below will clean
+				 * this approximation up.
+				 */
+				}
 #endif
-		  break;
+			}
 		}
-	    }
-	  else if (!(word0 (rv) & Bndry_mask) && !word1 (rv))
-	    {
-	    drop_down:
-	      /* boundary case -- decrement exponent */
+
+	/* Now the hard part -- adjusting rv to the correct value.*/
+
+	/* Put digits into bd: true value = bd * 10^e */
+
+	bd0 = s2b(ptr, s0, nd0, nd, y);
+
+	for(;;) {
+		bd = Balloc(ptr,bd0->_k);
+		Bcopy(bd, bd0);
+		bb = d2b(ptr,dval(rv), &bbe, &bbbits);	/* rv = bb * 2^bbe */
+		bs = i2b(ptr,1);
+
+		if (e >= 0) {
+			bb2 = bb5 = 0;
+			bd2 = bd5 = e;
+			}
+		else {
+			bb2 = bb5 = -e;
+			bd2 = bd5 = 0;
+			}
+		if (bbe >= 0)
+			bb2 += bbe;
+		else
+			bd2 -= bbe;
+		bs2 = bb2;
+#ifdef Honor_FLT_ROUNDS
+		if (rounding != 1)
+			bs2++;
+#endif
+#ifdef Avoid_Underflow
+		j = bbe - scale;
+		i = j + bbbits - 1;	/* logb(rv) */
+		if (i < Emin)	/* denormal */
+			j += P - Emin;
+		else
+			j = P + 1 - bbbits;
+#else /*Avoid_Underflow*/
 #ifdef Sudden_Underflow
-	      L = word0 (rv) & Exp_mask;
 #ifdef IBM
-	      if (L < Exp_msk1)
+		j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
 #else
-	      if (L <= Exp_msk1)
+		j = P + 1 - bbbits;
+#endif
+#else /*Sudden_Underflow*/
+		j = bbe;
+		i = j + bbbits - 1;	/* logb(rv) */
+		if (i < Emin)	/* denormal */
+			j += P - Emin;
+		else
+			j = P + 1 - bbbits;
+#endif /*Sudden_Underflow*/
+#endif /*Avoid_Underflow*/
+		bb2 += j;
+		bd2 += j;
+#ifdef Avoid_Underflow
+		bd2 += scale;
+#endif
+		i = bb2 < bd2 ? bb2 : bd2;
+		if (i > bs2)
+			i = bs2;
+		if (i > 0) {
+			bb2 -= i;
+			bd2 -= i;
+			bs2 -= i;
+			}
+		if (bb5 > 0) {
+			bs = pow5mult(ptr, bs, bb5);
+			bb1 = mult(ptr, bs, bb);
+			Bfree(ptr, bb);
+			bb = bb1;
+			}
+		if (bb2 > 0)
+			bb = lshift(ptr, bb, bb2);
+		if (bd5 > 0)
+			bd = pow5mult(ptr, bd, bd5);
+		if (bd2 > 0)
+			bd = lshift(ptr, bd, bd2);
+		if (bs2 > 0)
+			bs = lshift(ptr, bs, bs2);
+		delta = diff(ptr, bb, bd);
+		dsign = delta->_sign;
+		delta->_sign = 0;
+		i = cmp(delta, bs);
+#ifdef Honor_FLT_ROUNDS
+		if (rounding != 1) {
+			if (i < 0) {
+				/* Error is less than an ulp */
+				if (!delta->_x[0] && delta->_wds <= 1) {
+					/* exact */
+#ifdef SET_INEXACT
+					inexact = 0;
 #endif
-		goto undfl;
-	      L -= Exp_msk1;
+					break;
+					}
+				if (rounding) {
+					if (dsign) {
+						adj = 1.;
+						goto apply_adj;
+						}
+					}
+				else if (!dsign) {
+					adj = -1.;
+					if (!dword1(rv)
+					 && !(dword0(rv) & Frac_mask)) {
+						y = dword0(rv) & Exp_mask;
+#ifdef Avoid_Underflow
+						if (!scale || y > 2*P*Exp_msk1)
 #else
-	      L = (word0 (rv) & Exp_mask) - Exp_msk1;
+						if (y)
 #endif
-	      word0 (rv) = L | Bndry_mask1;
-#ifndef _DOUBLE_IS_32BITS
-	      word1 (rv) = 0xffffffff;
+						  {
+						  delta = lshift(ptr, delta,Log2P);
+						  if (cmp(delta, bs) <= 0)
+							adj = -0.5;
+						  }
+						}
+ apply_adj:
+#ifdef Avoid_Underflow
+					if (scale && (y = dword0(rv) & Exp_mask)
+						<= 2*P*Exp_msk1)
+					  dword0(adj) += (2*P+1)*Exp_msk1 - y;
+#else
+#ifdef Sudden_Underflow
+					if ((dword0(rv) & Exp_mask) <=
+							P*Exp_msk1) {
+						dword0(rv) += P*Exp_msk1;
+						dval(rv) += adj*ulp(dval(rv));
+						dword0(rv) -= P*Exp_msk1;
+						}
+					else
+#endif /*Sudden_Underflow*/
+#endif /*Avoid_Underflow*/
+					dval(rv) += adj*ulp(dval(rv));
+					}
+				break;
+				}
+			adj = ratio(delta, bs);
+			if (adj < 1.)
+				adj = 1.;
+			if (adj <= 0x7ffffffe) {
+				/* adj = rounding ? ceil(adj) : floor(adj); */
+				y = adj;
+				if (y != adj) {
+					if (!((rounding>>1) ^ dsign))
+						y++;
+					adj = y;
+					}
+				}
+#ifdef Avoid_Underflow
+			if (scale && (y = dword0(rv) & Exp_mask) <= 2*P*Exp_msk1)
+				dword0(adj) += (2*P+1)*Exp_msk1 - y;
+#else
+#ifdef Sudden_Underflow
+			if ((dword0(rv) & Exp_mask) <= P*Exp_msk1) {
+				dword0(rv) += P*Exp_msk1;
+				adj *= ulp(dval(rv));
+				if (dsign)
+					dval(rv) += adj;
+				else
+					dval(rv) -= adj;
+				dword0(rv) -= P*Exp_msk1;
+				goto cont;
+				}
+#endif /*Sudden_Underflow*/
+#endif /*Avoid_Underflow*/
+			adj *= ulp(dval(rv));
+			if (dsign)
+				dval(rv) += adj;
+			else
+				dval(rv) -= adj;
+			goto cont;
+			}
+#endif /*Honor_FLT_ROUNDS*/
+
+		if (i < 0) {
+			/* Error is less than half an ulp -- check for
+			 * special case of mantissa a power of two.
+			 */
+			if (dsign || dword1(rv) || dword0(rv) & Bndry_mask
+#ifdef IEEE_Arith
+#ifdef Avoid_Underflow
+			 || (dword0(rv) & Exp_mask) <= (2*P+1)*Exp_msk1
+#else
+			 || (dword0(rv) & Exp_mask) <= Exp_msk1
+#endif
+#endif
+				) {
+#ifdef SET_INEXACT
+				if (!delta->x[0] && delta->wds <= 1)
+					inexact = 0;
+#endif
+				break;
+				}
+			if (!delta->_x[0] && delta->_wds <= 1) {
+				/* exact result */
+#ifdef SET_INEXACT
+				inexact = 0;
+#endif
+				break;
+				}
+			delta = lshift(ptr,delta,Log2P);
+			if (cmp(delta, bs) > 0)
+				goto drop_down;
+			break;
+			}
+		if (i == 0) {
+			/* exactly half-way between */
+			if (dsign) {
+				if ((dword0(rv) & Bndry_mask1) == Bndry_mask1
+				 &&  dword1(rv) == (
+#ifdef Avoid_Underflow
+			(scale && (y = dword0(rv) & Exp_mask) <= 2*P*Exp_msk1)
+		? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
+#endif
+						   0xffffffff)) {
+					/*boundary case -- increment exponent*/
+					dword0(rv) = (dword0(rv) & Exp_mask)
+						+ Exp_msk1
+#ifdef IBM
+						| Exp_msk1 >> 4
+#endif
+						;
+					dword1(rv) = 0;
+#ifdef Avoid_Underflow
+					dsign = 0;
 #endif
+					break;
+					}
+				}
+			else if (!(dword0(rv) & Bndry_mask) && !dword1(rv)) {
+ drop_down:
+				/* boundary case -- decrement exponent */
+#ifdef Sudden_Underflow /*{{*/
+				L = dword0(rv) & Exp_mask;
+#ifdef IBM
+				if (L <  Exp_msk1)
+#else
+#ifdef Avoid_Underflow
+				if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1))
+#else
+				if (L <= Exp_msk1)
+#endif /*Avoid_Underflow*/
+#endif /*IBM*/
+					goto undfl;
+				L -= Exp_msk1;
+#else /*Sudden_Underflow}{*/
+#ifdef Avoid_Underflow
+				if (scale) {
+					L = dword0(rv) & Exp_mask;
+					if (L <= (2*P+1)*Exp_msk1) {
+						if (L > (P+2)*Exp_msk1)
+							/* round even ==> */
+							/* accept rv */
+							break;
+						/* rv = smallest denormal */
+						goto undfl;
+						}
+					}
+#endif /*Avoid_Underflow*/
+				L = (dword0(rv) & Exp_mask) - Exp_msk1;
+#endif /*Sudden_Underflow}*/
+				dword0(rv) = L | Bndry_mask1;
+				dword1(rv) = 0xffffffff;
 #ifdef IBM
-	      goto cont;
+				goto cont;
 #else
-	      break;
+				break;
 #endif
-	    }
+				}
 #ifndef ROUND_BIASED
-	  if (!(word1 (rv) & LSB))
-	    break;
+			if (!(dword1(rv) & LSB))
+				break;
 #endif
-	  if (dsign)
-	    rv.d += ulp (rv.d);
+			if (dsign)
+				dval(rv) += ulp(dval(rv));
 #ifndef ROUND_BIASED
-	  else
-	    {
-	      rv.d -= ulp (rv.d);
+			else {
+				dval(rv) -= ulp(dval(rv));
 #ifndef Sudden_Underflow
-	      if (!rv.d)
-		goto undfl;
-#endif
-	    }
-#endif
-	  break;
-	}
-      if ((aadj = ratio (delta, bs)) <= 2.)
-	{
-	  if (dsign)
-	    aadj = aadj1 = 1.;
-	  else if (word1 (rv) || word0 (rv) & Bndry_mask)
-	    {
+				if (!dval(rv))
+					goto undfl;
+#endif
+				}
+#ifdef Avoid_Underflow
+			dsign = 1 - dsign;
+#endif
+#endif
+			break;
+			}
+		if ((aadj = ratio(delta, bs)) <= 2.) {
+			if (dsign)
+				aadj = aadj1 = 1.;
+			else if (dword1(rv) || dword0(rv) & Bndry_mask) {
 #ifndef Sudden_Underflow
-	      if (word1 (rv) == Tiny1 && !word0 (rv))
-		goto undfl;
-#endif
-	      aadj = 1.;
-	      aadj1 = -1.;
-	    }
-	  else
-	    {
-	      /* special case -- power of FLT_RADIX to be */
-	      /* rounded down... */
-
-	      if (aadj < 2. / FLT_RADIX)
-		aadj = 1. / FLT_RADIX;
-	      else
-		aadj *= 0.5;
-	      aadj1 = -aadj;
-	    }
-	}
-      else
-	{
-	  aadj *= 0.5;
-	  aadj1 = dsign ? aadj : -aadj;
+				if (dword1(rv) == Tiny1 && !dword0(rv))
+					goto undfl;
+#endif
+				aadj = 1.;
+				aadj1 = -1.;
+				}
+			else {
+				/* special case -- power of FLT_RADIX to be */
+				/* rounded down... */
+
+				if (aadj < 2./FLT_RADIX)
+					aadj = 1./FLT_RADIX;
+				else
+					aadj *= 0.5;
+				aadj1 = -aadj;
+				}
+			}
+		else {
+			aadj *= 0.5;
+			aadj1 = dsign ? aadj : -aadj;
 #ifdef Check_FLT_ROUNDS
-	  switch (FLT_ROUNDS)
-	    {
-	    case 2:		/* towards +infinity */
-	      aadj1 -= 0.5;
-	      break;
-	    case 0:		/* towards 0 */
-	    case 3:		/* towards -infinity */
-	      aadj1 += 0.5;
-	    }
+			switch(Rounding) {
+				case 2: /* towards +infinity */
+					aadj1 -= 0.5;
+					break;
+				case 0: /* towards 0 */
+				case 3: /* towards -infinity */
+					aadj1 += 0.5;
+				}
 #else
-	  if (FLT_ROUNDS == 0)
-	    aadj1 += 0.5;
-#endif
-	}
-      y = word0 (rv) & Exp_mask;
-
-      /* Check for overflow */
-
-      if (y == Exp_msk1 * (DBL_MAX_EXP + Bias - 1))
-	{
-	  rv0.d = rv.d;
-	  word0 (rv) -= P * Exp_msk1;
-	  adj = aadj1 * ulp (rv.d);
-	  rv.d += adj;
-	  if ((word0 (rv) & Exp_mask) >=
-	      Exp_msk1 * (DBL_MAX_EXP + Bias - P))
-	    {
-	      if (word0 (rv0) == Big0 && word1 (rv0) == Big1)
-		goto ovfl;
-#ifdef _DOUBLE_IS_32BITS
-	      word0 (rv) = Big1;
+			if (Flt_Rounds == 0)
+				aadj1 += 0.5;
+#endif /*Check_FLT_ROUNDS*/
+			}
+		y = dword0(rv) & Exp_mask;
+
+		/* Check for overflow */
+
+		if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
+			dval(rv0) = dval(rv);
+			dword0(rv) -= P*Exp_msk1;
+			adj = aadj1 * ulp(dval(rv));
+			dval(rv) += adj;
+			if ((dword0(rv) & Exp_mask) >=
+					Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
+				if (dword0(rv0) == Big0 && dword1(rv0) == Big1)
+					goto ovfl;
+				dword0(rv) = Big0;
+				dword1(rv) = Big1;
+				goto cont;
+				}
+			else
+				dword0(rv) += P*Exp_msk1;
+			}
+		else {
+#ifdef Avoid_Underflow
+			if (scale && y <= 2*P*Exp_msk1) {
+				if (aadj <= 0x7fffffff) {
+					if ((z = aadj) <= 0)
+						z = 1;
+					aadj = z;
+					aadj1 = dsign ? aadj : -aadj;
+					}
+				dword0(aadj1) += (2*P+1)*Exp_msk1 - y;
+				}
+			adj = aadj1 * ulp(dval(rv));
+			dval(rv) += adj;
 #else
-	      word0 (rv) = Big0;
-	      word1 (rv) = Big1;
-#endif
-	      goto cont;
-	    }
-	  else
-	    word0 (rv) += P * Exp_msk1;
-	}
-      else
-	{
 #ifdef Sudden_Underflow
-	  if ((word0 (rv) & Exp_mask) <= P * Exp_msk1)
-	    {
-	      rv0.d = rv.d;
-	      word0 (rv) += P * Exp_msk1;
-	      adj = aadj1 * ulp (rv.d);
-	      rv.d += adj;
+			if ((dword0(rv) & Exp_mask) <= P*Exp_msk1) {
+				dval(rv0) = dval(rv);
+				dword0(rv) += P*Exp_msk1;
+				adj = aadj1 * ulp(dval(rv));
+				dval(rv) += adj;
 #ifdef IBM
-	      if ((word0 (rv) & Exp_mask) < P * Exp_msk1)
+				if ((dword0(rv) & Exp_mask) <  P*Exp_msk1)
 #else
-	      if ((word0 (rv) & Exp_mask) <= P * Exp_msk1)
+				if ((dword0(rv) & Exp_mask) <= P*Exp_msk1)
 #endif
-		{
-		  if (word0 (rv0) == Tiny0
-		      && word1 (rv0) == Tiny1)
-		    goto undfl;
-		  word0 (rv) = Tiny0;
-		  word1 (rv) = Tiny1;
-		  goto cont;
+					{
+					if (dword0(rv0) == Tiny0
+					 && dword1(rv0) == Tiny1)
+						goto undfl;
+					dword0(rv) = Tiny0;
+					dword1(rv) = Tiny1;
+					goto cont;
+					}
+				else
+					dword0(rv) -= P*Exp_msk1;
+				}
+			else {
+				adj = aadj1 * ulp(dval(rv));
+				dval(rv) += adj;
+				}
+#else /*Sudden_Underflow*/
+			/* Compute adj so that the IEEE rounding rules will
+			 * correctly round rv + adj in some half-way cases.
+			 * 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);
+				if (!dsign)
+					aadj1 = -aadj1;
+				}
+			adj = aadj1 * ulp(dval(rv));
+			dval(rv) += adj;
+#endif /*Sudden_Underflow*/
+#endif /*Avoid_Underflow*/
+			}
+		z = dword0(rv) & Exp_mask;
+#ifndef SET_INEXACT
+#ifdef Avoid_Underflow
+		if (!scale)
+#endif
+		if (y == z) {
+			/* Can we stop now? */
+			L = (Long)aadj;
+			aadj -= L;
+			/* The tolerances below are conservative. */
+			if (dsign || dword1(rv) || dword0(rv) & Bndry_mask) {
+				if (aadj < .4999999 || aadj > .5000001)
+					break;
+				}
+			else if (aadj < .4999999/FLT_RADIX)
+				break;
+			}
+#endif
+ cont:
+		Bfree(ptr,bb);
+		Bfree(ptr,bd);
+		Bfree(ptr,bs);
+		Bfree(ptr,delta);
 		}
-	      else
-		word0 (rv) -= P * Exp_msk1;
-	    }
-	  else
-	    {
-	      adj = aadj1 * ulp (rv.d);
-	      rv.d += adj;
-	    }
-#else
-	  /* Compute adj so that the IEEE rounding rules will
-	   * correctly round rv.d + adj in some half-way cases.
-	   * If rv.d * ulp(rv.d) 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);
-	      if (!dsign)
-		aadj1 = -aadj1;
-	    }
-	  adj = aadj1 * ulp (rv.d);
-	  rv.d += adj;
-#endif
-	}
-      z = word0 (rv) & Exp_mask;
-      if (y == z)
-	{
-	  /* Can we stop now? */
-	  L = aadj;
-	  aadj -= L;
-	  /* The tolerances below are conservative. */
-	  if (dsign || word1 (rv) || word0 (rv) & Bndry_mask)
-	    {
-	      if (aadj < .4999999 || aadj > .5000001)
-		break;
-	    }
-	  else if (aadj < .4999999 / FLT_RADIX)
-	    break;
-	}
-    cont:
-      Bfree (ptr, bb);
-      Bfree (ptr, bd);
-      Bfree (ptr, bs);
-      Bfree (ptr, delta);
-    }
-retfree:
-  Bfree (ptr, bb);
-  Bfree (ptr, bd);
-  Bfree (ptr, bs);
-  Bfree (ptr, bd0);
-  Bfree (ptr, delta);
-ret:
-  if (se)
-    *se = (char *) s;
-
-  if (nanflag)
-    return nan (NULL);
-  return (sign && (s != s00)) ? -rv.d : rv.d;
+#ifdef SET_INEXACT
+	if (inexact) {
+		if (!oldinexact) {
+			dword0(rv0) = Exp_1 + (70 << Exp_shift);
+			dword1(rv0) = 0;
+			dval(rv0) += 1.;
+			}
+		}
+	else if (!oldinexact)
+		clear_inexact();
+#endif
+#ifdef Avoid_Underflow
+	if (scale) {
+		dword0(rv0) = Exp_1 - 2*P*Exp_msk1;
+		dword1(rv0) = 0;
+		dval(rv) *= dval(rv0);
+#ifndef NO_ERRNO
+		/* try to avoid the bug of testing an 8087 register value */
+		if (dword0(rv) == 0 && dword1(rv) == 0)
+			ptr->_errno = ERANGE;
+#endif
+		}
+#endif /* Avoid_Underflow */
+#ifdef SET_INEXACT
+	if (inexact && !(dword0(rv) & Exp_mask)) {
+		/* set underflow bit */
+		dval(rv0) = 1e-300;
+		dval(rv0) *= dval(rv0);
+		}
+#endif
+ retfree:
+	Bfree(ptr,bb);
+	Bfree(ptr,bd);
+	Bfree(ptr,bs);
+	Bfree(ptr,bd0);
+	Bfree(ptr,delta);
+ ret:
+	if (se)
+		*se = (char *)s;
+	return sign ? -dval(rv) : dval(rv);
 }
 
 #ifndef NO_REENT