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.

1 files changed, 991 insertions, 606 deletions

diff --git a/newlib/libc/stdlib/strtod.c b/newlib/libc/stdlib/strtod.c
index 9b70dfc..57297e0 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