TIP#237 IMPLEMENTATION

        [kennykb-numerics-branch] Resynchronized with the HEAD; at this
         checkpoint [-rkennykb-numerics-branch-20051008], the HEAD and
        kennykb-numerics-branch contain identical code.

1 files changed, 1919 insertions, 643 deletions

diff --git a/generic/tclStrToD.c b/generic/tclStrToD.c
index 68f6c30..7f8ecd5 100755
--- a/generic/tclStrToD.c
+++ b/generic/tclStrToD.c
@@ -1,21 +1,20 @@
 /*
  *----------------------------------------------------------------------
  *
- * tclStrToD.c --
+ * tclDouble.c --
  *
- *	This file contains a TclStrToD procedure that handles conversion of
- *	string to double, with correct rounding even where extended precision
- *	is needed to achieve that.  It also contains a TclDoubleDigits
- *	procedure that handles conversion of double to string (at least the
- *	significand), and several utility functions for interconverting
- *	'double' and the integer types.
+ *	This file contains a collection of procedures for managing
+ *	conversions to/from floating-point in Tcl.  They include
+ *	TclParseNumber, which parses numbers from strings; TclDoubleDigits,
+ *	which formats numbers into strings of digits, and procedures for
+ *	interconversion among 'double' and 'mp_int' types.
  *
  * Copyright (c) 2005 by Kevin B. Kenny.  All rights reserved.
  *
  * See the file "license.terms" for information on usage and redistribution
  * of this file, and for a DISCLAIMER OF ALL WARRANTIES.
  *
- * RCS: @(#) $Id: tclStrToD.c,v 1.8 2005/08/24 15:15:45 kennykb Exp $
+ * RCS: @(#) $Id: tclStrToD.c,v 1.9 2005/10/08 14:42:45 dgp Exp $
  *
  *----------------------------------------------------------------------
  */
@@ -30,34 +29,42 @@
 #include <tommath.h>
 
 /*
- * The stuff below is a bit of a hack so that this file can be used in
- * environments that include no UNIX, i.e. no errno: just arrange to use the
- * errno from tclExecute.c here.
+ * Define TIP_114_FORMATS to accept 0b and 0o for binary and octal strings.
+ * Define KILL_OCTAL as well as TIP_114_FORMATS to suppress interpretation
+ * of numbers with leading zero as octal. (Ceterum censeo: numeros octonarios
+ * delendos esse.)
  */
 
-#ifdef TCL_GENERIC_ONLY
-#   define NO_ERRNO_H
-#endif
+#define	TIP_114_FORMATS
+#undef	KILL_OCTAL
 
-#ifdef NO_ERRNO_H
-extern int errno;			/* Use errno from tclExecute.c. */
-#   define ERANGE 34
+#ifndef TIP_114_FORMATS
+#undef	KILL_OCTAL
 #endif
 
+/*
+ * This code supports (at least hypothetically), IBM, Cray, VAX and
+ * IEEE-754 floating point; of these, only IEEE-754 can represent NaN.
+ * IEEE-754 can be uniquely determined by radix and by the widths of
+ * significand and exponent.
+ */
+
 #if (FLT_RADIX == 2) && (DBL_MANT_DIG == 53) && (DBL_MAX_EXP == 1024)
 #   define IEEE_FLOATING_POINT
 #endif
 
 /*
- * gcc on x86 needs access to rounding controls.  It is tempting to include
- * fpu_control.h, but that file exists only on Linux; it is missing on Cygwin
- * and MinGW. Most gcc-isms and ix86-isms are factored out here.
+ * gcc on x86 needs access to rounding controls, because of a questionable
+ * feature where it retains intermediate results as IEEE 'long double' values
+ * somewhat unpredictably.  It is tempting to include fpu_control.h, but
+ * that file exists only on Linux; it is missing on Cygwin and MinGW. Most 
+ * gcc-isms and ix86-isms are factored out here.
  */
 
 #if defined(__GNUC__) && defined(__i386)
 typedef unsigned int fpu_control_t __attribute__ ((__mode__ (__HI__)));
-#   define _FPU_GETCW(cw)	__asm__ ("fnstcw %0" : "=m" (*&cw))
-#   define _FPU_SETCW(cw)	__asm__ ("fldcw %0" : : "m" (*&cw))
+#define _FPU_GETCW(cw) __asm__ __volatile__ ("fnstcw %0" : "=m" (*&cw))
+#define _FPU_SETCW(cw) __asm__ __volatile__ ("fldcw %0" : : "m" (*&cw))
 #   define FPU_IEEE_ROUNDING	0x027f
 #   define ADJUST_FPU_CONTROL_WORD
 #endif
@@ -75,13 +82,26 @@ typedef unsigned int fpu_control_t __attribute__ ((__mode__ (__HI__)));
 #   define NAN_MASK (((Tcl_WideUInt) 1) << 51)
 #endif
 
-/*
- * There now follows a lot of static variables that are shared across all
- * threads but which are not guarded by mutexes. This is OK, because they are
- * only ever assigned _ONCE_ during Tcl's library initialization sequence.
- */
+/* The powers of ten that can be represented exactly as wide integers */
+
+static int maxpow10_wide;
+static Tcl_WideUInt *pow10_wide;
+
+/* The number of decimal digits that fit in an mp_digit */
+
+static int log10_DIGIT_MAX;
+
+/* The powers of ten that can be represented exactly as IEEE754 doubles. */
+
+#define MAXPOW 22
+static double pow10 [MAXPOW+1];
+
+static int mmaxpow;		/* Largest power of ten that can be
+				 * represented exactly in a 'double'. */
 
-static const double pow_10_2_n[] = {	/* Inexact higher powers of ten */
+/* Inexact higher powers of ten */
+
+static CONST double pow_10_2_n [] = {
     1.0,
     100.0,
     10000.0,
@@ -92,440 +112,1414 @@ static const double pow_10_2_n[] = {	/* Inexact higher powers of ten */
     1.0e+128,
     1.0e+256
 };
-#define MAXPOW	22		/* Num of exactly representable powers of 10 */
-static double pow10[MAXPOW+1];	/* The powers of ten that can be represented
-				 * exactly as IEEE754 doubles. */
-static int mmaxpow;		/* Largest power of ten that can be
-				 * represented exactly in a 'double'. */
-static int log2FLT_RADIX;	/* Logarithm of the floating point radix. */
-static int mantBits;		/* Number of bits in a double's significand. */
-static mp_int pow5[9];		/* Table of powers of 5**(2**n), up to
-				 * 5**256. */
-static double tiny;		/* The smallest representable double. */
-static int maxDigits;		/* The maximum number of digits to the left of
-				 * the decimal point of a double. */
-static int minDigits;		/* The maximum number of digits to the right
-				 * of the decimal point in a double. */
-static int mantDIGIT;		/* Number of mp_digit's needed to hold the
-				 * significand of a double. */
+
+ /* Logarithm of the floating point radix. */
+
+static int log2FLT_RADIX;
+
+/* Number of bits in a double's significand */
+
+static int mantBits;
+
+/* Table of powers of 5**(2**n), up to 5**256 */
+
+static mp_int pow5[9];
+
+/* The smallest representable double */
+
+static double tiny;
+
+/* The maximum number of digits to the left of the decimal point of a
+ * double. */
+
+static int maxDigits;
+
+/* The maximum number of digits to the right of the decimal point in a
+ * double. */
+
+static int minDigits;
+
+/* Number of mp_digit's needed to hold the significand of a double */
+
+static int mantDIGIT;
 
 /* Static functions defined in this file */
 
-static double	RefineResult(double approx, CONST char *start, int nDigits,
-		    long exponent);
-static double	ParseNaN(int signum, CONST char **end);
-static double	SafeLdExp(double fraction, int exponent);
+static int AccumulateDecimalDigit _ANSI_ARGS_((unsigned, int, 
+					       Tcl_WideUInt*, mp_int*, int));
+static double MakeLowPrecisionDouble _ANSI_ARGS_((int signum,
+						  Tcl_WideUInt significand,
+						  int nSigDigs,
+						  int exponent));
+static double MakeHighPrecisionDouble _ANSI_ARGS_((int signum,
+						   mp_int* significand,
+						   int nSigDigs,
+						   int exponent));
+static double MakeNaN _ANSI_ARGS_(( int signum, Tcl_WideUInt tag ));
+static double RefineApproximation _ANSI_ARGS_((double approx,
+					       mp_int* exactSignificand,
+					       int exponent));
+static double AbsoluteValue(double v, int* signum);
+static int GetIntegerTimesPower(double v, mp_int* r, int* e);
+static double BignumToBiasedFrExp _ANSI_ARGS_(( mp_int* big, int* machexp ));
+static double Pow10TimesFrExp _ANSI_ARGS_(( int exponent,
+					    double fraction,
+					    int* machexp ));
+static double SafeLdExp _ANSI_ARGS_(( double fraction, int exponent ));
+
 
 /*
  *----------------------------------------------------------------------
  *
- * TclStrToD --
+ * TclParseNumber --
  *
- *	Scans a double from a string.
+ *	Place a "numeric" internal representation on a Tcl object.
  *
  * Results:
- *	Returns the scanned number. In the case of underflow, returns an
- *	appropriately signed zero; in the case of overflow, returns an
- *	appropriately signed HUGE_VAL.
+ *	Returns a standard Tcl result.
  *
  * Side effects:
- *	Stores a pointer to the end of the scanned number in '*endPtr', if
- *	endPtr is not NULL. If '*endPtr' is equal to 's' on return from this
- *	function, it indicates that the input string could not be recognized
- *	as a number.  In the case of underflow or overflow, 'errno' is set to
- *	ERANGE.
+ *	Stores an internal representation appropriate to the string.
+ *	The internal representation may be an integer, a wide integer,
+ *	a bignum, or a double.
+ *
+ * TclMakeObjNumeric is called as a common scanner in routines
+ * that expect numbers in Tcl_Obj's.  It scans the string representation
+ * of a given Tcl_Obj and stores an internal rep that represents
+ * a "canonical" version of its numeric value.  The value of the
+ * canonicalization is that a routine can determine simply by
+ * examining the type pointer whether an object LooksLikeInt,
+ * what size of integer is needed to hold it, and similar questions,
+ * and never needs to refer back to the string representation, even
+ * for "impure" objects.
+ *
+ * The 'strPtr' and 'endPtrPtr' arguments allow for recognizing a number
+ * that is in a substring of a Tcl_Obj, for example a screen metric or
+ * "end-" index.  If 'strPtr' is not NULL, it designates where the
+ * number begins within the string.  (The default is the start of
+ * objPtr's string rep, which will be constructed if necessary.)
  *
- *------------------------------------------------------------------------
+ * If 'strPtr' is supplied, 'objPtr' may be NULL.  In this case,
+ * no internal representation will be generated; instead, the routine
+ * will simply check for a syntactically correct number, returning
+ * TCL_OK or TCL_ERROR as appropriate, and setting *endPtrPtr if
+ * necessary.
+ *
+ * If 'endPtrPtr' is not NULL, it designates the first character
+ * after the scanned number.  In this case, successfully recognizing
+ * any digits will yield a return code of TCL_OK.  Only in the case
+ * where no leading string of 'strPtr' (or of objPtr's internal rep)
+ * represents a number will TCL_ERROR be returned.
+ *
+ * When only a partial string is being recognized, it is the caller's
+ * responsibility to destroy the internal representation, or at
+ * least change its type.  Failure to do so will lead to subsequent
+ * problems where a string that does not represent a number will
+ * be recognized as one because it has a numeric internal representation.
+ *
+ * When the 'flags' word includes TCL_PARSE_DECIMAL_ONLY, only decimal
+ * numbers are recognized; leading 0 has no special interpretation as
+ * octal and leading '0x' is forbidden.
+ *
+ *----------------------------------------------------------------------
  */
 
-double
-TclStrToD(CONST char *s,	/* String to scan. */
-	  CONST char **endPtr)	/* Pointer to the end of the scanned number. */
+int
+TclParseNumber( Tcl_Interp* interp,
+				/* Tcl interpreter for error reporting.
+				 * May be NULL */
+		Tcl_Obj* objPtr,
+				/* Object to receive the internal rep */
+		CONST char* type,
+				/* Type of number being parsed ("integer",
+				 * "wide integer", etc. */
+		CONST char* string,
+				/* Pointer to the start of the string to
+				 * scan, see above */
+		size_t length,	/* Maximum length of the string to scan,
+				 * see above. */
+		CONST char** endPtrPtr,
+				/* (Output) pointer to the end of the
+				 * scanned number, see above */
+		int flags)	/* Flags governing the parse */
 {
-    const char *p = s;
-    const char *startOfSignificand = NULL;
-				/* Start of the significand in the string. */
-    int signum = 0;		/* Sign of the significand. */
-    double exactSignificand = 0.0;
-				/* Significand, represented exactly as a
-				 * floating-point number. */
-    int seenDigit = 0;		/* Flag == 1 if a digit has been seen. */
-    int nSigDigs = 0;		/* Number of significant digits presented. */
-    int nDigitsAfterDp = 0;	/* Number of digits after the decimal point. */
-    int nTrailZero = 0;		/* Number of trailing zeros in the
-				 * significand. */
-    long exponent = 0;		/* Exponent. */
-    int seenDp = 0;		/* Flag == 1 if decimal point has been seen. */
-    char c;			/* One character extracted from the input. */
-    volatile double v;		/* Scanned value; must be 'volatile double' on
-				 * gc-ix86 to force correct rounding to IEEE
-				 * double and not Intel double-extended. */
-    int machexp;		/* Exponent of the machine rep of the scanned
-				 * value. */
-    int expt2;			/* Exponent for computing first approximation
-				 * to the true value. */
-    int i, j;
 
-    /*
-     * With gcc on x86, the floating point rounding mode is double-extended.
-     * This causes the result of double-precision calculations to be rounded
-     * twice: once to the precision of double-extended and then again to the
-     * precision of double.  Double-rounding introduces gratuitous errors of
-     * one ulp, so we need to change rounding mode to 53-bits.
-     */
-
-#ifdef ADJUST_FPU_CONTROL_WORD
-    fpu_control_t roundTo53Bits = FPU_IEEE_ROUNDING;
-    fpu_control_t oldRoundingMode;
-    _FPU_GETCW(oldRoundingMode);
-    _FPU_SETCW(roundTo53Bits);
-#   define RestoreRoundingMode()	_FPU_SETCW(oldRoundingMode)
-#else
-#   define RestoreRoundingMode()	(void) 0 /* Do nothing */
+    enum State {
+	INITIAL, SIGNUM, ZERO, ZERO_X, 
+#ifdef TIP_114_FORMATS
+	ZERO_O, ZERO_B, BINARY,
+#endif
+	HEXADECIMAL, OCTAL, BAD_OCTAL, DECIMAL, 
+	LEADING_RADIX_POINT, FRACTION,
+	EXPONENT_START, EXPONENT_SIGNUM, EXPONENT,
+	sI, sIN, sINF, sINFI, sINFIN, sINFINI, sINFINIT, sINFINITY
+#ifdef IEEE_FLOATING_POINT
+	, sN, sNA, sNAN, sNANPAREN, sNANHEX, sNANFINISH
+#endif
+    } state = INITIAL;
+    enum State acceptState = INITIAL;
+
+    int signum = 0;		/* Sign of the number being parsed */
+    Tcl_WideUInt significandWide = 0;
+				/* Significand of the number being
+				 * parsed (if no overflow) */
+    mp_int significandBig;	/* Significand of the number being
+				 * parsed (if it overflows significandWide) */
+    int significandOverflow = 0;
+				/* Flag==1 iff significandBig is used */
+    Tcl_WideUInt octalSignificandWide = 0;
+				/* Significand of an octal number; needed
+				 * because we don't know whether a number
+				 * with a leading zero is octal or decimal
+				 * until we've scanned forward to a '.' or
+				 * 'e' */
+    mp_int octalSignificandBig;	/* Significand of octal number once
+				 * octalSignificandWide overflows */
+    int octalSignificandOverflow = 0;
+				/* Flag==1 if octalSignificandBig is used */
+    int numSigDigs = 0;		/* Number of significant digits in the
+				 * decimal significand */
+    int numTrailZeros = 0;	/* Number of trailing zeroes at the
+				 * current point in the parse. */
+    int numDigitsAfterDp = 0;	/* Number of digits scanned after the
+				 * decimal point */
+    int exponentSignum = 0;	/* Signum of the exponent of a floating
+				 * point number */
+    long exponent = 0;		/* Exponent of a floating point number */
+    CONST char* p;		/* Pointer to next character to scan */
+    size_t len;			/* Number of characters remaining after p */
+    CONST char* acceptPoint;	/* Pointer to position after last character
+				 * in an acceptable number */
+    size_t acceptLen;		/* Number of characters following that point */
+    int status = TCL_OK;	/* Status to return to caller */
+    char d;			/* Last hexadecimal digit scanned */
+    int shift = 0;		/* Amount to shift when accumulating binary */
+#ifdef TIP_114_FORMATS
+    int explicitOctal = 0;
 #endif
 
-    /*
-     * Discard leading whitespace from input.
+    /* 
+     * Initialize string to start of the object's string rep if
+     * the caller didn't pass anything else.
      */
 
-    while (isspace(UCHAR(*p))) {
-	++p;
+    if ( string == NULL ) {
+	string = Tcl_GetStringFromObj( objPtr, NULL );
     }
 
-    /*
-     * Determine the sign of the significand.
-     */
+    p = string;
+    len = length;
+    acceptPoint = p;
+    acceptLen = len;
+    while ( 1 ) {
+	char c = len ? *p : '\0';
+	switch (state) {
 
-    switch (*p) {
-	case '-':
-	    signum = 1;
+	case INITIAL:
+	    /*
+	     * Initial state. Acceptable characters are +, -, digits,
+	     * period, I, N, and whitespace.
+	     */
+	    if (isspace(UCHAR(c))) {
+		break;
+	    } else if (c == '+') {
+		state = SIGNUM;
+		break;
+	    } else if (c == '-') {
+		signum = 1;
+		state = SIGNUM;
+		break;
+	    } 
 	    /* FALLTHROUGH */
-	case '+':
-	    ++p;
-    }
-
-    /*
-     * Discard leading zeroes from input.
-     */
+	    
+	case SIGNUM:
+	    /*
+	     * Scanned a leading + or -. Acceptable characters are
+	     * digits, period, I, and N.
+	     */
+	    if (c == '0') {
+		if (flags & TCL_PARSE_DECIMAL_ONLY) {
+		    state = DECIMAL;
+		} else {
+		    state = ZERO;
+		}
+		break;
+	    } else if (flags & TCL_PARSE_HEXADECIMAL_ONLY) {
+		goto zerox;
+	    } else if (flags & TCL_PARSE_OCTAL_ONLY) {
+		goto zeroo;
+	    } else if (isdigit(UCHAR(c))) {
+		significandWide = c - '0';
+		numSigDigs = 1;
+		state = DECIMAL;
+		break;
+	    } else if (flags & TCL_PARSE_INTEGER_ONLY) {
+		goto endgame;
+	    } else if (c == '.') {
+		state = LEADING_RADIX_POINT;
+		break;
+	    } else if (c == 'I' || c == 'i') {
+		state = sI;
+		break;
+#ifdef IEEE_FLOATING_POINT
+	    } else if (c == 'N' || c == 'n') {
+		state = sN;
+		break;
+#endif
+	    }
+	    goto endgame;
 
-    while (*p == '0') {
-	seenDigit = 1;
-	++p;
-    }
+	case ZERO:
+	    /*
+	     * Scanned a leading zero (perhaps with a + or -).
+	     * Acceptable inputs are digits, period, X, and E.
+	     * If 8 or 9 is encountered, the number can't be
+	     * octal.  This state and the OCTAL state differ only
+	     * in whether they recognize 'X'.
+	     */
+	    acceptState = state;
+	    acceptPoint = p;
+	    acceptLen = len;
+	    if (c == 'x' || c == 'X') {
+		state = ZERO_X;
+		break;
+	    }
+	    if (flags & TCL_PARSE_HEXADECIMAL_ONLY) {
+		goto zerox;
+	    }
+#ifdef TIP_114_FORMATS 
+	    if (flags & TCL_PARSE_SCAN_PREFIXES) {
+		goto zeroo;
+	    }
+	    if (c == 'b' || c == 'B') {
+		state = ZERO_B;
+		break;
+	    }
+	    if (c == 'o' || c == 'O') {
+		explicitOctal = 1;
+		state = ZERO_O;
+		break;
+	    }
+#ifdef KILL_OCTAL
+	    goto decimal;
+#endif
+#endif
+	    /* FALLTHROUGH */
 
-    /*
-     * Scan digits from the significand. Simultaneously, keep track of the
-     * number of digits after the decimal point. Maintain a pointer to the
-     * start of the significand. Keep "exactSignificand" equal to the
-     * conversion of the DBL_DIG most significant digits.
-     */
+	case OCTAL:
+	    /*
+	     * Scanned an optional + or -, followed by a string of
+	     * octal digits.  Acceptable inputs are more digits,
+	     * period, or E.  If 8 or 9 is encountered, commit to
+	     * floating point.
+	     */
+	    acceptState = state;
+	    acceptPoint = p;
+	    acceptLen = len;
+#ifdef TIP_114_FORMATS
+	    /* FALLTHROUGH */
+	case ZERO_O:
+#endif
+	zeroo:
+	    if (c == '0') {
+		++numTrailZeros;
+		state = OCTAL;
+		break;
+	    } else if (c >= '1' && c <= '7') {
+		if (objPtr != NULL) {
+		    shift = 3 * (numTrailZeros + 1);
+		    significandOverflow = 
+			AccumulateDecimalDigit((unsigned)(c-'0'),
+					       numTrailZeros, 
+					       &significandWide,
+					       &significandBig, 
+					       significandOverflow);
+		    
+		    if (!octalSignificandOverflow) {
+			/*
+			 * Shifting by more bits than are in the value being
+			 * shifted is at least de facto nonportable.  Check
+			 * for too large shifts first.
+			 */
+			if ((octalSignificandWide != 0) 
+				&& ((shift >= CHAR_BIT*sizeof(Tcl_WideUInt)) 
+				|| (octalSignificandWide 
+				> (~(Tcl_WideUInt)0 >> shift)))) {
+			    octalSignificandOverflow = 1;
+			    TclBNInitBignumFromWideUInt(&octalSignificandBig,
+							octalSignificandWide);
+			}
+		    }
+		    if (!octalSignificandOverflow) {
+			octalSignificandWide
+			    = (octalSignificandWide << shift) + (c - '0');
+		    } else {
+			mp_mul_2d(&octalSignificandBig, shift,
+				  &octalSignificandBig);
+			mp_add_d(&octalSignificandBig, (mp_digit)(c - '0'),
+				 &octalSignificandBig);
+		    }
+		}
+		if ( numSigDigs != 0 ) {
+		    numSigDigs += ( numTrailZeros + 1 );
+		} else {
+		    numSigDigs = 1;
+		}
+		numTrailZeros = 0;
+		state = OCTAL;
+		break;
+	    }
+	    /* FALLTHROUGH */
 
-    for (;;) {
-	c = *p;
-	if (c == '.' && !seenDp) {
-	    seenDp = 1;
-	    ++p;
-	} else if (isdigit(UCHAR(c))) {
+	case BAD_OCTAL:
+#ifdef TIP_114_FORMATS
+	    if (explicitOctal) {
+		/* No forgiveness for bad digits in explicitly octal numbers */
+		goto endgame;
+	    }
+#endif
+	    if (flags & TCL_PARSE_INTEGER_ONLY) {
+		/* No seeking floating point when parsing only integer */
+		goto endgame;
+	    }
+#ifndef KILL_OCTAL
+	    /*
+	     * Scanned a number with a leading zero that contains an
+	     * 8, 9, radix point or E.  This is an invalid octal number,
+	     * but might still be floating point.  
+	     */
 	    if (c == '0') {
-		if (startOfSignificand != NULL) {
-		    ++nTrailZero;
+		++numTrailZeros;
+		state = BAD_OCTAL;
+		break;
+	    } else if (isdigit(UCHAR(c))) {
+		if (objPtr != NULL) {
+		    significandOverflow = 
+			AccumulateDecimalDigit((unsigned)(c-'0'),
+					       numTrailZeros, 
+					       &significandWide,
+					       &significandBig, 
+					       significandOverflow);
+		}
+		if ( numSigDigs != 0 ) {
+		    numSigDigs += ( numTrailZeros + 1 );
+		} else {
+		    numSigDigs = 1;
 		}
+		numTrailZeros = 0;
+		state = BAD_OCTAL;
+		break;
+	    } else if (c == '.') {
+		state = FRACTION;
+		break;
+	    } else if (c == 'E' || c == 'e') {
+		state = EXPONENT_START;
+		break;
+	    } 
+#endif
+	    goto endgame;
+
+	    /*
+	     * Scanned 0x.  If state is HEXADECIMAL, scanned at least
+	     * one character following the 0x.  The only acceptable
+	     * inputs are hexadecimal digits.
+	     */
+	case HEXADECIMAL:
+	    acceptState = state;
+	    acceptPoint = p;
+	    acceptLen = len;
+	    /* FALLTHROUGH */
+	case ZERO_X:
+	zerox:
+	    if (c == '0') {
+		++numTrailZeros;
+		state = HEXADECIMAL;
+		break;
+	    } else if (isdigit(UCHAR(c))) {
+		d = (c-'0');
+	    } else if (c >= 'A' && c <= 'F') {
+		d = (c-'A'+10);
+	    } else if (c >= 'a' && c <= 'f') {
+		d = (c-'a'+10);
 	    } else {
-		if (startOfSignificand == NULL) {
-		    startOfSignificand = p;
-		} else if (nTrailZero) {
-		    if (nTrailZero + nSigDigs < DBL_DIG) {
-			exactSignificand *= pow10[nTrailZero];
-		    } else if (nSigDigs < DBL_DIG) {
-			exactSignificand *= pow10[DBL_DIG - nSigDigs];
+		goto endgame;
+	    }
+	    if (objPtr != NULL) {
+		shift = 4 * (numTrailZeros + 1);
+		if (!significandOverflow) {
+		    /*
+		     * Shifting by more bits than are in the value being
+		     * shifted is at least de facto nonportable.  Check
+		     * for too large shifts first.
+		     */
+		    if (significandWide != 0 
+			    && (shift >= CHAR_BIT*sizeof(Tcl_WideUInt)
+			    || significandWide > (~(Tcl_WideUInt)0 >> shift))) {
+			significandOverflow = 1;
+			TclBNInitBignumFromWideUInt(&significandBig,
+						    significandWide);
 		    }
-		    nSigDigs += nTrailZero;
 		}
-		if (nSigDigs < DBL_DIG) {
-		    exactSignificand = 10. * exactSignificand + (c - '0');
+		if (!significandOverflow) {
+		    significandWide
+			= (significandWide << shift) + d;
+		} else {
+		    mp_mul_2d(&significandBig, shift,
+			      &significandBig);
+		    mp_add_d(&significandBig, (mp_digit) d,
+			     &significandBig);
 		}
-		++nSigDigs;
-		nTrailZero = 0;
 	    }
-	    if (seenDp) {
-		++nDigitsAfterDp;
+	    numTrailZeros = 0;
+	    state = HEXADECIMAL;
+	    break;
+
+#ifdef TIP_114_FORMATS
+	case BINARY:
+	    acceptState = state;
+	    acceptPoint = p;
+	    acceptLen = len;
+	case ZERO_B:
+	    if (c == '0') {
+		++numTrailZeros;
+		state = BINARY;
+		break;
+	    } else if (c != '1') {
+		goto endgame;
 	    }
-	    seenDigit = 1;
-	    ++p;
-	} else {
+	    if (objPtr != NULL) {
+		shift = numTrailZeros + 1;
+		if (!significandOverflow) {
+		    /*
+		     * Shifting by more bits than are in the value being
+		     * shifted is at least de facto nonportable.  Check
+		     * for too large shifts first.
+		     */
+		    if (significandWide != 0 
+			    && (shift >= CHAR_BIT*sizeof(Tcl_WideUInt)
+			    || significandWide > (~(Tcl_WideUInt)0 >> shift))) {
+			significandOverflow = 1;
+			TclBNInitBignumFromWideUInt(&significandBig,
+						    significandWide);
+		    }
+		}
+		if (!significandOverflow) {
+		    significandWide
+			= (significandWide << shift) + 1;
+		} else {
+		    mp_mul_2d(&significandBig, shift,
+			      &significandBig);
+		    mp_add_d(&significandBig, (mp_digit) 1,
+			     &significandBig);
+		}
+	    }
+	    numTrailZeros = 0;
+	    state = BINARY;
 	    break;
-	}
-    }
+#endif
 
-    /*
-     * At this point, we've scanned the significand, and p points to the
-     * character beyond it. "startOfSignificand" is the first non-zero
-     * character in the significand. "nSigDigs" is the number of significant
-     * digits of the significand, not including any trailing zeroes.
-     * "exactSignificand" is a floating point number that represents, without
-     * loss of precision, the first min(DBL_DIG,n) digits of the significand.
-     * "nDigitsAfterDp" is the number of digits after the decimal point, again
-     * excluding trailing zeroes.
-     *
-     * Now scan 'E' format
-     */
+	case DECIMAL:
+	    /*
+	     * Scanned an optional + or - followed by a string of
+	     * decimal digits.
+	     */
+#ifdef KILL_OCTAL
+	decimal:
+#endif
+	    acceptState = state;
+	    acceptPoint = p;
+	    acceptLen = len;
+	    if (c == '0') {
+		++numTrailZeros;
+		state = DECIMAL;
+		break;
+	    } else if (isdigit(UCHAR(c))) {
+		if (objPtr != NULL) {
+		    significandOverflow = 
+			AccumulateDecimalDigit((unsigned)(c - '0'), 
+					       numTrailZeros, 
+					       &significandWide,
+					       &significandBig, 
+					       significandOverflow);
+		}
+		numSigDigs += ( numTrailZeros + 1 );
+		numTrailZeros = 0;
+		state = DECIMAL;
+		break;
+	    } else if (flags & TCL_PARSE_INTEGER_ONLY) {
+		goto endgame;
+	    } else if (c == '.') {
+		state = FRACTION;
+		break;
+	    } else if (c == 'E' || c == 'e') {
+		state = EXPONENT_START;
+		break;
+	    }
+	    goto endgame;
 
-    exponent = 0;
-    if (seenDigit && (*p == 'e' || *p == 'E')) {
-	const char* stringSave = p;
-	++p;
-	c = *p;
-	if (isdigit(UCHAR(c)) || c == '+' || c == '-') {
-	    errno = 0;
-	    exponent = strtol(p, (char**)&p, 10);
-	    if (errno == ERANGE) {
-		if (exponent > 0) {
-		    v = HUGE_VAL;
+	    /* 
+	     * Found a decimal point.  If no digits have yet been scanned,
+	     * E is not allowed; otherwise, it introduces the exponent.
+	     * If at least one digit has been found, we have a possible
+	     * complete number.
+	     */
+	case FRACTION:
+	    acceptState = state;
+	    acceptPoint = p;
+	    acceptLen = len;
+	    if (c == 'E' || c=='e') {
+		state = EXPONENT_START;
+		break;
+	    }
+	    /* FALLTHROUGH */
+	case LEADING_RADIX_POINT:
+	    if (c == '0') {
+		++numDigitsAfterDp;
+		++numTrailZeros;
+		state = FRACTION;
+		break;
+	    } else if (isdigit(UCHAR(c))) {
+		++numDigitsAfterDp;
+		if (objPtr != NULL) {
+		    significandOverflow = 
+			AccumulateDecimalDigit((unsigned)(c-'0'),
+					       numTrailZeros, 
+					       &significandWide,
+					       &significandBig, 
+					       significandOverflow);
+		}
+		if ( numSigDigs != 0 ) {
+		    numSigDigs += ( numTrailZeros + 1 );
 		} else {
-		    v = 0.0;
+		    numSigDigs = 1;
 		}
-		*endPtr = p;
-		goto returnValue;
+		numTrailZeros = 0;
+		state = FRACTION;
+		break;
 	    }
-	}
-	if (p == stringSave+1) {
-	    p = stringSave;
-	    exponent = 0;
-	}
-    }
-    exponent += nTrailZero - nDigitsAfterDp;
+	    goto endgame;
 
-    /*
-     * If we come here with no significant digits, we might still be looking
-     * at Inf or NaN.  Go parse them.
-     */
+	case EXPONENT_START:
+	    /* 
+	     * Scanned the E at the start of an exponent. Make sure
+	     * a legal character follows before using the C library
+	     * strtol routine, which allows whitespace.
+	     */
+	    if (c == '+') {
+		state = EXPONENT_SIGNUM;
+		break;
+	    } else if (c == '-') {
+		exponentSignum = 1;
+		state = EXPONENT_SIGNUM;
+		break;
+	    }
+	    /* FALLTHROUGH */
 
-    if (!seenDigit) {
-	/*
-	 * Test for Inf or Infinity (in any case).
-	 */
+	case EXPONENT_SIGNUM:
+	    /*
+	     * Found the E at the start of the exponent, followed by
+	     * a sign character.
+	     */
+	    if (isdigit(UCHAR(c))) {
+		exponent = c - '0';
+		state = EXPONENT;
+		break;
+	    }
+	    goto endgame;
 
-	if (c == 'I' || c == 'i') {
-	    if ((p[1] == 'N' || p[1] == 'n')
-		    && (p[2] == 'F' || p[2] == 'f')) {
-		p += 3;
-		if ((p[0] == 'I' || p[0] == 'i')
-			&& (p[1] == 'N' || p[1] == 'n')
-			&& (p[2] == 'I' || p[2] == 'i')
-			&& (p[3] == 'T' || p[3] == 't')
-			&& (p[4] == 'Y' || p[1] == 'y')) {
-		    p += 5;
-		}
-		errno = ERANGE;
-		v = HUGE_VAL;
-		if (endPtr != NULL) {
-		    *endPtr = p;
+	case EXPONENT:
+	    /* 
+	     * Found an exponent with at least one digit. 
+	     * Accumulate it, making sure to hard-pin it to LONG_MAX
+	     * on overflow.
+	     */
+	    acceptState = state;
+	    acceptPoint = p;
+	    acceptLen = len;
+	    if (isdigit(UCHAR(c))) {
+		if (exponent < (LONG_MAX - 9) / 10) {
+		    exponent = 10 * exponent + (c - '0');
+		} else {
+		    exponent = LONG_MAX;
 		}
-		goto returnValue;
+		state = EXPONENT;
+		break;
 	    }
+	    goto endgame;
 
-#ifdef IEEE_FLOATING_POINT
 	    /*
-	     * Only IEEE floating point supports NaN
+	     * Parse out INFINITY by simply spelling it out.
+	     * INF is accepted as an abbreviation; other prefices are
+	     * not.
 	     */
-	} else if ((c == 'N' || c == 'n')
-		&& (sizeof(Tcl_WideUInt) == sizeof(double))) {
-	    if ((p[1] == 'A' || p[1] == 'a')
-		    && (p[2] == 'N' || p[2] == 'n')) {
-		p += 3;
-
-		if (endPtr != NULL) {
-		    *endPtr = p;
-		}
+	case sI:
+	    if ( c == 'n' || c == 'N' ) {
+		state = sIN;
+		break;
+	    } 
+	    goto endgame;
+	case sIN:
+	    if ( c == 'f' || c == 'F' ) {
+		state = sINF;
+		break;
+	    } 
+	    goto endgame;
+	case sINF:
+	    acceptState = state;
+	    acceptPoint = p;
+	    acceptLen = len;
+	    if ( c == 'i' || c == 'I' ) {
+		state = sINFI;
+		break;
+	    }
+	    goto endgame;
+	case sINFI:
+	    if ( c == 'n' || c == 'N' ) {
+		state = sINFIN;
+		break;
+	    }
+	    goto endgame;
+	case sINFIN:
+	    if ( c == 'i' || c == 'I' ) {
+		state = sINFINI;
+		break;
+	    }
+	    goto endgame;
+	case sINFINI:
+	    if ( c == 't' || c == 'T' ) {
+		state = sINFINIT;
+		break;
+	    }
+	    goto endgame;
+	case sINFINIT:
+	    if ( c == 'y' || c == 'Y' ) {
+		state = sINFINITY;
+		break;
+	    }
+	    goto endgame;
 
-		/*
-		 * Restore FPU mode word.
-		 */
+	    /*
+	     * Parse NaN's.
+	     */
+#ifdef IEEE_FLOATING_POINT
+	case sN:
+	    if ( c == 'a' || c == 'A' ) {
+		state = sNA;
+		break;
+	    }
+	    goto endgame;
+	case sNA:
+	    if ( c == 'n' || c == 'N' ) {
+		state = sNAN;
+		break;
+	    }
+	case sNAN:
+	    acceptState = state;
+	    acceptPoint = p;
+	    acceptLen = len;
+	    if ( c == '(' ) {
+		state = sNANPAREN;
+		break;
+	    }
+	    goto endgame;
 
-		RestoreRoundingMode();
-		return ParseNaN(signum, endPtr);
+	    /*
+	     * Parse NaN(hexdigits)
+	     */
+	case sNANHEX:
+	    if ( c == ')' ) {
+		state = sNANFINISH;
+		break;
 	    }
+	    /* FALLTHROUGH */
+	case sNANPAREN:
+	    if ( isspace(UCHAR(c)) ) {
+		break;
+	    }
+	    if ( numSigDigs < 13 ) {
+		if ( c >= '0' && c <= '9' ) {
+		    d = c - '0';
+		} else if ( c >= 'a' && c <= 'f' ) {
+		    d = 10 + c - 'a';
+		} else if ( c >= 'A' && c <= 'F' ) {
+		    d = 10 + c - 'A';
+		}
+		significandWide = (significandWide << 4) + d;
+		state = sNANHEX;
+		break;
+	    }
+	    goto endgame;
+	case sNANFINISH:
 #endif
+	case sINFINITY:
+	    acceptState = state;
+	    acceptPoint = p;
+	    acceptLen = len;
+	    goto endgame;
 	}
+	++p; 
+	--len;
+    }
+
+  endgame:
 
-	goto error;
+    /* Back up to the last accepting state in the lexer */
+
+    if (acceptState == INITIAL) {
+	status = TCL_ERROR;
     }
+    p = acceptPoint;
+    len = acceptLen;
 
-    /*
-     * We've successfully scanned; update the end-of-element pointer.
-     */
+    /* Skip past trailing whitespace */
 
-    if (endPtr != NULL) {
-	*endPtr = p;
+    if (endPtrPtr != NULL) {
+	*endPtrPtr = p;
     }
 
-    /*
-     * Test for zero.
-     */
+    while (len > 0 && isspace(UCHAR(*p))) {
+	++p;
+	--len;
+    }
 
-    if (nSigDigs == 0) {
-	v = 0.0;
-	goto returnValue;
+    /* Determine whether a partial string is acceptable. */
+
+    if (endPtrPtr == NULL && len != 0 && *p != '\0') {
+	status = TCL_ERROR;
     }
 
-    /*
-     * The easy cases are where we have an exact significand and the exponent
-     * is small enough that we can compute the value with only one roundoff.
-     * In addition to the cases where we can multiply or divide an
-     * exact-integer significand by an exact-integer power of 10, there is
-     * also David Gay's case where we can scale the significand by a power of
-     * 10 (still keeping it exact) and then multiply by an exact power of 10.
-     * The last case enables combinations like 83e25 that would otherwise
-     * require high precision arithmetic.
-     */
+    /* Generate and store the appropriate internal rep */
 
-    if (nSigDigs <= DBL_DIG) {
-	if (exponent >= 0) {
-	    if (exponent <= mmaxpow) {
-		v = exactSignificand * pow10[exponent];
-		goto returnValue;
-	    } else {
-		int diff = DBL_DIG - nSigDigs;
-		if (exponent - diff <= mmaxpow) {
-		    volatile double factor = exactSignificand * pow10[diff];
-		    v = factor * pow10[exponent - diff];
-		    goto returnValue;
+    if (status == TCL_OK && objPtr != NULL) {
+	if ( acceptState != INITIAL ) {
+	    TclFreeIntRep( objPtr );
+	}
+	switch (acceptState) {
+
+	case INITIAL:
+	    status = TCL_ERROR;
+	    break;
+
+	case SIGNUM:
+	case BAD_OCTAL:
+	case ZERO_X:
+#ifdef TIP_114_FORMATS
+	case ZERO_O:
+	case ZERO_B:
+#endif
+	case LEADING_RADIX_POINT:
+	case EXPONENT_START:
+	case EXPONENT_SIGNUM:
+	case sI:
+	case sIN:
+	case sINFI:
+	case sINFIN:
+	case sINFINI:
+	case sINFINIT:
+	case sN:
+	case sNA:
+	case sNANPAREN:
+	case sNANHEX:
+	    panic("in TclParseNumber: bad acceptState, can't happen.");
+
+#ifdef TIP_114_FORMATS
+	case BINARY:
+	    shift = numTrailZeros;
+	    if (!significandOverflow) {
+		if (significandWide !=0 
+			&& (shift >= CHAR_BIT*sizeof(Tcl_WideUInt) 
+			|| significandWide 
+			> (((~(Tcl_WideUInt)0) >> 1) + signum) >> shift )) {
+		    significandOverflow = 1;
+		    TclBNInitBignumFromWideUInt(&significandBig,
+						significandWide);
 		}
 	    }
-	} else if (exponent >= -mmaxpow) {
-	    v = exactSignificand / pow10[-exponent];
-	    goto returnValue;
+	    if (shift) {
+		if ( !significandOverflow ) {
+		    significandWide <<= shift;
+		} else {
+		    mp_mul_2d( &significandBig, shift, &significandBig );
+		}
+	    }
+	    goto returnInteger;
+#endif
+	case HEXADECIMAL:
+	    /* Returning a hex integer.  Final scaling step */
+	    shift = 4 * numTrailZeros;
+	    if (!significandOverflow) {
+		if (significandWide !=0 
+			&& (shift >= CHAR_BIT*sizeof(Tcl_WideUInt) 
+			|| significandWide 
+			> (((~(Tcl_WideUInt)0) >> 1) + signum) >> shift )) {
+		    significandOverflow = 1;
+		    TclBNInitBignumFromWideUInt(&significandBig,
+						significandWide);
+		}
+	    }
+	    if (shift) {
+		if ( !significandOverflow ) {
+		    significandWide <<= shift;
+		} else {
+		    mp_mul_2d( &significandBig, shift, &significandBig );
+		}
+	    }
+	    goto returnInteger;
+
+	case OCTAL:
+	    /* Returning an octal integer.  Final scaling step */
+	    shift = 3 * numTrailZeros;
+	    if (!octalSignificandOverflow) {
+		if (octalSignificandWide != 0
+			&& (shift >= CHAR_BIT*sizeof(Tcl_WideUInt) 
+			|| octalSignificandWide 
+			> (((~(Tcl_WideUInt)0) >> 1) + signum) >> shift )) {
+		    octalSignificandOverflow = 1;
+		    TclBNInitBignumFromWideUInt(&octalSignificandBig,
+						octalSignificandWide);
+		}
+	    }
+	    if ( shift ) {
+		if ( !octalSignificandOverflow ) {
+		    octalSignificandWide <<= shift;
+		} else {
+		    mp_mul_2d( &octalSignificandBig, shift,
+			       &octalSignificandBig );
+		}
+	    }
+	    if (!octalSignificandOverflow) {
+		if (octalSignificandWide > 
+			(Tcl_WideUInt)(((~(unsigned long)0) >> 1) + signum)) {
+#ifndef NO_WIDE_TYPE
+		    if (octalSignificandWide 
+			    <= (((~(Tcl_WideUInt)0) >> 1) + signum)) {
+			objPtr->typePtr = &tclWideIntType;
+			if (signum) {
+			    objPtr->internalRep.wideValue =
+				    - (Tcl_WideInt) octalSignificandWide;
+			} else {
+			    objPtr->internalRep.wideValue =
+				    (Tcl_WideInt) octalSignificandWide;
+			}
+			break;
+		    }
+#endif
+		    TclBNInitBignumFromWideUInt(&octalSignificandBig,
+						octalSignificandWide);
+		    octalSignificandOverflow = 1;
+		} else {
+		    objPtr->typePtr = &tclIntType;
+		    if (signum) {
+			objPtr->internalRep.longValue =
+				- (long) octalSignificandWide;
+		    } else {
+			objPtr->internalRep.longValue =
+				(long) octalSignificandWide;
+		    }
+		}
+	    }
+	    if (octalSignificandOverflow) {
+		if (signum) {
+		    mp_neg(&octalSignificandBig, &octalSignificandBig);
+		}
+		TclSetBignumIntRep(objPtr, &octalSignificandBig);
+	    }		
+	    break;
+
+	case ZERO:
+	case DECIMAL:
+	    significandOverflow =
+		AccumulateDecimalDigit( 0, numTrailZeros-1,
+					&significandWide, &significandBig,
+					significandOverflow );
+	    if (!significandOverflow
+		&& (significandWide
+		    > (((~(Tcl_WideUInt)0) >> 1) + signum))) {
+		significandOverflow = 1;
+		TclBNInitBignumFromWideUInt(&significandBig,
+					    significandWide);
+	    }
+	  returnInteger:
+	    if (!significandOverflow) {
+		if (significandWide > 
+			(Tcl_WideUInt)(((~(unsigned long)0) >> 1) + signum)) {
+#ifndef NO_WIDE_TYPE
+		    if (significandWide 
+			    <= (((~(Tcl_WideUInt)0) >> 1) + signum)) {
+			objPtr->typePtr = &tclWideIntType;
+			if (signum) {
+			    objPtr->internalRep.wideValue =
+				    - (Tcl_WideInt) significandWide;
+			} else {
+			    objPtr->internalRep.wideValue =
+				    (Tcl_WideInt) significandWide;
+			}
+			break;
+		    }
+#endif
+		    TclBNInitBignumFromWideUInt(&significandBig,
+						significandWide);
+		    significandOverflow = 1;
+		} else {
+		    objPtr->typePtr = &tclIntType;
+		    if (signum) {
+			objPtr->internalRep.longValue =
+				- (long) significandWide;
+		    } else {
+			objPtr->internalRep.longValue =
+				(long) significandWide;
+		    }
+		}
+	    }
+	    if (significandOverflow) {
+		if (signum) {
+		    mp_neg(&significandBig, &significandBig);
+		}
+		TclSetBignumIntRep(objPtr, &significandBig);
+	    }
+	    break;
+
+	case FRACTION:
+	case EXPONENT:
+
+	    /*
+	     * Here, we're parsing a floating-point number.
+	     * 'significandWide' or 'significandBig' contains the
+	     * exact significand, according to whether
+	     * 'significandOverflow' is set.  The desired floating
+	     * point value is significand * 10**k, where
+	     * k = numTrailZeros+exponent-numDigitsAfterDp.
+	     */
+
+	    objPtr->typePtr = &tclDoubleType;
+	    if ( exponentSignum ) {
+		exponent = - exponent;
+	    }
+	    if ( !significandOverflow ) {
+		objPtr->internalRep.doubleValue = 
+		    MakeLowPrecisionDouble( signum, 
+					    significandWide,
+					    numSigDigs,
+					    ( numTrailZeros
+					      + exponent
+					      - numDigitsAfterDp ) );
+	    } else {
+		objPtr->internalRep.doubleValue =
+		    MakeHighPrecisionDouble( signum,
+					     &significandBig,
+					     numSigDigs,
+					     ( numTrailZeros
+					       + exponent
+					       - numDigitsAfterDp ) );
+	    }
+	    break;
+
+	case sINF:
+	case sINFINITY:
+	    if ( signum ) {
+		objPtr->internalRep.doubleValue = -HUGE_VAL;
+	    } else {
+		objPtr->internalRep.doubleValue = HUGE_VAL;
+	    }
+	    objPtr->typePtr = &tclDoubleType;
+	    break;
+
+	case sNAN:
+	case sNANFINISH:
+	    objPtr->internalRep.doubleValue
+		= MakeNaN( signum, significandWide );
+	    objPtr->typePtr = &tclDoubleType;
+	    break;
+	    
 	}
     }
 
-    /*
-     * We don't have one of the easy cases, so we can't compute the scanned
-     * number exactly, and have to do it in multiple precision.  Begin by
-     * testing for obvious overflows and underflows.
-     */
+    /* Format an error message when an invalid number is encountered. */
+
+    if ( status != TCL_OK ) {
+	if ( interp != NULL ) {
+	    Tcl_Obj *msg = Tcl_NewStringObj( "expected ", -1 );
+	    Tcl_AppendToObj( msg, type, -1 );
+	    Tcl_AppendToObj( msg, " but got \"", -1 );
+	    TclAppendLimitedToObj( msg, string, length, 50, "" );
+	    Tcl_AppendToObj( msg, "\"", -1 );
+	    if ( state == BAD_OCTAL ) {
+		Tcl_AppendToObj( msg, " (looks like invalid octal number)",
+				 -1 );
+	    }
+	    Tcl_SetObjResult( interp, msg );
+	}
+    }
 
-    if (nSigDigs + exponent - 1 > maxDigits) {
-	v = HUGE_VAL;
-	errno = ERANGE;
-	goto returnValue;
+    /* Free memory */
+
+    if (octalSignificandOverflow) {
+	mp_clear(&octalSignificandBig);
     }
-    if (nSigDigs + exponent - 1 < minDigits) {
-	errno = ERANGE;
-	v = 0.;
-	goto returnValue;
+    if (significandOverflow) {
+	mp_clear(&significandBig);
     }
+    return status;
+}
+
+/*
+ *----------------------------------------------------------------------
+ *
+ * AccumulateDecimalDigit --
+ *
+ *	Consume a decimal digit in a number being scanned.
+ *
+ * Results:
+ *	Returns 1 if the number has overflowed to a bignum, 0 if it
+ *	still fits in a wide integer.
+ *
+ * Side effects:
+ *	Updates either the wide or bignum representation.
+ *
+ *----------------------------------------------------------------------
+ */
 
-    /*
-     * Nothing exceeds the boundaries of the tables, at least.  Compute an
-     * approximate value for the number, with no possibility of overflow
-     * because we manage the exponent separately.
-     */
-
-    if (nSigDigs > DBL_DIG) {
-	expt2 = exponent + nSigDigs - DBL_DIG;
-    } else {
-	expt2 = exponent;
-    }
-    v = frexp(exactSignificand, &machexp);
-    if (expt2 > 0) {
-	v = frexp(v * pow10[expt2 & 0xf], &j);
-	machexp += j;
-	for (i=4 ; i<9 ; ++i) {
-	    if (expt2 & (1 << i)) {
-		v = frexp(v * pow_10_2_n[i], &j);
-		machexp += j;
+static int
+AccumulateDecimalDigit( unsigned digit,
+				/* Digit being scanned */
+			int numZeros,
+				/* Count of zero digits preceding the
+				 * digit being scanned */
+			Tcl_WideUInt* wideRepPtr,
+				/* Representation of the partial number
+				 * as a wide integer */
+			mp_int* bignumRepPtr,
+				/* Representation of the partial number
+				 * as a bignum */
+			int bignumFlag )
+				/* Flag == 1 if the number overflowed
+				 * previous to this digit. */
+{
+    int i, n;
+
+    /* Check if the number still fits in a wide */
+
+    if (!bignumFlag) {
+	if (*wideRepPtr != 0) {
+	    if ((numZeros >= maxpow10_wide)
+		|| (*wideRepPtr > (((~(Tcl_WideUInt)0) - digit)
+				   / pow10_wide[numZeros+1]))) {
+		/* Oops, it's overflowed, have to allocate a bignum */
+		TclBNInitBignumFromWideUInt (bignumRepPtr, *wideRepPtr);
+		bignumFlag = 1;
 	    }
 	}
+    }
+
+    /* Multiply the number by 10**numZeros+1 and add in the new digit. */
+
+    if (!bignumFlag) {
+
+	/* Wide multiplication */
+
+	*wideRepPtr = *wideRepPtr * pow10_wide[numZeros+1] + digit;
+    } else if (numZeros < log10_DIGIT_MAX ) {
+
+	/* Up to about 8 zeros - single digit multiplication */
+
+	mp_mul_d (bignumRepPtr, (mp_digit) pow10_wide[numZeros+1],
+		  bignumRepPtr);
+	mp_add_d (bignumRepPtr, (mp_digit) digit, bignumRepPtr);
+
     } else {
-	v = frexp(v / pow10[(-expt2) & 0xf], &j);
-	machexp += j;
-	for (i=4 ; i<9 ; ++i) {
-	    if ((-expt2) & (1 << i)) {
-		v = frexp(v / pow_10_2_n[i], &j);
-		machexp += j;
+
+	/* 
+	 * More than single digit multiplication. Multiply by the appropriate
+	 * small powers of 5, and then shift.  Large strings of zeroes are
+	 * eaten 256 at a time; this is less efficient than it could be,
+	 * but seems implausible.  We presume that DIGIT_BIT is at least 27.
+	 * The first multiplication, by up to 10**7, is done with a
+	 * one-DIGIT multiply (this presumes that DIGIT_BIT >= 24).
+	 */
+
+	n = numZeros + 1;
+	mp_mul_d (bignumRepPtr, (mp_digit) pow10_wide[n&0x7], bignumRepPtr);
+	for (i = 3; i <= 7; ++i) {
+	    if (n & (1 << i)) {
+		mp_mul (bignumRepPtr, pow5+i, bignumRepPtr);
 	    }
 	}
+	while (n >= 256) {
+	    mp_mul (bignumRepPtr, pow5+8, bignumRepPtr);
+	    n -= 256;
+	}
+	mp_mul_2d (bignumRepPtr, (int)(numZeros+1)&~0x7, bignumRepPtr);
     }
 
+    return bignumFlag;
+}
+
+/*
+ *----------------------------------------------------------------------
+ *
+ * MakeLowPrecisionDouble --
+ *
+ *	Makes the double precision number, signum*significand*10**exponent.
+ *
+ * Results:
+ *	Returns the constructed number.
+ *
+ * Common cases, where there are few enough digits that the number can
+ * be represented with at most roundoff, are handled specially here.
+ * If the number requires more than one rounded operation to compute,
+ * the code promotes the significand to a bignum and calls
+ * MakeHighPrecisionDouble to do it instead.
+ *
+ *----------------------------------------------------------------------
+ */
+
+static double
+MakeLowPrecisionDouble( int signum,
+				/* 1 if the number is negative, 0 otherwise */
+			Tcl_WideUInt significand,
+				/* Significand of the number */
+			int numSigDigs,
+				/* Number of digits in the significand */
+			int exponent )
+				/* Power of ten */
+{
+    double retval;		/* Value of the number */
+    mp_int significandBig;	/* Significand expressed as a bignum */
+
     /*
-     * A first approximation is that the result will be v * 2 ** machexp. v is
-     * greater than or equal to 0.5 and less than 1. If machexp >
-     * DBL_MAX_EXP*log2(FLT_RADIX), there is an overflow. Constrain the result
-     * to the smallest representible number to avoid premature underflow.
+     * With gcc on x86, the floating point rounding mode is double-extended.
+     * This causes the result of double-precision calculations to be rounded
+     * twice: once to the precision of double-extended and then again to the
+     * precision of double.  Double-rounding introduces gratuitous errors of
+     * 1 ulp, so we need to change rounding mode to 53-bits.
      */
 
-    if (machexp > DBL_MAX_EXP * log2FLT_RADIX) {
-	v = HUGE_VAL;
-	errno = ERANGE;
-	goto returnValue;
-    }
+#if defined(__GNUC__) && defined(__i386)
+    fpu_control_t roundTo53Bits = 0x027f;
+    fpu_control_t oldRoundingMode;
+    _FPU_GETCW( oldRoundingMode );
+    _FPU_SETCW( roundTo53Bits );
+#endif
 
-    v = SafeLdExp(v, machexp);
-    if (v < tiny) {
-	v = tiny;
-    }
+    /* Test for the easy cases */
 
-    /*
-     * We have a first approximation in v. Now we need to refine it.
-     */
+    if ( numSigDigs <= DBL_DIG ) {
+	if ( exponent >= 0 ) {
+	    if ( exponent <= mmaxpow ) {
 
-    v = RefineResult(v, startOfSignificand, nSigDigs, exponent);
+		/* 
+		 * The significand is an exact integer, and so is 
+		 * 10**exponent. The product will be correct to within
+		 * 1/2 ulp without special handling.
+		 */
 
-    /*
-     * In a very few cases, a second iteration is needed. e.g., 457e-102
-     */
+		retval = (double)(Tcl_WideInt)significand * pow10[ exponent ];
+		goto returnValue;
 
-    v = RefineResult(v, startOfSignificand, nSigDigs, exponent);
+	    } else {
+		int diff = DBL_DIG - numSigDigs;
+		if ( exponent-diff <= mmaxpow ) {
+
+		    /* 
+		     * 10**exponent is not an exact integer, but
+		     * 10**(exponent-diff) is exact, and so is
+		     * significand*10**diff, so we can still compute
+		     * the value with only one roundoff.
+		     */
+		    volatile double factor
+			= (double)(Tcl_WideInt)significand * pow10[diff];
+		    retval = factor * pow10[exponent-diff];
+		    goto returnValue;
+		}
+	    }
+	} else {
+	    if ( exponent >= -mmaxpow ) {
+
+		/* 
+		 * 10**-exponent is an exact integer, and so is the
+		 * significand. Compute the result by one division,
+		 * again with only one rounding.
+		 */
+
+		retval = (double)(Tcl_WideInt)significand / pow10[-exponent];
+		goto returnValue;
+	    }
+	}
+    }
 
     /*
-     * Handle underflow.
+     * All the easy cases have failed.  Promote ths significand
+     * to bignum and call MakeHighPrecisionDouble to do it the hard way.
      */
 
+    TclBNInitBignumFromWideUInt (&significandBig, significand);
+    retval = MakeHighPrecisionDouble( 0, &significandBig, numSigDigs,
+				      exponent );
+    
+    /* Come here to return the computed value */
+
   returnValue:
-    if (nSigDigs != 0 && v == 0.0) {
-	errno = ERANGE;
+
+    if ( signum ) {
+	retval = -retval;
     }
 
+    /* On gcc on x86, restore the floating point mode word. */
+
+#if defined(__GNUC__) && defined(__i386)
+    _FPU_SETCW( oldRoundingMode );
+#endif
+
+    return retval;
+}
+
+/*
+ *----------------------------------------------------------------------
+ *
+ * MakeHighPrecisionDouble --
+ *
+ *	Makes the double precision number, signum*significand*10**exponent.
+ *
+ * Results:
+ *	Returns the constructed number.
+ *
+ * MakeHighPrecisionDouble is used when arbitrary-precision arithmetic
+ * is needed to ensure correct rounding.  It begins by calculating a
+ * low-precision approximation to the desired number, and then refines
+ * the answer in high precision.
+ *
+ *----------------------------------------------------------------------
+ */
+
+static double
+MakeHighPrecisionDouble( int signum,
+				/* 1=negative, 0=nonnegative */
+			 mp_int* significand,
+				/* Exact significand of the number */
+			 int numSigDigs,
+				/* Number of significant digits */
+			 int exponent )
+				/* Power of 10 by which to multiply */
+{
+
+    double retval;
+    int machexp;		/* Machine exponent of a power of 10 */
+
     /*
-     * Return a number with correct sign.
+     * With gcc on x86, the floating point rounding mode is double-extended.
+     * This causes the result of double-precision calculations to be rounded
+     * twice: once to the precision of double-extended and then again to the
+     * precision of double.  Double-rounding introduces gratuitous errors of
+     * 1 ulp, so we need to change rounding mode to 53-bits.
      */
 
-    if (signum) {
-	v = -v;
+#if defined(__GNUC__) && defined(__i386)
+    fpu_control_t roundTo53Bits = 0x027f;
+    fpu_control_t oldRoundingMode;
+    _FPU_GETCW( oldRoundingMode );
+    _FPU_SETCW( roundTo53Bits );
+#endif
+
+    /* Quick checks for over/underflow */
+
+    if ( numSigDigs + exponent - 1 > maxDigits ) {
+	retval = HUGE_VAL;
+	goto returnValue;
+    }
+    if ( numSigDigs + exponent - 1 < minDigits ) {
+	retval = 0;
+	goto returnValue;
     }
 
-    /*
-     * Restore FPU mode word and return.
+    /* 
+     * Develop a first approximation to the significand. It is tempting
+     * simply to force bignum to double, but that will overflow on input
+     * numbers like 1.[string repeat 0 1000]1; while this is a not terribly
+     * likely scenario, we still have to deal with it. Use fraction and
+     * exponent instead.  Once we have the significand, multiply by
+     * 10**exponent. Test for overflow. Convert back to a double, and
+     * test for underflow.
      */
 
-    RestoreRoundingMode();
-    return v;
+    retval = BignumToBiasedFrExp( significand, &machexp );
+    retval = Pow10TimesFrExp( exponent, retval, &machexp );
+    if ( machexp > DBL_MAX_EXP * log2FLT_RADIX ) {
+	retval = HUGE_VAL;
+	goto returnValue;
+    }
+    retval = SafeLdExp( retval, machexp );
+    if ( retval < tiny ) {
+	retval = tiny;
+    }
 
-    /*
-     * Come here on an invalid input.
+    /* 
+     * Refine the result twice. (The second refinement should be 
+     * necessary only if the best approximation is a power of 2
+     * minus 1/2 ulp).
      */
 
-  error:
-    if (endPtr != NULL) {
-	*endPtr = s;
+    retval = RefineApproximation( retval, significand, exponent );
+    retval = RefineApproximation( retval, significand, exponent );
+
+    /* Come here to return the computed value */
+
+  returnValue:
+    if ( signum ) {
+	retval = -retval;
     }
 
-    /*
-     * Restore FPU mode word and return.
-     */
+    /* On gcc on x86, restore the floating point mode word. */
 
-    RestoreRoundingMode();
-    return 0.0;
+#if defined(__GNUC__) && defined(__i386)
+    _FPU_SETCW( oldRoundingMode );
+#endif
+    return retval;
 }
 
 /*
  *----------------------------------------------------------------------
  *
- * RefineResult --
+ * MakeNaN --
+ *
+ *	Makes a "Not a Number" given a set of bits to put in the
+ *	tag bits
  *
- *	Given a poor approximation to a floating point number, returns a
- *	better one. (The better approximation is correct to within 1 ulp, and
- *	is entirely correct if the poor approximation is correct to 1 ulp.)
+ * Note that a signalling NaN is never returned.
+ *
+ *----------------------------------------------------------------------
+ */
+
+#ifdef IEEE_FLOATING_POINT
+static double
+MakeNaN( int signum,		/* Sign bit (1=negative, 0=nonnegative */
+	 Tcl_WideUInt tags ) 	/* Tag bits to put in the NaN */
+{
+    union {
+	Tcl_WideUInt iv;
+	double dv;
+    } theNaN;
+
+    theNaN.iv = tags;
+    theNaN.iv &= ( ((Tcl_WideUInt) 1) << 51 ) - 1;
+    if ( signum ) {
+	theNaN.iv |= ((Tcl_WideUInt) (0x8000 | NAN_START)) << 48;
+    } else {
+	theNaN.iv |= ((Tcl_WideUInt) NAN_START) << 48;
+    }
+
+    return theNaN.dv;
+}
+#endif
+
+/*
+ *----------------------------------------------------------------------
+ *
+ * RefineApproximation --
+ *
+ *	Given a poor approximation to a floating point number, returns
+ *	a better one (The better approximation is correct to within
+ *	1 ulp, and is entirely correct if the poor approximation is
+ *	correct to 1 ulp.)
  *
  * Results:
  *	Returns the improved result.
@@ -534,119 +1528,109 @@ TclStrToD(CONST char *s,	/* String to scan. */
  */
 
 static double
-RefineResult(double approxResult, /* Approximate result of conversion. */
-	     CONST char* sigStart,
-				/* Pointer to start of significand in input
-				 * string. */
-	     int nSigDigs,	/* Number of significant digits. */
-	     long exponent)	/* Power of ten to multiply by significand. */
+RefineApproximation( double approxResult, 
+				/* Approximate result of conversion */
+		     mp_int* exactSignificand,
+				/* Integer significand */
+		     int exponent )
+				/* Power of 10 to multiply by significand */
 {
-    int M2, M5;			/* Powers of 2 and of 5 needed to put the
-				 * decimal and binary numbers over a common
-				 * denominator. */
-    double significand;		/* Sigificand of the binary number. */
-    int binExponent;		/* Exponent of the binary number. */
+
+    int M2, M5;			/*  Powers of 2 and of 5 needed to put
+				 * the decimal and binary numbers over
+				 * a common denominator. */
+    double significand;		/* Sigificand of the binary number */
+    int binExponent;		/* Exponent of the binary number */
+
     int msb;			/* Most significant bit position of an
-				 * intermediate result. */
+				 * intermediate result */
     int nDigits;		/* Number of mp_digit's in an intermediate
-				 * result. */
-    mp_int twoMv;		/* Approx binary value expressed as an exact
-				 * integer scaled by the multiplier 2M. */
-    mp_int twoMd;		/* Exact decimal value expressed as an exact
-				 * integer scaled by the multiplier 2M. */
-    int scale;			/* Scale factor for M. */
-    int multiplier;		/* Power of two to scale M. */
-    double num, den;		/* Numerator and denominator of the correction
-				 * term. */
-    double quot;		/* Correction term. */
-    double minincr;		/* Lower bound on the absolute value of the
-				 * correction term. */
+				 * result */
+    mp_int twoMv;		/* Approx binary value expressed as an
+				 * exact integer scaled by the multiplier 2M */
+    mp_int twoMd;		/* Exact decimal value expressed as an
+				 * exact integer scaled by the multiplier 2M */
+    int scale;			/* Scale factor for M */
+    int multiplier;		/* Power of two to scale M */
+    double num, den;		/* Numerator and denominator of the
+				 * correction term */
+    double quot;		/* Correction term */
+    double minincr;		/* Lower bound on the absolute value
+				 * of the correction term. */
     int i;
-    const char* p;
 
     /*
-     * The first approximation is always low. If we find that it's HUGE_VAL,
-     * we're done.
+     * The first approximation is always low.  If we find that
+     * it's HUGE_VAL, we're done.
      */
 
-    if (approxResult == HUGE_VAL) {
+    if ( approxResult == HUGE_VAL ) {
 	return approxResult;
     }
 
     /*
-     * Find a common denominator for the decimal and binary fractions. The
-     * common denominator will be 2**M2 + 5**M5.
+     * Find a common denominator for the decimal and binary fractions.
+     * The common denominator will be 2**M2 + 5**M5.
      */
 
-    significand = frexp(approxResult, &binExponent);
+    significand = frexp( approxResult, &binExponent );
     i = mantBits - binExponent;
-    if (i < 0) {
+    if ( i < 0 ) {
 	M2 = 0;
     } else {
 	M2 = i;
     }
-    if (exponent > 0) {
+    if ( exponent > 0 ) {
 	M5 = 0;
     } else {
 	M5 = -exponent;
-	if ((M5-1) > M2) {
+	if ( (M5-1) > M2 ) {
 	    M2 = M5-1;
 	}
     }
 
-    /*
-     * The floating point number is significand*2**binExponent. The 2**-1 bit
-     * of the significand (the most significant) corresponds to the
-     * 2**(binExponent+M2 + 1) bit of 2*M2*v. Allocate enough digits to hold
-     * that quantity, then convert the significand to a large integer, scaled
+    /* 
+     * The floating point number is significand*2**binExponent.
+     * Compute the large integer significand*2**(binExponent+M2+1)
+     * The 2**-1 bit of the significand (the most significant) 
+     * corresponds to the 2**(binExponent+M2 + 1) bit of 2*M2*v.
+     * Allocate enough digits to hold that quantity, then
+     * convert the significand to a large integer, scaled
      * appropriately. Then multiply by the appropriate power of 5.
      */
 
-    msb = binExponent + M2;	/* 1008 */
+    msb = binExponent + M2;  /* 1008 */
     nDigits = msb / DIGIT_BIT + 1;
-    mp_init_size(&twoMv, nDigits);
-    i = (msb % DIGIT_BIT + 1);
+    mp_init_size( &twoMv, nDigits );
+    i = ( msb % DIGIT_BIT + 1 ); 
     twoMv.used = nDigits;
-    significand *= SafeLdExp(1.0, i);
-    while (--nDigits >= 0) {
+    significand *= SafeLdExp( 1.0, i );
+    while ( -- nDigits >= 0 ) {
 	twoMv.dp[nDigits] = (mp_digit) significand;
 	significand -= (mp_digit) significand;
-	significand = SafeLdExp(significand, DIGIT_BIT);
+	significand = SafeLdExp( significand, DIGIT_BIT );
     }
-    for (i=0 ; i<=8 ; ++i) {
-	if (M5 & (1 << i)) {
-	    mp_mul(&twoMv, pow5+i, &twoMv);
+    for ( i = 0; i <= 8; ++i ) {
+	if ( M5 & ( 1 << i ) ) {
+	    mp_mul( &twoMv, pow5+i, &twoMv );
 	}
     }
-
-    /*
-     * Collect the decimal significand as a high precision integer. The least
-     * significant bit corresponds to bit M2+exponent+1 so it will need to be
-     * shifted left by that many bits after being multiplied by
-     * 5**(M5+exponent).
+    
+    /* 
+     * Collect the decimal significand as a high precision integer.
+     * The least significant bit corresponds to bit M2+exponent+1
+     * so it will need to be shifted left by that many bits after
+     * being multiplied by 5**(M5+exponent).
      */
 
-    mp_init(&twoMd);
-    mp_zero(&twoMd);
-    i = nSigDigs;
-    for (p=sigStart ;; ++p) {
-	char c = *p;
-	if (isdigit(UCHAR(c))) {
-	    mp_mul_d(&twoMd, (unsigned) 10, &twoMd);
-	    mp_add_d(&twoMd, (unsigned) (c - '0'), &twoMd);
-	    --i;
-	    if (i == 0) {
-		break;
-	    }
-	}
-    }
-    for (i=0 ; i<=8 ; ++i) {
-	if ((M5+exponent) & (1 << i)) {
-	    mp_mul(&twoMd, pow5+i, &twoMd);
+    mp_init_copy( &twoMd, exactSignificand );
+    for ( i = 0; i <= 8; ++i ) {
+	if ( (M5+exponent) & ( 1 << i ) ) {
+	    mp_mul( &twoMd, pow5+i, &twoMd );
 	}
     }
-    mp_mul_2d(&twoMd, M2+exponent+1, &twoMd);
-    mp_sub(&twoMd, &twoMv, &twoMd);
+    mp_mul_2d( &twoMd, M2+exponent+1, &twoMd );
+    mp_sub( &twoMd, &twoMv, &twoMd );
 
     /*
      * The result, 2Mv-2Md, needs to be divided by 2M to yield a correction
@@ -656,135 +1640,55 @@ RefineResult(double approxResult, /* Approximate result of conversion. */
 
     scale = binExponent - mantBits - 1;
 
-    mp_set(&twoMv, 1);
-    for (i=0 ; i<=8 ; ++i) {
-	if (M5 & (1 << i)) {
-	    mp_mul(&twoMv, pow5+i, &twoMv);
+    mp_set( &twoMv, 1 );
+    for ( i = 0; i <= 8; ++i ) {
+	if ( M5 & ( 1 << i ) ) {
+	    mp_mul( &twoMv, pow5+i, &twoMv );
 	}
     }
     multiplier = M2 + scale + 1;
-    if (multiplier > 0) {
-	mp_mul_2d(&twoMv, multiplier, &twoMv);
-    } else if (multiplier < 0) {
-	mp_div_2d(&twoMv, -multiplier, &twoMv, NULL);
+    if ( multiplier > 0 ) {
+	mp_mul_2d( &twoMv, multiplier, &twoMv );
+    } else if ( multiplier < 0 ) {
+	mp_div_2d( &twoMv, -multiplier, &twoMv, NULL );
     }
 
     /*
-     * If the result is less than unity, the error is less than 1/2 unit in
-     * the last place, so there's no correction to make.
+     * If the result is less than unity, the error is less than 1/2 unit
+     * in the last place, so there's no correction to make.
      */
 
-    if (mp_cmp_mag(&twoMd, &twoMv) == MP_LT) {
-	mp_clear(&twoMd);
-	mp_clear(&twoMv);
+    if ( mp_cmp_mag( &twoMd, &twoMv ) == MP_LT ) {
 	return approxResult;
     }
 
-    /*
-     * Convert the numerator and denominator of the corrector term accurately
-     * to floating point numbers.
+    /* 
+     * Convert the numerator and denominator of the corrector term
+     * accurately to floating point numbers.
      */
 
-    num = TclBignumToDouble(&twoMd);
-    den = TclBignumToDouble(&twoMv);
+    num = TclBignumToDouble( &twoMd );
+    den = TclBignumToDouble( &twoMv );
 
-    quot = SafeLdExp(num/den, scale);
-    minincr = SafeLdExp(1.0, binExponent - mantBits);
+    quot = SafeLdExp( num/den, scale );
+    minincr = SafeLdExp( 1.0, binExponent - mantBits );
 
-    if (quot<0. && quot>-minincr) {
+    if ( quot < 0. && quot > -minincr ) {
 	quot = -minincr;
-    } else if (quot>0. && quot<minincr) {
+    } else if ( quot > 0. && quot < minincr ) {
 	quot = minincr;
     }
 
-    mp_clear(&twoMd);
-    mp_clear(&twoMv);
+    mp_clear( &twoMd );
+    mp_clear( &twoMv );
 
+    
     return approxResult + quot;
 }
 
 /*
  *----------------------------------------------------------------------
  *
- * ParseNaN --
- *
- *	Parses a "not a number" from an input string, and returns the double
- *	precision NaN corresponding to it.
- *
- * Side effects:
- *	Advances endPtr to follow any (hex) in the input string.
- *
- *	If the NaN is followed by a left paren, a string of spaes and
- *	hexadecimal digits, and a right paren, endPtr is advanced to follow
- *	it.
- *
- *	The string of hexadecimal digits is OR'ed into the resulting NaN, and
- *	the signum is set as well.  Note that a signalling NaN is never
- *	returned.
- *
- *----------------------------------------------------------------------
- */
-
-static double
-ParseNaN(int signum,		/* Flag == 1 if minus sign has been seen in
-				 * front of NaN. */
-	 CONST char** endPtr)	/* Pointer-to-pointer to char following "NaN"
-				 * in the input string. */
-{
-    const char* p = *endPtr;
-    char c;
-    union {
-	Tcl_WideUInt iv;
-	double dv;
-    } theNaN;
-
-    /*
-     * Scan off a hex number in parentheses.  Embedded blanks are ok.
-     */
-
-    theNaN.iv = 0;
-    if (*p == '(') {
-	++p;
-	for (;;) {
-	    c = *p++;
-	    if (isspace(UCHAR(c))) {
-		continue;
-	    } else if (c == ')') {
-		*endPtr = p;
-		break;
-	    } else if (isdigit(UCHAR(c))) {
-		c -= '0';
-	    } else if (c >= 'A' && c <= 'F') {
-		c -= 'A' + 10;
-	    } else if (c >= 'a' && c <= 'f') {
-		c -= 'a' + 10;
-	    } else {
-		theNaN.iv = (((Tcl_WideUInt) NAN_START) << 48)
-			| (((Tcl_WideUInt) signum) << 63);
-		return theNaN.dv;
-	    }
-	    theNaN.iv = (theNaN.iv << 4) | c;
-	}
-    }
-
-    /*
-     * Mask the hex number down to the least significant 51 bits.
-     */
-
-    theNaN.iv &= (((Tcl_WideUInt) 1) << 51) - 1;
-    if (signum) {
-	theNaN.iv |= ((Tcl_WideUInt) 0xfff8) << 48;
-    } else {
-	theNaN.iv |= ((Tcl_WideUInt) NAN_START) << 48;
-    }
-
-    *endPtr = p;
-    return theNaN.dv;
-}
-
-/*
- *----------------------------------------------------------------------
- *
  * TclDoubleDigits --
  *
  *	Converts a double to a string of digits.
@@ -803,20 +1707,21 @@ ParseNaN(int signum,		/* Flag == 1 if minus sign has been seen in
  */
 
 int
-TclDoubleDigits(char * strPtr,	/* Buffer in which to store the result, must
-				 * have at least 18 chars. */
-		double v,	/* Number to convert. Must be finite, and not
-				 * NaN. */
-		int *signum)	/* Output: 1 if the number is negative.
-				 * Should handle -0 correctly on the IEEE
-				 * architecture. */
+TclDoubleDigits( char * string,	/* Buffer in which to store the result,
+				 * must have at least 18 chars */
+		 double v,	/* Number to convert. Must be
+				 * finite, and not NaN */
+		 int *signum )	/* Output: 1 if the number is negative.
+				 * Should handle -0 correctly on the
+				 * IEEE architecture. */
 {
-    double f;			/* Significand of v. */
     int e;			/* Power of FLT_RADIX that satisfies
 				 * v = f * FLT_RADIX**e */
     int lowOK, highOK;
     mp_int r;			/* Scaled significand. */
     mp_int s;			/* Divisor such that v = r / s */
+    int smallestSig;		/* Flag == 1 iff v's significand is
+				 * the smallest that can be represented. */
     mp_int mplus;		/* Scaled epsilon: (r + 2* mplus) == v(+)
 				 * where v(+) is the floating point successor
 				 * of v. */
@@ -830,103 +1735,31 @@ TclDoubleDigits(char * strPtr,	/* Buffer in which to store the result, must
     int sfac5 = 0;
     int mplusfac2 = 0;
     int mminusfac2 = 0;
-    double a;
     char c;
     int i, k, n;
 
-    /*
-     * Take the absolute value of the number, and report the number's sign.
-     * Take special steps to preserve signed zeroes in IEEE floating point.
-     * (We can't use fpclassify, because that's a C9x feature and we still
-     * have to build on C89 compilers.)
-     */
+    /* Split the number into absolute value and signum. */
 
-#ifndef IEEE_FLOATING_POINT
-    if (v >= 0.0) {
-	*signum = 0;
-    } else {
-	*signum = 1;
-	v = -v;
-    }
-#else
-    union {
-	Tcl_WideUInt iv;
-	double dv;
-    } bitwhack;
-    bitwhack.dv = v;
-    if (bitwhack.iv & ((Tcl_WideUInt) 1 << 63)) {
-	*signum = 1;
-	bitwhack.iv &= ~((Tcl_WideUInt) 1 << 63);
-	v = bitwhack.dv;
-    } else {
-	*signum = 0;
-    }
-#endif
+    v = AbsoluteValue(v, signum);
 
     /*
      * Handle zero specially.
      */
 
-    if (v == 0.0) {
-	*strPtr++ = '0';
-	*strPtr++ = '\0';
+    if ( v == 0.0 ) {
+	*string++ = '0';
+	*string++ = '\0';
 	return 1;
     }
 
-    /*
-     * Develop f and e such that v = f * FLT_RADIX**e, with
-     * 1.0/FLT_RADIX <= f < 1.
+    /* 
+     * Find a large integer r, and integer e, such that 
+     *         v = r * FLT_RADIX**e
+     * and r is as small as possible.  Also determine whether the
+     * significand is the smallest possible.
      */
 
-    f = frexp(v, &e);
-    n = e % log2FLT_RADIX;
-    if (n > 0) {
-	n -= log2FLT_RADIX;
-	e += 1;
-    }
-    f *= ldexp(1.0, n);
-    e = (e - n) / log2FLT_RADIX;
-    if (f == 1.0) {
-	f = 1.0 / FLT_RADIX;
-	e += 1;
-    }
-
-    /*
-     * If the original number was denormalized, adjust e and f to be denormal
-     * as well.
-     */
-
-    if (e < DBL_MIN_EXP) {
-	n = mantBits + (e - DBL_MIN_EXP)*log2FLT_RADIX;
-	f = ldexp(f, (e - DBL_MIN_EXP)*log2FLT_RADIX);
-	e = DBL_MIN_EXP;
-	n = (n + DIGIT_BIT - 1) / DIGIT_BIT;
-    } else {
-	n = mantDIGIT;
-    }
-
-    /*
-     * Now extract the base-2**DIGIT_BIT digits of f into a multi-precision
-     * integer r. Preserve the invariant v = r * 2**rfac2 * FLT_RADIX**e by
-     * adjusting e.
-     */
-
-    a = f;
-    n = mantDIGIT;
-    mp_init_size(&r, n);
-    r.used = n;
-    r.sign = MP_ZPOS;
-    i = (mantBits % DIGIT_BIT);
-    if (i == 0) {
-	i = DIGIT_BIT;
-    }
-    while (n > 0) {
-	a *= ldexp(1.0, i);
-	i = DIGIT_BIT;
-	r.dp[--n] = (mp_digit) a;
-	a -= (mp_digit) a;
-    }
-    e -= DBL_MANT_DIG;
+    smallestSig = GetIntegerTimesPower(v, &r, &e);
 
     lowOK = highOK = (mp_iseven(&r));
 
@@ -943,12 +1776,12 @@ TclDoubleDigits(char * strPtr,	/* Buffer in which to store the result, must
     if (e >= 0) {
 	int bits = e * log2FLT_RADIX;
 
-	if (f != 1.0/FLT_RADIX) {
+	if (!smallestSig) {
 	    /*
 	     * Normal case, m+ and m- are both FLT_RADIX**e
 	     */
 
-	    rfac2 += bits + 1;
+	    rfac2 = bits + 1;
 	    sfac2 = 1;
 	    mplusfac2 = bits;
 	    mminusfac2 = bits;
@@ -959,7 +1792,7 @@ TclDoubleDigits(char * strPtr,	/* Buffer in which to store the result, must
 	     * smaller exponent when going to e's predecessor.
 	     */
 
-	    rfac2 += bits + log2FLT_RADIX + 1;
+	    rfac2 = bits + log2FLT_RADIX + 1;
 	    sfac2 = 1 + log2FLT_RADIX;
 	    mplusfac2 = bits + log2FLT_RADIX;
 	    mminusfac2 = bits;
@@ -969,13 +1802,13 @@ TclDoubleDigits(char * strPtr,	/* Buffer in which to store the result, must
 	 * v has digits after the binary point
 	 */
 
-	if (e <= DBL_MIN_EXP-DBL_MANT_DIG || f != 1.0/FLT_RADIX) {
+	if (e <= DBL_MIN_EXP-DBL_MANT_DIG || !smallestSig) {
 	    /*
 	     * Either f isn't the smallest significand or e is the smallest
 	     * exponent.  mplus and mminus will both be 1.
 	     */
 
-	    rfac2 += 1;
+	    rfac2 = 1;
 	    sfac2 = 1 - e * log2FLT_RADIX;
 	    mplusfac2 = 0;
 	    mminusfac2 = 0;
@@ -986,7 +1819,7 @@ TclDoubleDigits(char * strPtr,	/* Buffer in which to store the result, must
 	     * fact that v's predecessor has a smaller exponent.
 	     */
 
-	    rfac2 += 1 + log2FLT_RADIX;
+	    rfac2 = 1 + log2FLT_RADIX;
 	    sfac2 = 1 + log2FLT_RADIX * (1 - e);
 	    mplusfac2 = FLT_RADIX;
 	    mminusfac2 = 0;
@@ -1081,9 +1914,9 @@ TclDoubleDigits(char * strPtr,	/* Buffer in which to store the result, must
 	} else {
 	    tc2= (tc2 > 0);
 	}
-	if (!tc1) {
-	    if (!tc2) {
-		*strPtr++ = '0' + i;
+	if ( ! tc1 ) {
+	    if ( !tc2 ) {
+		*string++ = '0' + i;
 	    } else {
 		c = (char) (i + '1');
 		break;
@@ -1103,8 +1936,8 @@ TclDoubleDigits(char * strPtr,	/* Buffer in which to store the result, must
 	    break;
 	}
     };
-    *strPtr++ = c;
-    *strPtr++ = '\0';
+    *string++ = c;
+    *string++ = '\0';
 
     /*
      * Free memory, and return.
@@ -1117,6 +1950,148 @@ TclDoubleDigits(char * strPtr,	/* Buffer in which to store the result, must
 /*
  *----------------------------------------------------------------------
  *
+ * AbsoluteValue --
+ *
+ *	Splits a 'double' into its absolute value and sign.
+ *
+ * Results:
+ *	Returns the absolute value.
+ *
+ * Side effects:
+ *	Stores the signum in '*signum'.
+ *
+ *----------------------------------------------------------------------
+ */
+
+static double
+AbsoluteValue (double v,	/* Number to split */
+	       int* signum)	/* (Output) Sign of the number 1=-, 0=+ */
+{
+    /*
+     * Take the absolute value of the number, and report the number's sign.
+     * Take special steps to preserve signed zeroes in IEEE floating point.
+     * (We can't use fpclassify, because that's a C9x feature and we still
+     * have to build on C89 compilers.)
+     */
+
+#ifndef IEEE_FLOATING_POINT
+    if (v >= 0.0) {
+	*signum = 0;
+    } else {
+	*signum = 1;
+	v = -v;
+    }
+#else
+    union {
+	Tcl_WideUInt iv;
+	double dv;
+    } bitwhack;
+    bitwhack.dv = v;
+    if (bitwhack.iv & ((Tcl_WideUInt) 1 << 63)) {
+	*signum = 1;
+	bitwhack.iv &= ~((Tcl_WideUInt) 1 << 63);
+	v = bitwhack.dv;
+    } else {
+	*signum = 0;
+    }
+#endif
+    return v;
+}
+
+/*
+ *----------------------------------------------------------------------
+ *
+ * GetIntegerTimesPower --
+ *
+ *	Converts a floating point number to an exact integer times a
+ *	power of the floating point radix.
+ *
+ * Results:
+ *	Returns 1 if it converted the smallest significand, 0 otherwise.
+ *
+ * Side effects:
+ *	Initializes the integer value (does not just assign it),
+ *	and stores the exponent.
+ *
+ *----------------------------------------------------------------------
+ */
+
+static int
+GetIntegerTimesPower(double v,	/* Value to convert */
+		     mp_int* rPtr,
+				/* (Output) Integer value */
+		     int* ePtr)	/* (Output) Power of FLT_RADIX by which
+				 * r must be multiplied to yield v*/
+{
+
+    double a;
+    double f;
+    int e;
+    int i;
+    int n;
+
+    /*
+     * Develop f and e such that v = f * FLT_RADIX**e, with
+     * 1.0/FLT_RADIX <= f < 1.
+     */
+
+    f = frexp(v, &e);
+#if FLT_RADIX > 2
+    n = e % log2FLT_RADIX;
+    if (n > 0) {
+	n -= log2FLT_RADIX;
+	e += 1;
+	f *= ldexp(1.0, n);
+    }
+    e = (e - n) / log2FLT_RADIX;
+#endif
+    if (f == 1.0) {
+	f = 1.0 / FLT_RADIX;
+	e += 1;
+    }
+
+    /*
+     * If the original number was denormalized, adjust e and f to be denormal
+     * as well.
+     */
+
+    if (e < DBL_MIN_EXP) {
+	n = mantBits + (e - DBL_MIN_EXP)*log2FLT_RADIX;
+	f = ldexp(f, (e - DBL_MIN_EXP)*log2FLT_RADIX);
+	e = DBL_MIN_EXP;
+	n = (n + DIGIT_BIT - 1) / DIGIT_BIT;
+    } else {
+	n = mantDIGIT;
+    }
+
+    /*
+     * Now extract the base-2**DIGIT_BIT digits of f into a multi-precision
+     * integer r. Preserve the invariant v = r * 2**rfac2 * FLT_RADIX**e by
+     * adjusting e.
+     */
+
+    a = f;
+    n = mantDIGIT;
+    mp_init_size(rPtr, n);
+    rPtr->used = n;
+    rPtr->sign = MP_ZPOS;
+    i = (mantBits % DIGIT_BIT);
+    if (i == 0) {
+	i = DIGIT_BIT;
+    }
+    while (n > 0) {
+	a *= ldexp(1.0, i);
+	i = DIGIT_BIT;
+	rPtr->dp[--n] = (mp_digit) a;
+	a -= (mp_digit) a;
+    }
+    *ePtr = e - DBL_MANT_DIG;
+    return (f == 1.0 / FLT_RADIX);
+}
+
+/*
+ *----------------------------------------------------------------------
+ *
  * TclInitDoubleConversion --
  *
  *	Initializes constants that are needed for conversions to and from
@@ -1138,16 +2113,43 @@ TclInitDoubleConversion(void)
 {
     int i;
     int x;
+    Tcl_WideUInt u;
     double d;
 
-    if (frexp((double) FLT_RADIX, &log2FLT_RADIX) != 0.5) {
-	Tcl_Panic("This code doesn't work on a decimal machine!");
+    /*
+     * Initialize table of powers of 10 expressed as wide integers.
+     */
+
+    maxpow10_wide =
+	(int) floor(sizeof (Tcl_WideUInt) * CHAR_BIT * log (2.) / log (10.));
+    pow10_wide = (Tcl_WideUInt*) Tcl_Alloc ((maxpow10_wide + 1)
+					    * sizeof (Tcl_WideUInt));
+    u = 1;
+    for (i = 0; i < maxpow10_wide; ++i) {
+	pow10_wide[i] = u;
+	u *= 10;
+    }
+    pow10_wide[i] = u;
+
+    /*
+     * Determine how many bits of precision a double has, and how many
+     * decimal digits that represents.
+     */
+
+    if ( frexp( (double) FLT_RADIX, &log2FLT_RADIX ) != 0.5 ) {
+	Tcl_Panic( "This code doesn't work on a decimal machine!" );
     }
     --log2FLT_RADIX;
     mantBits = DBL_MANT_DIG * log2FLT_RADIX;
     d = 1.0;
-    x = (int) (DBL_MANT_DIG * log((double) FLT_RADIX) / log(5.0));
-    if (x < MAXPOW) {
+
+    /* 
+     * Initialize a table of powers of ten that can be exactly represented
+     * in a double.
+     */
+
+    x = (int) (DBL_MANT_DIG * log((double) FLT_RADIX) / log( 5.0 ));
+    if ( x < MAXPOW ) {
 	mmaxpow = x;
     } else {
 	mmaxpow = MAXPOW;
@@ -1156,19 +2158,32 @@ TclInitDoubleConversion(void)
 	pow10[i] = d;
 	d *= 10.0;
     }
-    for (i=0 ; i<9 ; ++i) {
-	mp_init(pow5 + i);
+
+    /* Initialize a table of large powers of five. */
+
+    for ( i = 0; i < 9; ++i ) {
+	mp_init( pow5 + i );
     }
-    mp_set(pow5, 5);
-    for (i=0 ; i<8 ; ++i) {
-	mp_sqr(pow5+i, pow5+i+1);
+    mp_set( pow5, 5 );
+    for ( i = 0; i < 8; ++i ) {
+	mp_sqr( pow5+i, pow5+i+1 );
     }
-    tiny = SafeLdExp(1.0, DBL_MIN_EXP * log2FLT_RADIX - mantBits);
-    maxDigits = (int)
-	    ((DBL_MAX_EXP * log((double) FLT_RADIX) + log(10.)/2) / log(10.));
-    minDigits = (int)
-	    floor((DBL_MIN_EXP-DBL_MANT_DIG)*log((double)FLT_RADIX)/log(10.));
-    mantDIGIT = (mantBits + DIGIT_BIT - 1) / DIGIT_BIT;
+
+    /* 
+     * Determine the number of decimal digits to the left and right of the
+     * decimal point in the largest and smallest double, the smallest double
+     * that differs from zero, and the number of mp_digits needed to represent 
+     * the significand of a double.
+     */
+
+    tiny = SafeLdExp( 1.0, DBL_MIN_EXP * log2FLT_RADIX - mantBits );
+    maxDigits = (int) ((DBL_MAX_EXP * log((double) FLT_RADIX)
+			+ 0.5 * log(10.))
+		       / log( 10. ));
+    minDigits = (int) floor ( ( DBL_MIN_EXP - DBL_MANT_DIG )
+			      * log( (double) FLT_RADIX ) / log( 10. ) );
+    mantDIGIT = ( mantBits + DIGIT_BIT - 1 ) / DIGIT_BIT;
+    log10_DIGIT_MAX = (int) floor (DIGIT_BIT * log(2.) / log (10.));
 }
 
 /*
@@ -1191,9 +2206,62 @@ void
 TclFinalizeDoubleConversion()
 {
     int i;
-    for (i=0 ; i<9 ; ++i) {
-	mp_clear(pow5 + i);
+    Tcl_Free ((char*)pow10_wide);
+    for ( i = 0; i < 9; ++i ) {
+	mp_clear( pow5 + i );
+    }
+}
+
+/*
+ *----------------------------------------------------------------------
+ *
+ * TclInitBignumFromDouble --
+ *
+ *	Extracts the integer part of a double and converts it to
+ *	an arbitrary precision integer.
+ *
+ * Results:
+ *	None.
+ *
+ * Side effects:
+ *	Initializes the bignum supplied, and stores the converted number
+ *      in it.
+ *
+ *----------------------------------------------------------------------
+ */
+
+int
+TclInitBignumFromDouble(Tcl_Interp *interp,	/* For error message */
+			double d,	/* Number to convert */
+			mp_int* b)	/* Place to store the result */
+{
+    double fract;
+    int expt;
+
+    /* Infinite values can't convert to bignum */
+    if (TclIsInfinite(d)) {
+	if (interp != NULL) {
+	    char *s = "integer value too large to represent";
+	    Tcl_SetObjResult(interp, Tcl_NewStringObj(s, -1));
+	    Tcl_SetErrorCode(interp, "ARITH", "IOVERFLOW", s, NULL);
+	}
+	return TCL_ERROR;
+    }
+    fract = frexp(d,&expt);
+    if (expt <= 0) {
+	mp_init(b);
+	mp_zero(b);
+    } else {
+	Tcl_WideInt w = (Tcl_WideInt) ldexp(fract, mantBits);
+	int shift = expt - mantBits;
+	TclBNInitBignumFromWideInt(b, w);
+	if (shift < 0) {
+	    mp_div_2d(b, -shift, b, NULL);
+	} else if (shift > 0) {
+	    mp_mul_2d(b, shift, b);
+	}
     }
+    return TCL_OK;
 }
 
 /*
@@ -1228,7 +2296,11 @@ TclBignumToDouble(mp_int *a)	/* Integer to convert. */
     bits = mp_count_bits(a);
     if (bits > DBL_MAX_EXP*log2FLT_RADIX) {
 	errno = ERANGE;
-	return HUGE_VAL;
+	if (a->sign == MP_ZPOS) {
+	    return HUGE_VAL;
+	} else {
+	    return -HUGE_VAL;
+	}
     }
     shift = mantBits + 1 - bits;
     mp_init(&b);
@@ -1268,6 +2340,210 @@ TclBignumToDouble(mp_int *a)	/* Integer to convert. */
 	return -r;
     }
 }
+
+double
+TclCeil(mp_int *a)	/* Integer to convert. */
+{
+    double r = 0.0;
+    mp_int b;
+
+    mp_init(&b);
+    if (mp_cmp_d(a, 0) == MP_LT) {
+	mp_neg(a, &b);
+	r = -TclFloor(&b);
+    } else {
+	int bits = mp_count_bits(a);
+
+	if (bits > DBL_MAX_EXP*log2FLT_RADIX) {
+	    r = HUGE_VAL;
+	} else {
+	    int i, exact = 1, shift = mantBits - bits;
+
+	    if (shift > 0) {
+		mp_mul_2d(a, shift, &b);
+	    } else if (shift < 0) {
+		mp_int d;
+		mp_init(&d);
+		mp_div_2d(a, -shift, &b, &d);
+		exact = mp_iszero(&d);
+		mp_clear(&d);
+	    } else {
+		mp_copy(a, &b);
+	    }
+	    if (!exact) {
+		mp_add_d(&b, 1, &b);
+	    }
+	    for (i=b.used-1 ; i>=0 ; --i) {
+		r = ldexp(r, DIGIT_BIT) + b.dp[i];
+	    }
+	    r = ldexp(r, bits - mantBits);
+	}
+    }
+    mp_clear(&b);
+    return r;
+}
+
+double
+TclFloor(mp_int *a)	/* Integer to convert. */
+{
+    double r = 0.0;
+    mp_int b;
+
+    mp_init(&b);
+    if (mp_cmp_d(a, 0) == MP_LT) {
+	mp_neg(a, &b);
+	r = -TclCeil(&b);
+    } else {
+	int bits = mp_count_bits(a);
+
+	if (bits > DBL_MAX_EXP*log2FLT_RADIX) {
+	    r = DBL_MAX;
+	} else {
+	    int i, shift = mantBits - bits;
+
+	    if (shift > 0) {
+		mp_mul_2d(a, shift, &b);
+	    } else if (shift < 0) {
+		mp_div_2d(a, -shift, &b, NULL);
+	    } else {
+		mp_copy(a, &b);
+	    }
+	    for (i=b.used-1 ; i>=0 ; --i) {
+		r = ldexp(r, DIGIT_BIT) + b.dp[i];
+	    }
+	    r = ldexp(r, bits - mantBits);
+	}
+    }
+    mp_clear(&b);
+    return r;
+}
+
+/*
+ *----------------------------------------------------------------------
+ *
+ * BignumToBiasedFrExp --
+ *
+ *	Convert an arbitrary-precision integer to a native floating
+ *	point number in the range [0.5,1) times a power of two.
+ *	NOTE: Intentionally converts to a number that's a few
+ *	ulp too small, so that RefineApproximation will not overflow
+ *	near the high end of the machine's arithmetic range.
+ *
+ * Results:
+ *	Returns the converted number.
+ *
+ * Side effects:
+ *	Stores the exponent of two in 'machexp'.
+ *
+ *----------------------------------------------------------------------
+ */
+
+static double
+BignumToBiasedFrExp( mp_int* a,
+				/* Integer to convert */
+		     int* machexp )
+				/* Power of two */
+{
+    mp_int b;
+    int bits;
+    int shift;
+    int i;
+    double r;
+
+    /* Determine how many bits we need, and extract that many from 
+     * the input. Round to nearest unit in the last place. */
+
+    bits = mp_count_bits( a );
+    shift = mantBits - 2 - bits;
+    mp_init( &b );
+    if ( shift > 0 ) {
+	mp_mul_2d( a, shift, &b );
+    } else if ( shift < 0 ) {
+	mp_div_2d( a, -shift, &b, NULL );
+    } else {
+	mp_copy( a, &b );
+    }
+
+    /* Accumulate the result, one mp_digit at a time */
+
+    r = 0.0;
+    for ( i = b.used-1; i >= 0; --i ) {
+	r = ldexp( r, DIGIT_BIT ) + b.dp[i];
+    }
+    mp_clear( &b );
+
+    /* Return the result with the appropriate sign. */
+
+    *machexp = bits - mantBits + 2;
+    if ( a->sign == MP_ZPOS ) {
+	return r;
+    } else {
+	return -r;
+    }
+}
+
+/*
+ *----------------------------------------------------------------------
+ *
+ * Pow10TimesFrExp --
+ *
+ *	Multiply a power of ten by a number expressed as fraction and
+ *	exponent.
+ *
+ * Results:
+ *	Returns the significand of the result.
+ *
+ * Side effects:
+ *	Overwrites the 'machexp' parameter with the exponent of the
+ *	result.
+ *
+ * Assumes that 'exponent' is such that 10**exponent would be a double,
+ * even though 'fraction*10**(machexp+exponent)' might overflow.
+ *
+ *----------------------------------------------------------------------
+ */
+
+static double
+Pow10TimesFrExp( int exponent, 	/* Power of 10 to multiply by */
+		 double fraction,
+				/* Significand of multiplicand */
+		 int* machexp )	/* On input, exponent of multiplicand.
+				 * On output, exponent of result. */
+{
+    int i, j;
+    int expt = *machexp;
+    double retval = fraction;
+
+    if ( exponent > 0 ) {
+
+	/* Multiply by 10**exponent */
+	
+	retval = frexp( retval * pow10[ exponent & 0xf ], &j );
+	expt += j;
+	for ( i = 4; i < 9; ++i ) {
+	    if ( exponent & (1<<i) ) {
+		retval = frexp( retval * pow_10_2_n[ i ], &j );
+		expt += j;
+	    }
+	}
+    } else if ( exponent < 0 ) {
+
+	/* Divide by 10**-exponent */
+
+	retval = frexp( retval / pow10[ (-exponent) & 0xf ], &j );
+	expt += j;
+	for ( i = 4; i < 9; ++i ) {
+	    if ( (-exponent) & (1<<i) ) {
+		retval = frexp( retval / pow_10_2_n[ i ], &j );
+		expt += j;
+	    }
+	}
+    }
+
+    *machexp = expt;
+    return retval;
+    
+}
 
 /*
  *----------------------------------------------------------------------