2010-11-29 Kevin B. Kenny <kennykb@acm.org>

	* generic/tclInt.decls:
	* generic/tclInt.h:
	* generic/tclStrToD.c:
	* generic/tclTest.c:
	* generic/tclTomMath.decls:
	* generic/tclUtil.c:
	* tests/util.test:
	* unix/Makefile.in:
	* win/Makefile.in:
	* win/makefile.vc: Rewrite of Tcl_PrintDouble and TclDoubleDigits
	that (a) fixes a severe performance problem with floating point
	shimmering reported by Karl Lehenbauer, (b) allows TclDoubleDigits
	to generate the digit strings for 'e' and 'f' format, so that it
	can be used for tcl_precision != 0 (and possibly later for [format]),
	(c) fixes [Bug 3120139] by making TclPrintDouble inherently
	locale-independent, (d) adds test cases to util.test for
	correct rounding in difficult cases of TclDoubleDigits where fixed-
	precision results are requested. (e) adds test cases to util.test for
	the controversial aspects of [Bug 3105247]. As a side effect, two
	more modules from libtommath (bn_mp_set_int.c and bn_mp_init_set_int.c)
	are brought into the build, since the new code uses them.

14 files changed, 3559 insertions, 410 deletions

diff --git a/ChangeLog b/ChangeLog
index 2af363c..a0acecc 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,31 @@
+2010-11-29  Kevin B. Kenny  <kennykb@acm.org>
+
+	* generic/tclInt.decls:
+	* generic/tclInt.h:
+	* generic/tclStrToD.c:
+	* generic/tclTest.c:
+	* generic/tclTomMath.decls:
+	* generic/tclUtil.c:
+	* tests/util.test:
+	* unix/Makefile.in:
+	* win/Makefile.in:
+	* win/makefile.vc: Rewrite of Tcl_PrintDouble and TclDoubleDigits
+	that (a) fixes a severe performance problem with floating point
+	shimmering reported by Karl Lehenbauer, (b) allows TclDoubleDigits
+	to generate the digit strings for 'e' and 'f' format, so that it
+	can be used for tcl_precision != 0 (and possibly later for [format]),
+	(c) fixes [Bug 3120139] by making TclPrintDouble inherently
+	locale-independent, (d) adds test cases to util.test for
+	correct rounding in difficult cases of TclDoubleDigits where fixed-
+	precision results are requested. (e) adds test cases to util.test for
+	the controversial aspects of [Bug 3105247]. As a side effect, two
+	more modules from libtommath (bn_mp_set_int.c and bn_mp_init_set_int.c)
+	are brought into the build, since the new code uses them.
+
+	* generic/tclIntDecls.h:
+	* generic/tclStubInit.c:
+	* generic/tclTomMathDecls.h:	Regenerated.
+
 2010-11-24  Donal K. Fellows  <dkf@users.sf.net>
 
 	* tests/chanio.test, tests/iogt.test, tests/ioTrans.test: Convert more
diff --git a/generic/tclInt.decls b/generic/tclInt.decls
index 1a114f1..75780d5 100644
--- a/generic/tclInt.decls
+++ b/generic/tclInt.decls
@@ -13,7 +13,7 @@
 # See the file "license.terms" for information on usage and redistribution
 # of this file, and for a DISCLAIMER OF ALL WARRANTIES.
 #
-# RCS: @(#) $Id: tclInt.decls,v 1.150 2010/10/02 00:23:44 hobbs Exp $
+# RCS: @(#) $Id: tclInt.decls,v 1.151 2010/11/28 23:20:10 kennykb Exp $
 
 library tcl
 
@@ -996,6 +996,11 @@ declare 248 {
     int TclCopyChannel(Tcl_Interp *interp, Tcl_Channel inChan,
 	    Tcl_Channel outChan, Tcl_WideInt toRead, Tcl_Obj *cmdPtr)
 }
+
+declare 249 {
+    char* TclDoubleDigits(double dv, int ndigits, int flags,
+			  int* decpt, int* signum, char** endPtr)
+}
 
 ##############################################################################
 
diff --git a/generic/tclInt.h b/generic/tclInt.h
index 218c40b..8c5d863 100644
--- a/generic/tclInt.h
+++ b/generic/tclInt.h
@@ -15,7 +15,7 @@
  * See the file "license.terms" for information on usage and redistribution of
  * this file, and for a DISCLAIMER OF ALL WARRANTIES.
  *
- * RCS: @(#) $Id: tclInt.h,v 1.486 2010/11/15 21:34:54 andreas_kupries Exp $
+ * RCS: @(#) $Id: tclInt.h,v 1.487 2010/11/28 23:20:11 kennykb Exp $
  */
 
 #ifndef _TCLINT
@@ -2795,6 +2795,32 @@ struct Tcl_LoadHandle_ {
 				/* Procedure that unloads a loaded module */
 };
 
+/* Flags for conversion of doubles to digit strings */
+
+#define TCL_DD_SHORTEST 		0x4
+				/* Use the shortest possible string */
+#define TCL_DD_STEELE   		0x5
+				/* Use the original Steele&White algorithm */
+#define TCL_DD_E_FORMAT 		0x2
+				/* Use a fixed-length string of digits,
+				 * suitable for E format*/
+#define TCL_DD_F_FORMAT 		0x3
+				/* Use a fixed number of digits after the
+				 * decimal point, suitable for F format */
+
+#define TCL_DD_SHORTEN_FLAG 		0x4
+				/* Allow return of a shorter digit string
+				 * if it converts losslessly */
+#define TCL_DD_NO_QUICK 		0x8
+				/* Debug flag: forbid quick FP conversion */
+
+#define TCL_DD_CONVERSION_TYPE_MASK	0x3
+				/* Mask to isolate the conversion type */
+#define TCL_DD_STEELE0 			0x1
+				/* 'Steele&White' after masking */
+#define TCL_DD_SHORTEST0		0x0
+				/* 'Shortest possible' after masking */
+
 /*
  *----------------------------------------------------------------
  * Procedures shared among Tcl modules but not used by the outside world:
@@ -2843,7 +2869,6 @@ MODULE_SCOPE void	TclContinuationsEnterDerived(Tcl_Obj *objPtr,
 MODULE_SCOPE ContLineLoc *TclContinuationsGet(Tcl_Obj *objPtr);
 MODULE_SCOPE void	TclContinuationsCopy(Tcl_Obj *objPtr,
 			    Tcl_Obj *originObjPtr);
-MODULE_SCOPE int	TclDoubleDigits(char *buf, double value, int *signum);
 MODULE_SCOPE void	TclDeleteNamespaceVars(Namespace *nsPtr);
 /* TIP #280 - Modified token based evulation, with line information. */
 MODULE_SCOPE int	TclEvalEx(Tcl_Interp *interp, const char *script,
diff --git a/generic/tclIntDecls.h b/generic/tclIntDecls.h
index 259229b..ae4046d 100644
--- a/generic/tclIntDecls.h
+++ b/generic/tclIntDecls.h
@@ -11,7 +11,7 @@
  * See the file "license.terms" for information on usage and redistribution
  * of this file, and for a DISCLAIMER OF ALL WARRANTIES.
  *
- * RCS: @(#) $Id: tclIntDecls.h,v 1.144 2010/10/02 00:23:45 hobbs Exp $
+ * RCS: @(#) $Id: tclIntDecls.h,v 1.145 2010/11/28 23:20:11 kennykb Exp $
  */
 
 #ifndef _TCLINTDECLS
@@ -596,6 +596,9 @@ EXTERN void		TclResetRewriteEnsemble(Tcl_Interp *interp,
 EXTERN int		TclCopyChannel(Tcl_Interp *interp,
 				Tcl_Channel inChan, Tcl_Channel outChan,
 				Tcl_WideInt toRead, Tcl_Obj *cmdPtr);
+/* 249 */
+EXTERN char*		TclDoubleDigits(double dv, int ndigits, int flags,
+				int*decpt, int*signum, char**endPtr);
 
 typedef struct TclIntStubs {
     int magic;
@@ -850,6 +853,7 @@ typedef struct TclIntStubs {
     int (*tclInitRewriteEnsemble) (Tcl_Interp *interp, int numRemoved, int numInserted, Tcl_Obj *const *objv); /* 246 */
     void (*tclResetRewriteEnsemble) (Tcl_Interp *interp, int isRootEnsemble); /* 247 */
     int (*tclCopyChannel) (Tcl_Interp *interp, Tcl_Channel inChan, Tcl_Channel outChan, Tcl_WideInt toRead, Tcl_Obj *cmdPtr); /* 248 */
+    char* (*tclDoubleDigits) (double dv, int ndigits, int flags, int*decpt, int*signum, char**endPtr); /* 249 */
 } TclIntStubs;
 
 #ifdef __cplusplus
@@ -1269,6 +1273,8 @@ extern const TclIntStubs *tclIntStubsPtr;
 	(tclIntStubsPtr->tclResetRewriteEnsemble) /* 247 */
 #define TclCopyChannel \
 	(tclIntStubsPtr->tclCopyChannel) /* 248 */
+#define TclDoubleDigits \
+	(tclIntStubsPtr->tclDoubleDigits) /* 249 */
 
 #endif /* defined(USE_TCL_STUBS) */
 
diff --git a/generic/tclStrToD.c b/generic/tclStrToD.c
index d0a5345..8171df0 100755
--- a/generic/tclStrToD.c
+++ b/generic/tclStrToD.c
@@ -14,7 +14,7 @@
  * 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.46 2010/05/21 12:43:29 nijtmans Exp $
+ * RCS: @(#) $Id: tclStrToD.c,v 1.47 2010/11/28 23:20:11 kennykb Exp $
  *
  *----------------------------------------------------------------------
  */
@@ -90,6 +90,66 @@ typedef unsigned int fpu_control_t __attribute__ ((__mode__ (__HI__)));
  * runtime).
  */
 
+/* Magic constants */
+
+#define LOG10_2 0.3010299956639812
+#define TWO_OVER_3LOG10 0.28952965460216784
+#define LOG10_3HALVES_PLUS_FUDGE 0.1760912590558
+
+/* Definitions of the parts of an IEEE754-format floating point number */
+
+#define SIGN_BIT 		0x80000000
+				/* Mask for the sign bit in the first
+				 * word of a double */
+#define EXP_MASK	     	0x7ff00000
+				/* Mask for the exponent field in the
+				 * first word of a double */
+#define EXP_SHIFT		20
+				/* Shift count to make the exponent an 
+				 * integer */
+#define HIDDEN_BIT		(((Tcl_WideUInt) 0x00100000) << 32)
+				/* Hidden 1 bit for the significand */
+#define HI_ORDER_SIG_MASK	0x000fffff
+				/* Mask for the high-order part of the
+				 * significand in the first word of a
+				 * double */
+#define SIG_MASK		(((Tcl_WideUInt) HI_ORDER_SIG_MASK << 32) \
+				 | 0xffffffff)
+				/* Mask for the 52-bit significand. */
+#define FP_PRECISION		53
+				/* Number of bits of significand plus the
+				 * hidden bit */
+#define EXPONENT_BIAS		0x3ff
+				/* Bias of the exponent 0 */
+
+/* Derived quantities */
+
+#define TEN_PMAX		22
+				/* floor(FP_PRECISION*log(2)/log(5)) */
+#define QUICK_MAX		14
+				/* floor((FP_PRECISION-1)*log(2)/log(10)) - 1 */
+#define BLETCH			0x10
+				/* Highest power of two that is greater than
+				 * DBL_MAX_10_EXP, divided by 16 */
+#define DIGIT_GROUP		8
+				/* floor(DIGIT_BIT*log(2)/log(10)) */
+
+/* Union used to dismantle floating point numbers. */
+
+typedef union Double {
+    struct {
+#ifdef WORDS_BIGENDIAN
+	int word0;
+	int word1;
+#else
+	int word1;
+	int word0;
+#endif
+    } w;
+    double d;
+    Tcl_WideUInt q;
+} Double;
+
 static int maxpow10_wide;	/* The powers of ten that can be represented
 				 * exactly as wide integers. */
 static Tcl_WideUInt *pow10_wide;
@@ -123,6 +183,7 @@ static const double pow_10_2_n[] = {	/* Inexact higher powers of ten. */
     1.0e+128,
     1.0e+256
 };
+
 static int n770_fp;		/* Flag is 1 on Nokia N770 floating point.
 				 * Nokia's floating point has the words
 				 * reversed: if big-endian is 7654 3210,
@@ -131,27 +192,161 @@ static int n770_fp;		/* Flag is 1 on Nokia N770 floating point.
 				 * little-endian within the 32-bit words
 				 * but big-endian between them. */
 
+/* Table of powers of 5 that are small enough to fit in an mp_digit. */
+
+static const mp_digit dpow5[13] = {
+               1,              5,             25,            125,
+             625,           3125,          15625,          78125,
+          390625,        1953125,        9765625,       48828125,
+       244140625
+};
+
+/* Table of powers: pow5_13[n] = 5**(13*2**(n+1)) */
+static mp_int pow5_13[5];	/* Table of powers: 5**13, 5**26, 5**52,
+				 * 5**104, 5**208 */
+static const double tens[] = {
+    1e00, 1e01, 1e02, 1e03, 1e04, 1e05, 1e06, 1e07, 1e08, 1e09,
+    1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
+    1e20, 1e21, 1e22
+};
+
+static const int itens [] = {
+    1,
+    10,
+    100,
+    1000,
+    10000,
+    100000,
+    1000000,
+    10000000,
+    100000000
+};
+
+static const Tcl_WideUInt wtens[] = {
+    1, 10, 100, 1000, 10000, 100000, 1000000,
+    (Tcl_WideUInt) 1000000*10, 		(Tcl_WideUInt) 1000000*100,
+    (Tcl_WideUInt) 1000000*1000, 	(Tcl_WideUInt) 1000000*10000,
+    (Tcl_WideUInt) 1000000*100000, 	(Tcl_WideUInt) 1000000*1000000,
+    (Tcl_WideUInt) 1000000*1000000*10, 	(Tcl_WideUInt) 1000000*1000000*100,
+    (Tcl_WideUInt) 1000000*1000000*1000,(Tcl_WideUInt) 1000000*1000000*10000
+    
+};
+
+static const double bigtens[] = {
+    1e016, 1e032, 1e064, 1e128, 1e256
+};
+#define N_BIGTENS 5
+
+static const int log2pow5[27] = {
+    01,  3,  5,  7, 10, 12, 14, 17, 19, 21,
+    24, 26, 28, 31, 33, 35, 38, 40, 42, 45,
+    47, 49, 52, 54, 56, 59, 61
+};
+#define N_LOG2POW5 27
+
+static const Tcl_WideUInt wuipow5[27] = {
+    (Tcl_WideUInt) 1,		/* 5**0 */
+    (Tcl_WideUInt) 5,
+    (Tcl_WideUInt) 25,
+    (Tcl_WideUInt) 125,
+    (Tcl_WideUInt) 625,
+    (Tcl_WideUInt) 3125,	/* 5**5 */
+    (Tcl_WideUInt) 3125*5,
+    (Tcl_WideUInt) 3125*25,
+    (Tcl_WideUInt) 3125*125,
+    (Tcl_WideUInt) 3125*625,
+    (Tcl_WideUInt) 3125*3125,	/* 5**10 */
+    (Tcl_WideUInt) 3125*3125*5,
+    (Tcl_WideUInt) 3125*3125*25,
+    (Tcl_WideUInt) 3125*3125*125,
+    (Tcl_WideUInt) 3125*3125*625,
+    (Tcl_WideUInt) 3125*3125*3125, /* 5**15 */
+    (Tcl_WideUInt) 3125*3125*3125*5,
+    (Tcl_WideUInt) 3125*3125*3125*25,
+    (Tcl_WideUInt) 3125*3125*3125*125,
+    (Tcl_WideUInt) 3125*3125*3125*625,
+    (Tcl_WideUInt) 3125*3125*3125*3125,	/* 5**20 */
+    (Tcl_WideUInt) 3125*3125*3125*3125*5,
+    (Tcl_WideUInt) 3125*3125*3125*3125*25,
+    (Tcl_WideUInt) 3125*3125*3125*3125*125,
+    (Tcl_WideUInt) 3125*3125*3125*3125*625,
+    (Tcl_WideUInt) 3125*3125*3125*3125*3125,  /* 5**25 */
+    (Tcl_WideUInt) 3125*3125*3125*3125*3125*5 /* 5**26 */
+};
+
 /*
  * Static functions defined in this file.
  */
 
-static double		AbsoluteValue(double v, int *signum);
 static int		AccumulateDecimalDigit(unsigned, int, 
 			    Tcl_WideUInt *, mp_int *, int);
-static double		BignumToBiasedFrExp(const mp_int *big, int *machexp);
-static int		GetIntegerTimesPower(double v, mp_int *r, int *e);
 static double		MakeHighPrecisionDouble(int signum,
 			    mp_int *significand, int nSigDigs, int exponent);
 static double		MakeLowPrecisionDouble(int signum,
 			    Tcl_WideUInt significand, int nSigDigs,
 			    int exponent);
 static double		MakeNaN(int signum, Tcl_WideUInt tag);
-static Tcl_WideUInt	Nokia770Twiddle(Tcl_WideUInt w);
-static double		Pow10TimesFrExp(int exponent, double fraction,
-			    int *machexp);
 static double		RefineApproximation(double approx,
 			    mp_int *exactSignificand, int exponent);
+static void		MulPow5(mp_int*, unsigned, mp_int*);
+static int 		NormalizeRightward(Tcl_WideUInt*);
+static int		RequiredPrecision(Tcl_WideUInt);
+static void		DoubleToExpAndSig(double, Tcl_WideUInt*, int*, int*);
+static void		TakeAbsoluteValue(Double*, int*);
+static char*		FormatInfAndNaN(Double*, int*, char**);
+static char*		FormatZero(int*, char**);
+static int		ApproximateLog10(Tcl_WideUInt, int, int);
+static int		BetterLog10(double, int, int*);
+static void		ComputeScale(int, int, int*, int*, int*, int*);
+static void		SetPrecisionLimits(int, int, int*, int*, int*, int*);
+static char*		BumpUp(char*, char*, int*);
+static int		AdjustRange(double*, int);
+static char*		ShorteningQuickFormat(double, int, int, double, 
+			    char*, int*);
+static char*		StrictQuickFormat(double, int, int, double,
+			    char*, int*);
+static char*		QuickConversion(double, int, int, int, int, int, int,
+			    int*, char**);
+static void		CastOutPowersOf2(int*, int*, int*);
+static char*		ShorteningInt64Conversion(Double*, int, Tcl_WideUInt,
+			    int, int, int, int, int, int, int, int, int,
+			    int, int, int*, char**);
+static char*		StrictInt64Conversion(Double*, int, Tcl_WideUInt,
+			    int, int, int, int, int, int,
+			    int, int, int*, char**);
+static int		ShouldBankerRoundUpPowD(mp_int*, int, int);
+static int		ShouldBankerRoundUpToNextPowD(mp_int*, mp_int*,
+			    int, int, int, mp_int*);
+static char*		ShorteningBignumConversionPowD(Double* dPtr, 
+			    int convType, Tcl_WideUInt bw, int b2, int b5,
+			    int m2plus, int m2minus, int m5,
+			    int sd, int k, int len, 
+			    int ilim, int ilim1, int* decpt,
+			    char** endPtr);
+static char*		StrictBignumConversionPowD(Double* dPtr, int convType,
+			    Tcl_WideUInt bw, int b2, int b5,
+			    int sd, int k, int len, 
+			    int ilim, int ilim1, int* decpt,
+			    char** endPtr);
+static int		ShouldBankerRoundUp(mp_int*, mp_int*, int);
+static int		ShouldBankerRoundUpToNext(mp_int*, mp_int*, mp_int*,
+			    int, int, mp_int*);
+static char*		ShorteningBignumConversion(Double* dPtr, int convType,
+			    Tcl_WideUInt bw, int b2,
+			    int m2plus, int m2minus,
+			    int s2, int s5, int k, int len, 
+			    int ilim, int ilim1, int* decpt,
+			    char** endPtr);
+static char*		StrictBignumConversion(Double* dPtr, int convType,
+			    Tcl_WideUInt bw, int b2,
+			    int s2, int s5, int k, int len, 
+			    int ilim, int ilim1, int* decpt,
+			    char** endPtr);
+static double		BignumToBiasedFrExp(const mp_int *big, int *machexp);
+static double		Pow10TimesFrExp(int exponent, double fraction,
+			    int *machexp);
 static double		SafeLdExp(double fraction, int exponent);
+static Tcl_WideUInt	Nokia770Twiddle(Tcl_WideUInt w);
 
 /*
  *----------------------------------------------------------------------
@@ -1732,414 +1927,2362 @@ RefineApproximation(
 }
 
 /*
- *----------------------------------------------------------------------
+ *-----------------------------------------------------------------------------
  *
- * TclDoubleDigits --
+ * MultPow5 --
+ *
+ *	Multiply a bignum by a power of 5.
+ *
+ * Side effects:
+ *	Stores base*5**n in result
+ *
+ *-----------------------------------------------------------------------------
+ */
+
+inline static void
+MulPow5(mp_int* base, 		/* Number to multiply */
+	 unsigned n,		/* Power of 5 to multiply by */
+	 mp_int* result)	/* Place to store the result */
+{
+    mp_int* p = base;
+    int n13 = n / 13;
+    int r = n % 13;
+    if (r != 0) {
+	mp_mul_d(p, dpow5[r], result);
+	p = result;
+    }
+    r = 0;
+    while (n13 != 0) {
+	if (n13 & 1) {
+	    mp_mul(p, pow5_13+r, result);
+	    p = result;
+	}
+	n13 >>= 1;
+	++r;
+    }
+    if (p != result) {
+	mp_copy(p, result);
+    }
+}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * NormalizeRightward --
  *
- *	Converts a double to a string of digits.
+ *	Shifts a number rightward until it is odd (that is, until the
+ *	least significant bit is nonzero.
  *
  * Results:
- *	Returns the position of the character in the string after which the
- *	decimal point should appear. Since the string contains only
- *	significant digits, the position may be less than zero or greater than
- *	the length of the string.
+ *	Returns the number of bit positions by which the number was shifted.
  *
  * Side effects:
- *	Stores the digits in the given buffer and sets 'signum' according to
- *	the sign of the number.
+ *	Shifts the number in place; *wPtr is replaced by the shifted number.
  *
- *----------------------------------------------------------------------
+ *-----------------------------------------------------------------------------
+ */
 
+inline static int
+NormalizeRightward(Tcl_WideUInt* wPtr)
+				/* INOUT: Number to shift */
+{
+    int rv = 0;
+    Tcl_WideUInt w = *wPtr;
+    if (!(w & (Tcl_WideUInt) 0xffffffff)) {
+	w >>= 32; rv += 32;
+    } 
+    if (!(w & (Tcl_WideUInt) 0xffff)) {
+	w >>= 16; rv += 16;
+    }
+    if (!(w & (Tcl_WideUInt) 0xff)) {
+	w >>= 8; rv += 8;
+    }
+    if (!(w & (Tcl_WideUInt) 0xf)) {
+	w >>= 4; rv += 4;
+    }
+    if (!(w & 0x3)) {
+	w >>= 2; rv += 2;
+    }
+    if (!(w & 0x1)) {
+	w >>= 1; ++rv;
+    }
+    *wPtr = w;
+    return rv;
+}
+
+/*
+ *-----------------------------------------------------------------------------0
+ *
+ * RequiredPrecision --
+ *
+ *	Determines the number of bits needed to hold an intger.
+ *
+ * Results:
+ *	Returns the position of the most significant bit (0 - 63).
+ *	Returns 0 if the number is zero.
+ *
+ *----------------------------------------------------------------------------
  */
 
-int
-TclDoubleDigits(
-    char *buffer,		/* 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. */
+static int
+RequiredPrecision(Tcl_WideUInt w)
+				/* Number to interrogate */
 {
-    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. */
-    mp_int mminus;		/* Scaled epsilon: (r - 2*mminus) == v(-)
-				 * where v(-) is the floating point
-				 * predecessor of v. */
-    mp_int temp;
-    int rfac2 = 0;		/* Powers of 2 and 5 by which large */
-    int rfac5 = 0;		/* integers should be scaled.	    */
-    int sfac2 = 0;
-    int sfac5 = 0;
-    int mplusfac2 = 0;
-    int mminusfac2 = 0;
-    char c;
-    int i, k, n;
+    int rv;
+    unsigned long wi;
+    if (w & ((Tcl_WideUInt) 0xffffffff << 32)) {
+	wi = w >> 32; rv = 32;
+    } else {
+	wi = w; rv = 0;
+    }
+    if (wi & 0xffff0000) {
+	wi >>= 16; rv += 16;
+    }
+    if (wi & 0xff00) {
+	wi >>= 8; rv += 8;
+    }
+    if (wi & 0xf0) {
+	wi >>= 4; rv += 4;
+    }
+    if (wi & 0xc) {
+	wi >>= 2; rv += 2;
+    }
+    if (wi & 0x2) {
+	wi >>= 1; ++rv;
+    }
+    if (wi & 0x1) {
+	++rv;
+    }
+    return rv;
+}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * DoubleToExpAndSig --
+ *
+ *	Separates a 'double' into exponent and significand.
+ *
+ * Side effects:
+ *	Stores the significand in '*significand' and the exponent in
+ *	'*expon' so that dv == significand * 2.0**expon, and significand
+ *	is odd.  Also stores the position of the leftmost 1-bit in 'significand'
+ *	in 'bits'.
+ *
+ *-----------------------------------------------------------------------------
+ */
+
+inline static void
+DoubleToExpAndSig(double dv,	/* Number to convert */
+		  Tcl_WideUInt* significand,
+				/* OUTPUT: Significand of the number */
+		  int* expon,	/* OUTPUT: Exponent to multiply the number by */
+		  int* bits)	/* OUTPUT: Number of significant bits */
+{
+    Double d;			/* Number being converted */
+    Tcl_WideUInt z;		/* Significand under construction */
+    int de;			/* Exponent of the number */
+    int k;			/* Bit count */
+
+    d.d = dv;
+
+    /* Extract exponent and significand */
+
+    de = (d.w.word0 & EXP_MASK) >> EXP_SHIFT;
+    z = d.q & SIG_MASK;
+    if (de != 0) {
+	z |= HIDDEN_BIT;
+	k = NormalizeRightward(&z);
+	*bits = FP_PRECISION - k;
+	*expon = k + (de - EXPONENT_BIAS) - (FP_PRECISION-1);
+    } else {
+	k = NormalizeRightward(&z);
+	*expon = k + (de - EXPONENT_BIAS) - (FP_PRECISION-1) + 1;
+	*bits = RequiredPrecision(z);
+    }
+    *significand = z;
+}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * TakeAbsoluteValue --
+ *
+ *	Takes the absolute value of a 'double' including 0, Inf and NaN
+ *
+ * Side effects:
+ *	The 'double' in *d is replaced with its absolute value. The
+ *	signum is stored in 'sign': 1 for negative, 0 for nonnegative.
+ *
+ *-----------------------------------------------------------------------------
+ */
+
+inline static void
+TakeAbsoluteValue(Double* d,	/* Number to replace with absolute value */
+		  int* sign)	/* Place to put the signum */
+{
+    if (d->w.word0 & SIGN_BIT) {
+	*sign = 1;
+	d->w.word0 &= ~SIGN_BIT;
+    } else {
+	*sign = 0;
+    }
+}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * FormatInfAndNaN --
+ *
+ *	Bailout for formatting infinities and Not-A-Number.
+ *
+ * Results:
+ *	Returns one of the strings 'Infinity' and 'NaN'.
+ *
+ * Side effects:
+ *	Stores 9999 in *decpt, and sets '*endPtr' to designate the
+ *	terminating NUL byte of the string if 'endPtr' is not NULL. 
+ *
+ * The string returned must be freed by the caller using 'ckfree'.
+ *
+ *-----------------------------------------------------------------------------
+ */
+
+inline static char*
+FormatInfAndNaN(Double* d,	/* Exceptional number to format */
+		int* decpt,	/* Decimal point to set to a bogus value */
+		char** endPtr)	/* Pointer to the end of the formatted data */
+{
+    char* retval;
+    *decpt = 9999;
+    if (!(d->w.word1) && !(d->w.word0 & HI_ORDER_SIG_MASK)) {
+	retval = ckalloc(9);
+	strcpy(retval, "Infinity");
+	if (endPtr) {
+	    *endPtr = retval + 8;
+	}
+    } else {
+	retval = ckalloc(4);
+	strcpy(retval, "NaN");
+	if (endPtr) {
+	    *endPtr = retval + 3;
+	}
+    }
+    return retval;
+}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * FormatZero --
+ *
+ *	Bailout to format a zero floating-point number.
+ *
+ * Results:
+ *	Returns the constant string "0"
+ *
+ * Side effects:
+ *	Stores 1 in '*decpt' and puts a pointer to the NUL byte terminating
+ *	the string in '*endPtr' if 'endPtr' is not NULL.
+ *
+ *-----------------------------------------------------------------------------
+ */
+
+inline static char*
+FormatZero(int* decpt,		/* Location of the decimal point */
+	   char** endPtr)	/* Pointer to the end of the formatted data */
+{
+    char* retval = ckalloc(2);
+    strcpy(retval, "0");
+    if (endPtr) {
+	*endPtr = retval+1;
+    }
+    *decpt = 0;
+    return retval;
+}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * ApproximateLog10 --
+ *
+ *	Computes a two-term Taylor series approximation to the common
+ *	log of a number, and computes the number's binary log.
+ *
+ * Results:
+ *	Return an approximation to floor(log10(bw*2**be)) that is either
+ *	exact or 1 too high.
+ *
+ *-----------------------------------------------------------------------------
+ */
+
+inline static int
+ApproximateLog10(Tcl_WideUInt bw,
+				/* Integer significand of the number */
+		 int be,	/* Power of two to scale bw */
+		 int bbits)	/* Number of bits of precision in bw */
+{
+    int i;			/* Log base 2 of the number */
+    int k;			/* Floor(Log base 10 of the number) */
+    double ds;			/* Mantissa of the number */
+    Double d2;
 
     /*
-     * Split the number into absolute value and signum.
+     * Compute i and d2 such that d = d2*2**i, and 1 < d2 < 2.
+     * Compute an approximation to log10(d), 
+     *   log10(d) ~ log10(2) * i + log10(1.5) 
+     *            + (significand-1.5)/(1.5 * log(10))
      */
 
-    v = AbsoluteValue(v, signum);
+    d2.q = bw << (FP_PRECISION - bbits) & SIG_MASK;
+    d2.w.word0 |= (EXPONENT_BIAS) << EXP_SHIFT;
+    i = be + bbits - 1;
+    ds = (d2.d - 1.5) * TWO_OVER_3LOG10
+	+ LOG10_3HALVES_PLUS_FUDGE
+	+ LOG10_2 * i;
+    k = (int) ds;
+    if (k > ds) {
+	--k;
+    }
+    return k;
+}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * BetterLog10 --
+ *
+ *	Improves the result of ApproximateLog10 for numbers in the range
+ *	1 .. 10**(TEN_PMAX)-1
+ *
+ * Side effects:
+ *	Sets k_check to 0 if the new result is known to be exact, and to
+ *	1 if it may still be one too high.
+ *
+ * Results:
+ *	Returns the improved approximation to log10(d)
+ *
+ *-----------------------------------------------------------------------------
+ */
 
-    /*
-     * Handle zero specially.
+inline static int
+BetterLog10(double d,		/* Original number to format */
+	  int k,		/* Characteristic(Log base 10) of the number */
+	  int* k_check)		/* Flag == 1 if k is inexact */
+{
+    /* 
+     * Performance hack. If k is in the range 0..TEN_PMAX, then we can
+     * use a powers-of-ten table to check it.
      */
+    if (k >= 0 && k <= TEN_PMAX) {
+	if (d < tens[k]) {
+	    k--;
+	}
+	*k_check = 0;
+    } else {
+	*k_check = 1;
+    }
+    return k;
+}
+
+/* 
+ *-----------------------------------------------------------------------------
+ *
+ * ComputeScale --
+ *
+ *	Prepares to format a floating-point number as decimal.
+ *
+ * Parameters:
+ *	floor(log10*x) is k (or possibly k-1).  floor(log2(x) is i.
+ *	The significand of x requires bbits bits to represent.
+ *
+ * Results:
+ *	Determines integers b2, b5, s2, s5 so that sig*2**b2*5**b5/2**s2*2**s5
+ *	exactly represents the value of the x/10**k. This value will lie
+ *	in the range [1 .. 10), and allows for computing successive digits
+ *	by multiplying sig%10 by 10.
+ *
+ *-----------------------------------------------------------------------------
+ */
 
-    if (v == 0.0) {
-	*buffer++ = '0';
-	*buffer++ = '\0';
-	return 1;
+inline static void
+ComputeScale(int be,		/* Exponent part of number: d = bw * 2**be */
+	     int k,		/* Characteristic of log10(number) */
+	     int* b2,		/* OUTPUT: Power of 2 in the numerator */
+	     int* b5,		/* OUTPUT: Power of 5 in the numerator */
+	     int* s2,		/* OUTPUT: Power of 2 in the denominator */
+	     int* s5)		/* OUTPUT: Power of 5 in the denominator */
+{
+
+    /* 
+     * Scale numerator and denominator powers of 2 so that the
+     * input binary number is the ratio of integers
+     */
+    if (be <= 0) {
+	*b2 = 0;
+	*s2 = -be;
+    } else {
+	*b2 = be;
+	*s2 = 0;
     }
 
+    /* 
+     * Scale numerator and denominator so that the output decimal number
+     * is the ratio of integers
+     */
+    if (k >= 0) {
+	*b5 = 0;
+	*s5 = k;
+	*s2 += k;
+    } else {
+	*b2 -= k;
+	*b5 = -k;
+	*s5 = 0;
+    }
+}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * SetPrecisionLimits --
+ *
+ *	Determines how many digits of significance should be computed
+ *	(and, hence, how much memory need be allocated) for formatting a
+ *	floating point number.
+ *
+ * Given that 'k' is floor(log10(x)):
+ * if 'shortest' format is used, there will be at most 18 digits in the result.
+ * if 'F' format is used, there will be at most 'ndigits' + k + 1 digits
+ * if 'E' format is used, there will be exactly 'ndigits' digits.
+ *
+ * Side effects:
+ *	Adjusts '*ndigitsPtr' to have a valid value.
+ *	Stores the maximum memory allocation needed in *iPtr.
+ *	Sets '*iLimPtr' to the limiting number of digits to convert if k
+ *	has been guessed correctly, and '*iLim1Ptr' to the limiting number
+ *	of digits to convert if k has been guessed to be one too high.
+ *
+ *-----------------------------------------------------------------------------
+ */
+
+inline static void
+SetPrecisionLimits(int convType,
+				/* Type of conversion:
+				 *   TCL_DD_SHORTEST
+				 *   TCL_DD_STEELE0
+				 *   TCL_DD_E_FMT
+				 *   TCL_DD_F_FMT */
+		   int k,	/* Floor(log10(number to convert)) */
+		   int* ndigitsPtr,
+				/* IN/OUT: Number of digits requested
+				 *         (Will be adjusted if needed) */
+		   int* iPtr,	/* OUT: Maximum number of digits
+				 *      to return */
+		   int *iLimPtr,/* OUT: Number of digits of significance
+				 *      if the bignum method is used.*/
+		   int *iLim1Ptr)
+				/* OUT: Number of digits of significance
+				 *      if the quick method is used. */
+{
+    switch(convType) {
+    case TCL_DD_SHORTEST0:
+    case TCL_DD_STEELE0:
+	*iLimPtr = *iLim1Ptr = -1;
+	*iPtr = 18;
+	*ndigitsPtr = 0;
+	break;
+    case TCL_DD_E_FORMAT:
+	if (*ndigitsPtr <= 0) {
+	    *ndigitsPtr = 1;
+	}
+	*iLimPtr = *iLim1Ptr = *iPtr = *ndigitsPtr;
+	break;
+    case TCL_DD_F_FORMAT:
+	*iPtr = *ndigitsPtr + k + 1;
+	*iLimPtr = *iPtr;
+	*iLim1Ptr = *iPtr - 1;
+	if (*iPtr <= 0) {
+	    *iPtr = 1;
+	}
+	break;
+    }
+}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * BumpUp --
+ *
+ *	Increases a string of digits ending in a series of nines to
+ *	designate the next higher number.  xxxxb9999... -> xxxx(b+1)0000...
+ *
+ * Results:
+ *	Returns a pointer to the end of the adjusted string.
+ *
+ * Side effects:
+ *	In the case that the string consists solely of '999999', sets it
+ *	to "1" and moves the decimal point (*kPtr) one place to the right.
+ *
+ *-----------------------------------------------------------------------------
+ */
+
+
+inline static char*
+BumpUp(char* s,		    	/* Cursor pointing one past the end of the
+				 * string */ 
+       char* retval,		/* Start of the string of digits */
+       int* kPtr)		/* Position of the decimal point */
+{
+    while (*--s == '9') {
+	if (s == retval) {
+	    ++(*kPtr);
+	    *s = '1';
+	    return s+1;
+	}
+    }
+    ++*s;
+    ++s;
+    return s;
+}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * AdjustRange --
+ *
+ *	Rescales a 'double' in preparation for formatting it using the
+ *	'quick' double-to-string method.
+ *
+ * Results:
+ *	Returns the precision that has been lost in the prescaling as
+ *	a count of units in the least significant place.
+ *
+ *-----------------------------------------------------------------------------
+ */
+
+inline static int
+AdjustRange(double* dPtr,	/* INOUT: Number to adjust */
+	    int k)		/* IN: floor(log10(d)) */
+{
+    int ieps;			/* Number of roundoff errors that have
+				 * accumulated */
+    double d = *dPtr;		/* Number to adjust */
+    double ds;
+    int i, j, j1;
+
+    ieps = 2;
+
+    if (k > 0) {
+	/*
+	 * The number must be reduced to bring it into range.
+	 */
+	ds = tens[k & 0xf];
+	j = k >> 4;
+	if (j & BLETCH) {
+	    j &= (BLETCH-1);
+	    d /= bigtens[N_BIGTENS - 1];
+	    ieps++;
+	}
+	i = 0;
+	for (; j != 0; j>>=1) {
+	    if (j & 1) {
+		ds *= bigtens[i];
+		++ieps;
+	    }
+	    ++i;
+	}
+	d /= ds;
+    } else if ((j1 = -k) != 0) {
+	/*
+	 * The number must be increased to bring it into range
+	 */
+	d *= tens[j1 & 0xf];
+	i = 0;
+	for (j = j1>>4; j; j>>=1) {
+	    if (j & 1) {
+		ieps++;
+		d *= bigtens[i];
+	    }
+	    ++i;
+	}
+    }
+
+    *dPtr = d;
+    return ieps;
+}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * ShorteningQuickFormat --
+ *
+ *	Returns a 'quick' format of a double precision number to a string
+ *	of digits, preferring a shorter string of digits if the shorter
+ *	string is still within 1/2 ulp of the number. 
+ *
+ * Results:
+ *	Returns the string of digits. Returns NULL if the 'quick' method
+ *	fails and the bignum method must be used.
+ *
+ * Side effects:
+ *	Stores the position of the decimal point at '*kPtr'.
+ *
+ *-----------------------------------------------------------------------------
+ */
+
+inline static char*
+ShorteningQuickFormat(double d,	/* Number to convert */
+		      int k,	/* floor(log10(d)) */
+		      int ilim,	/* Number of significant digits to return */
+		      double eps,
+				/* Estimated roundoff error */
+		      char* retval,
+				/* Buffer to receive the digit string */
+		      int* kPtr)
+				/* Pointer to stash the position of
+				 * the decimal point */
+{
+    char* s = retval;		/* Cursor in the return value */
+    int digit;			/* Current digit */
+    int i;
+
+    eps = 0.5 / tens[ilim-1] - eps;
+    i = 0;
+    for (;;) {
+	/* Convert a digit */
+
+	digit = d;
+	d -= digit;
+	*s++ = '0' + digit;
+
+	/* 
+	 * Truncate the conversion if the string of digits is within
+	 * 1/2 ulp of the actual value.
+	 */
+
+	if (d < eps) {
+	    *kPtr = k;
+	    return s;
+	}
+	if ((1. - d) < eps) {
+	    *kPtr = k;
+	    return BumpUp(s, retval, kPtr);
+	}
+
+	/*
+	 * Bail out if the conversion fails to converge to a sufficiently
+	 * precise value
+	 */
+
+	if (++i >= ilim) {
+	    return NULL;
+	}
+
+	/*
+	 * Bring the next digit to the integer part.
+	 */
+
+	eps *= 10;
+	d *= 10.0;
+    }
+}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * StrictQuickFormat --
+ *
+ *	Convert a double precision number of a string of a precise number
+ *	of digits, using the 'quick' double precision method.
+ *
+ * Results:
+ *	Returns the digit string, or NULL if the bignum method must be
+ *	used to do the formatting.
+ *
+ * Side effects:
+ *	Stores the position of the decimal point in '*kPtr'.
+ *
+ *-----------------------------------------------------------------------------
+ */
+
+inline static char*
+StrictQuickFormat(double d,	/* Number to convert */
+		  int k,	/* floor(log10(d)) */
+		  int ilim,	/* Number of significant digits to return */
+		  double eps,	/* Estimated roundoff error */
+		  char* retval,	/* Start of the digit string */
+		  int* kPtr)	/* Pointer to stash the position of
+				 * the decimal point */
+{
+    char* s = retval;		/* Cursor in the return value */
+    int digit;			/* Current digit of the answer */
+    int i;
+
+    eps *= tens[ilim-1];
+    i = 1;
+    for (;;) {
+	/* Extract a digit */
+	digit = d;
+	d -= digit;
+	if (d == 0.0) {
+	    ilim = i;
+	}
+	*s++ = '0' + digit;
+
+	/* 
+	 * When the given digit count is reached, handle trailing strings
+	 * of 0 and 9.
+	 */
+	if (i == ilim) {
+	    if (d > 0.5 + eps) {
+		*kPtr = k;
+		return BumpUp(s, retval, kPtr);
+	    } else if (d < 0.5 - eps) {
+		while (*--s == '0') {
+		    /* do nothing */
+		}
+		s++;
+		*kPtr = k;
+		return s;
+	    } else {
+		return NULL;
+	    }
+	}
+
+	/* Advance to the next digit */
+	++i;
+	d *= 10.0;
+    }
+}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * QuickConversion --
+ *
+ *	Converts a floating point number the 'quick' way, when only a limited
+ *	number of digits is required and floating point arithmetic can
+ *	therefore be used for the intermediate results.
+ *
+ * Results:
+ *	Returns the converted string, or NULL if the bignum method must
+ *	be used.
+ *
+ *-----------------------------------------------------------------------------
+ */
+
+inline static char*
+QuickConversion(double d,	/* Number to format */
+		int k,		/* floor(log10(d)), approximately */
+		int k_check,	/* 0 if k is exact, 1 if it may be too high */
+		int flags,	/* Flags passed to dtoa:
+				 *    TCL_DD_SHORTEN_FLAG */
+		int len,	/* Length of the return value */
+		int ilim,	/* Number of digits to store */
+		int ilim1,	/* Number of digits to store if we
+				 * musguessed k */
+		int* decpt,	/* OUTPUT: Location of the decimal point */
+		char** endPtr)	/* OUTPUT: Pointer to the terminal null byte */
+{
+    int ieps;			/* Number of 1-ulp roundoff errors that have
+				 * accumulated in the calculation*/
+    Double eps;			/* Estimated roundoff error */
+    char* retval;		/* Returned string */
+    char* end;			/* Pointer to the terminal null byte in the
+				 * returned string */
+
     /*
-     * 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.
+     * Bring d into the range [1 .. 10)
      */
+    ieps = AdjustRange(&d, k);
 
-    smallestSig = GetIntegerTimesPower(v, &r, &e);
+    /*
+     * If the guessed value of k didn't get d into range, adjust it
+     * by one. If that leaves us outside the range in which quick format
+     * is accurate, bail out.
+     */
+    if (k_check && d < 1. && ilim > 0) {
+	if (ilim1 < 0) {
+	    return NULL;
+	}
+	ilim = ilim1;
+	--k;
+	d *= 10.0;
+	++ieps;
+    }
 
-    lowOK = highOK = (mp_iseven(&r));
+    /*
+     * Compute estimated roundoff error
+     */
+    eps.d = ieps * d + 7.;
+    eps.w.word0 -= (FP_PRECISION-1) << EXP_SHIFT;
 
     /*
-     * We are going to want to develop integers r, s, mplus, and mminus such
-     * that v = r / s, v(+)-v / 2 = mplus / s; v-v(-) / 2 = mminus / s and
-     * then scale either s or r, mplus, mminus by an appropriate power of ten.
-     *
-     * We actually do this by keeping track of the powers of 2 and 5 by which
-     * f is multiplied to yield v and by which 1 is multiplied to yield s,
-     * mplus, and mminus.
+     * Handle the peculiar case where the result has no significant
+     * digits.
      */
+    retval = ckalloc(len + 1);
+    if (ilim == 0) {
+	d -= 5.;
+	if (d > eps.d) {
+	    *retval = '1';
+	    *decpt = k;
+	    return retval;
+	} else if (d < -eps.d) {
+	    *decpt = k;
+	    return retval;
+	} else {
+	    ckfree(retval);
+	    return NULL;
+	}
+    }
 
-    if (e >= 0) {
-	int bits = e * log2FLT_RADIX;
+    /* Format the digit string */
 
-	if (!smallestSig) {
-	    /*
-	     * Normal case, m+ and m- are both FLT_RADIX**e
-	     */
+    if (flags & TCL_DD_SHORTEN_FLAG) {
+	end = ShorteningQuickFormat(d, k, ilim, eps.d, retval, decpt);
+    } else {
+	end = StrictQuickFormat(d, k, ilim, eps.d, retval, decpt);
+    }
+    if (end == NULL) {
+	ckfree(retval);
+	return NULL;
+    }
+    *end = '\0';
+    if (endPtr != NULL) {
+	*endPtr = end;
+    }
+    return retval;
+}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * CastOutPowersOf2 --
+ *
+ *	Adjust the factors 'b2', 'm2', and 's2' to cast out common powers
+ *	of 2 from numerator and denominator in preparation for the 'bignum'
+ *	method of floating point conversion.
+ *
+ *-----------------------------------------------------------------------------
+ */
 
-	    rfac2 = bits + 1;
-	    sfac2 = 1;
-	    mplusfac2 = bits;
-	    mminusfac2 = bits;
+inline static void
+CastOutPowersOf2(int* b2,	/* Power of 2 to multiply the significand  */
+		 int* m2,	/* Power of 2 to multiply 1/2 ulp  */
+		 int* s2)	/* Power of 2 to multiply the common
+				 * denominator */
+{
+    int i;
+    if (*m2 > 0 && *s2 > 0) {	/* Find the smallest power of 2 in the
+				 * numerator */
+        if (*m2 < *s2) {	/* Find the lowest common denominatorr */
+	    i = *m2;
 	} else {
+	    i = *s2;
+	}
+	*b2 -= i;		/* Reduce to lowest terms */
+	*m2 -= i;
+	*s2 -= i;
+    }
+}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * ShorteningInt64Conversion --
+ *
+ *	Converts a double-precision number to the shortest string of
+ *	digits that reconverts exactly to the given number, or to
+ *	'ilim' digits if that will yield a shorter result. The numerator and
+ *	denominator in David Gay's conversion algorithm are known to fit
+ *	in Tcl_WideUInt, giving considerably faster arithmetic than mp_int's.
+ *
+ * Results:
+ *	Returns the string of significant decimal digits, in newly
+ *	allocated memory
+ *
+ * Side effects:
+ *	Stores the location of the decimal point in '*decpt' and the
+ *      location of the terminal null byte in '*endPtr'.
+ *
+ *-----------------------------------------------------------------------------
+ */
+
+inline static char*
+ShorteningInt64Conversion(Double* dPtr,
+				/* Original number to convert */
+			  int convType,
+				/* Type of conversion (shortest, Steele,
+				   E format, F format) */
+			  Tcl_WideUInt bw,
+				/* Integer significand */
+			  int b2, int b5,
+				/* Scale factor for the significand
+				 * in the numerator */
+			  int m2plus, int m2minus, int m5,
+			  	/* Scale factors for 1/2 ulp in
+				 * the numerator (will be different if
+				 * bw == 1 */
+			  int s2, int s5,
+				/* Scale factors for the denominator */
+			  int k,
+				/* Number of output digits before the decimal
+				 * point */
+			  int len,
+				/* Number of digits to allocate */
+			  int ilim,
+				/* Number of digits to convert if b >= s */
+			  int ilim1,
+				/* Number of digits to convert if b < s */
+			  int* decpt,
+				/* OUTPUT: Position of the decimal point */
+			  char** endPtr)
+				/* OUTPUT: Position of the terminal '\0'
+				 *         at the end of the returned string */
+{
+    
+    char* retval = ckalloc(len + 1);
+				/* Output buffer */
+    Tcl_WideUInt b = (bw * wuipow5[b5]) << b2;
+				/* Numerator of the fraction being converted */
+    Tcl_WideUInt S = wuipow5[s5] << s2;
+				/* Denominator of the fraction being 
+				 * converted */
+    Tcl_WideUInt mplus, mminus;	/* Ranges for testing whether the result
+				 * is within roundoff of being exact */
+    int digit;			/* Current output digit */
+    char* s = retval;		/* Cursor in the output buffer */
+    int i;			/* Current position in the output buffer */
+
+    /* Adjust if the logarithm was guessed wrong */
+
+    if (b < S) {
+	b = 10 * b;
+	++m2plus; ++m2minus; ++m5;
+	ilim = ilim1;
+	--k;
+    }
+
+    /* Compute roundoff ranges */
+
+    mplus = wuipow5[m5] << m2plus;
+    mminus = wuipow5[m5] << m2minus;
+
+    /* Loop through the digits */
+
+    i = 1;
+    for (;;) {
+	digit = (int)(b / S);
+	if (digit > 10) {
+	    Tcl_Panic("wrong digit!");
+	}
+	b = b % S;
+
+	/* 
+	 * Does the current digit put us on the low side of the exact value
+	 * but within within roundoff of being exact?
+	 */
+	if (b < mplus
+	    || (b == mplus
+		&& convType != TCL_DD_STEELE0
+		&& (dPtr->w.word1 & 1) == 0)) {
 	    /*
-	     * If f is equal to the smallest significand, then we need another
-	     * factor of FLT_RADIX in s to cope with stepping to the next
-	     * smaller exponent when going to e's predecessor.
+	     * Make sure we shouldn't be rounding *up* instead,
+	     * in case the next number above is closer
 	     */
+	    if (2 * b > S
+		|| (2 * b == S
+		    && (digit & 1) != 0)) {
+		++digit;
+		if (digit == 10) {
+		    *s++ = '9';
+		    s = BumpUp(s, retval, &k);
+		    break;
+		}
+	    }
 
-	    rfac2 = bits + log2FLT_RADIX + 1;
-	    sfac2 = 1 + log2FLT_RADIX;
-	    mplusfac2 = bits + log2FLT_RADIX;
-	    mminusfac2 = bits;
+	    /* Stash the current digit */
+
+	    *s++ = '0' + digit;
+	    break;
 	}
-    } else {
+
 	/*
-	 * v has digits after the binary point
+	 * Does one plus the current digit put us within roundoff of the
+	 * number?
 	 */
+	if (b > S - mminus
+	    || (b == S - mminus
+		&& convType != TCL_DD_STEELE0
+		&& (dPtr->w.word1 & 1) == 0)) {
+	    if (digit == 9) {
+		*s++ = '9';
+		s = BumpUp(s, retval, &k);
+		break;
+	    }
+	    ++digit;
+	    *s++ = '0' + digit;
+	    break;
+	}
 
-	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.
-	     */
+	/*
+	 * Have we converted all the requested digits?
+	 */
+	*s++ = '0' + digit;
+	if (i == ilim) {
+	    if (2*b > S
+		|| (2*b == S && (digit & 1) != 0)) {
+		s = BumpUp(s, retval, &k);
+	    }
+	    break;
+	}
+	
+	/* Advance to the next digit */
+	
+	b = 10 * b;
+	mplus = 10 * mplus;
+	mminus = 10 * mminus;
+	++i;
+    }
+
+    /* 
+     * Endgame - store the location of the decimal point and the end of the
+     * string.
+     */
+    *s = '\0';
+    *decpt = k;
+    if (endPtr) {
+	*endPtr = s;
+    }
+    return retval;
+}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * StrictInt64Conversion --
+ *
+ *	Converts a double-precision number to a fixed-length string of
+ *	'ilim' digits that reconverts exactly to the given number.
+ *	('ilim' should be replaced with 'ilim1' in the case where
+ *	log10(d) has been overestimated).  The numerator and
+ *	denominator in David Gay's conversion algorithm are known to fit
+ *	in Tcl_WideUInt, giving considerably faster arithmetic than mp_int's.
+ *
+ * Results:
+ *	Returns the string of significant decimal digits, in newly
+ *	allocated memory
+ *
+ * Side effects:
+ *	Stores the location of the decimal point in '*decpt' and the
+ *      location of the terminal null byte in '*endPtr'.
+ *
+ *-----------------------------------------------------------------------------
+ */
+
+inline static char*
+StrictInt64Conversion(Double* dPtr,
+				/* Original number to convert */
+		      int convType,
+				/* Type of conversion (shortest, Steele,
+				   E format, F format) */
+		      Tcl_WideUInt bw,
+				/* Integer significand */
+		      int b2, int b5,
+				/* Scale factor for the significand
+				 * in the numerator */
+		      int s2, int s5,
+				/* Scale factors for the denominator */
+		      int k,
+				/* Number of output digits before the decimal
+				 * point */
+		      int len,
+				/* Number of digits to allocate */
+		      int ilim,
+				/* Number of digits to convert if b >= s */
+		      int ilim1,
+				/* Number of digits to convert if b < s */
+		      int* decpt,
+				/* OUTPUT: Position of the decimal point */
+		      char** endPtr)
+				/* OUTPUT: Position of the terminal '\0'
+				 *         at the end of the returned string */
+{
+    
+    char* retval = ckalloc(len + 1);
+				/* Output buffer */
+    Tcl_WideUInt b = (bw * wuipow5[b5]) << b2;
+				/* Numerator of the fraction being converted */
+    Tcl_WideUInt S = wuipow5[s5] << s2;
+				/* Denominator of the fraction being 
+				 * converted */
+    int digit;			/* Current output digit */
+    char* s = retval;		/* Cursor in the output buffer */
+    int i;			/* Current position in the output buffer */
+
+    /* Adjust if the logarithm was guessed wrong */
+
+    if (b < S) {
+	b = 10 * b;
+	ilim = ilim1;
+	--k;
+    }
+
+    /* Loop through the digits */
+
+    i = 1;
+    for (;;) {
+	digit = (int)(b / S);
+	if (digit > 10) {
+	    Tcl_Panic("wrong digit!");
+	}
+	b = b % S;
+
+	/*
+	 * Have we converted all the requested digits?
+	 */
+	*s++ = '0' + digit;
+	if (i == ilim) {
+	    if (2*b > S
+		|| (2*b == S && (digit & 1) != 0)) {
+		s = BumpUp(s, retval, &k);
+	    }
+	    break;
+	}
+	
+	/* Advance to the next digit */
+	
+	b = 10 * b;
+	++i;
+    }
+
+    /* 
+     * Endgame - store the location of the decimal point and the end of the
+     * string.
+     */
+    *s = '\0';
+    *decpt = k;
+    if (endPtr) {
+	*endPtr = s;
+    }
+    return retval;
+}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * ShouldBankerRoundUpPowD --
+ *
+ *	Test whether bankers' rounding should round a digit up. Assumption
+ *	is made that the denominator of the fraction being tested is
+ *	a power of 2**DIGIT_BIT.
+ *
+ * Results:
+ *	Returns 1 iff the fraction is more than 1/2, or if the fraction
+ *	is exactly 1/2 and the digit is odd.
+ *
+ *-----------------------------------------------------------------------------
+ */
+
+inline static int
+ShouldBankerRoundUpPowD(mp_int* b,
+				/* Numerator of the fraction */
+			int sd,	/* Denominator is 2**(sd*DIGIT_BIT) */
+			int isodd)
+				/* 1 if the digit is odd, 0 if even */
+{
+    int i;
+    const static mp_digit topbit = (1<<(DIGIT_BIT-1));
+    if (b->used < sd || (b->dp[sd-1] & topbit) == 0) {
+	return 0;
+    }
+    if (b->dp[sd-1] != topbit) {
+	return 1;
+    }
+    for (i = sd-2; i >= 0; --i) {
+	if (b->dp[i] != 0) {
+	    return 1;
+	}
+    }
+    return isodd;
+}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * ShouldBankerRoundUpToNextPowD --
+ *
+ *	Tests whether bankers' rounding will round down in the
+ *	"denominator is a power of 2**MP_DIGIT" case.
+ *
+ * Results:
+ *	Returns 1 if the rounding will be performed - which increases the
+ *	digit by one - and 0 otherwise.
+ *
+ *-----------------------------------------------------------------------------
+ */
+
+inline static int
+ShouldBankerRoundUpToNextPowD(mp_int* b,
+				/* Numerator of the fraction */
+			      mp_int* m,
+				/* Numerator of the rounding tolerance */
+			      int sd, 
+				/* Common denominator is 2**(sd*DIGIT_BIT) */
+			      int convType,
+				/* Conversion type: STEELE defeats 
+				 * round-to-even (Not sure why one wants to
+				 * do this; I copied it from Gay) FIXME */
+			      int isodd,
+				/* 1 if the integer significand is odd */
+			      mp_int* temp)
+				/* Work area for the calculation */
+{
+    int i;
+
+    /* 
+     * Compare B and S-m -- which is the same as comparing B+m and S --
+     * which we do by computing b+m and doing a bitwhack compare against
+     * 2**(DIGIT_BIT*sd)
+     */
+    mp_add(b, m, temp);
+    if (temp->used <= sd) {	/* too few digits to be > S */
+	return 0;
+    }
+    if (temp->used > sd+1 || temp->dp[sd] > 1) {
+				/* >= 2s */
+	return 1;
+    }
+    for (i = sd-1; i >= 0; --i) {
+				/* check for ==s */
+	if (temp->dp[i] != 0) {	/* > s */
+	    return 1;
+	}
+    }
+    if (convType == TCL_DD_STEELE0) {
+				/* biased rounding */
+	return 0;
+    }
+    return isodd;
+}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * ShorteningBignumConversionPowD --
+ *
+ *	Converts a double-precision number to the shortest string of
+ *	digits that reconverts exactly to the given number, or to
+ *	'ilim' digits if that will yield a shorter result. The denominator
+ *	in David Gay's conversion algorithm is known to be a power of
+ *	2**DIGIT_BIT, and hence the division in the main loop may be replaced
+ *	by a digit shift and mask.
+ *
+ * Results:
+ *	Returns the string of significant decimal digits, in newly
+ *	allocated memory
+ *
+ * Side effects:
+ *	Stores the location of the decimal point in '*decpt' and the
+ *      location of the terminal null byte in '*endPtr'.
+ *
+ *-----------------------------------------------------------------------------
+ */
+
+inline static char*
+ShorteningBignumConversionPowD(Double* dPtr,
+				/* Original number to convert */
+			       int convType,
+				/* Type of conversion (shortest, Steele,
+				   E format, F format) */
+			       Tcl_WideUInt bw,
+				/* Integer significand */
+			       int b2, int b5,
+				/* Scale factor for the significand
+				 * in the numerator */
+			       int m2plus, int m2minus, int m5,
+			  	/* Scale factors for 1/2 ulp in
+				 * the numerator (will be different if
+				 * bw == 1 */
+			       int sd,
+				/* Scale factor for the denominator */
+			       int k,
+				/* Number of output digits before the decimal
+				 * point */
+			       int len,
+				/* Number of digits to allocate */
+			       int ilim,
+				/* Number of digits to convert if b >= s */
+			       int ilim1,
+				/* Number of digits to convert if b < s */
+			       int* decpt,
+				/* OUTPUT: Position of the decimal point */
+			       char** endPtr)
+				/* OUTPUT: Position of the terminal '\0'
+				 *         at the end of the returned string */
+{
+    
+    char* retval = ckalloc(len + 1);
+				/* Output buffer */
+    mp_int b;			/* Numerator of the fraction being converted */
+    mp_int mplus, mminus;	/* Bounds for roundoff */
+    mp_digit digit;		/* Current output digit */
+    char* s = retval;		/* Cursor in the output buffer */
+    int i;			/* Index in the output buffer */
+    mp_int temp;
+    int r1;
+
+    /* 
+     * b = bw * 2**b2 * 5**b5
+     * mminus = 5**m5
+     */
+
+    TclBNInitBignumFromWideUInt(&b, bw);
+    mp_init_set_int(&mminus, 1);
+    MulPow5(&b, b5, &b);
+    mp_mul_2d(&b, b2, &b);
+
+    /* Adjust if the logarithm was guessed wrong */
+
+    if (b.used <= sd) {
+	mp_mul_d(&b, 10, &b);
+	++m2plus; ++m2minus; ++m5;
+	ilim = ilim1;
+	--k;
+    }
+
+    /*
+     * mminus = 5**m5 * 2**m2minus
+     * mplus = 5**m5 * 2**m2plus
+     */
+
+    mp_mul_2d(&mminus, m2minus, &mminus);
+    MulPow5(&mminus, m5, &mminus);
+    if (m2plus > m2minus) {
+	mp_init_copy(&mplus, &mminus);
+	mp_mul_2d(&mplus, m2plus-m2minus, &mplus);
+    }
+    mp_init(&temp);
+
+    /* Loop through the digits. Do division and mod by s == 2**(sd*DIGIT_BIT)
+     * by mp_digit extraction */
 
-	    rfac2 = 1;
-	    sfac2 = 1 - e * log2FLT_RADIX;
-	    mplusfac2 = 0;
-	    mminusfac2 = 0;
+    i = 0;
+    for (;;) {
+	if (b.used <= sd) {
+	    digit = 0;
 	} else {
+	    digit = b.dp[sd];
+	    if (b.used > sd+1 || digit >= 10) {
+		Tcl_Panic("wrong digit!");
+	    }
+	    --b.used; mp_clamp(&b);
+	}
+
+	/* 
+	 * Does the current digit put us on the low side of the exact value
+	 * but within within roundoff of being exact?
+	 */
+	
+	r1 = mp_cmp_mag(&b, (m2plus > m2minus)? &mplus : &mminus);
+	if (r1 == MP_LT
+	    || (r1 == MP_EQ
+		&& convType != TCL_DD_STEELE0
+		&& (dPtr->w.word1 & 1) == 0)) {
 	    /*
-	     * f is the smallest significand, but e is not the smallest
-	     * exponent. We need to scale by FLT_RADIX again to cope with the
-	     * fact that v's predecessor has a smaller exponent.
+	     * Make sure we shouldn't be rounding *up* instead,
+	     * in case the next number above is closer
 	     */
+	    if (ShouldBankerRoundUpPowD(&b, sd, digit&1)) {
+		++digit;
+		if (digit == 10) {
+		    *s++ = '9';
+		    s = BumpUp(s, retval, &k);
+		    break;
+		}
+	    }
+
+	    /* Stash the last digit */
+
+	    *s++ = '0' + digit;
+	    break;
+	}
+
+	/*
+	 * Does one plus the current digit put us within roundoff of the
+	 * number?
+	 */
+	
+	if (ShouldBankerRoundUpToNextPowD(&b, &mminus, sd, 
+					   convType, dPtr->w.word1 & 1,
+					   &temp)) {
+	    if (digit == 9) {
+		*s++ = '9';
+		s = BumpUp(s, retval, &k);
+		break;
+	    }
+	    ++digit;
+	    *s++ = '0' + digit;
+	    break;
+	}
 
-	    rfac2 = 1 + log2FLT_RADIX;
-	    sfac2 = 1 + log2FLT_RADIX * (1 - e);
-	    mplusfac2 = FLT_RADIX;
-	    mminusfac2 = 0;
+	/*
+	 * Have we converted all the requested digits?
+	 */
+	*s++ = '0' + digit;
+	if (i == ilim) {
+	    if (ShouldBankerRoundUpPowD(&b, sd, digit&1)) {
+		s = BumpUp(s, retval, &k);
+	    }
+	    break;
+	}
+	
+	/* Advance to the next digit */
+	
+	mp_mul_d(&b, 10, &b);
+	mp_mul_d(&mminus, 10, &mminus);
+	if (m2plus > m2minus) {
+	    mp_mul_2d(&mminus, m2plus-m2minus, &mplus);
 	}
+	++i;
     }
 
-    /*
-     * Estimate the highest power of ten that will be needed to hold the
-     * result.
+    /* 
+     * Endgame - store the location of the decimal point and the end of the
+     * string.
+     */
+    if (m2plus > m2minus) {
+	mp_clear(&mplus);
+    }
+    mp_clear_multi(&b, &mminus, &temp, NULL);
+    *s = '\0';
+    *decpt = k;
+    if (endPtr) {
+	*endPtr = s;
+    }
+    return retval;
+}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * StrictBignumConversionPowD --
+ *
+ *	Converts a double-precision number to a fixed-lengt string of
+ *	'ilim' digits (or 'ilim1' if log10(d) has been overestimated.)
+ *	The denominator in David Gay's conversion algorithm is known to
+ *	be a power of 2**DIGIT_BIT, and hence the division in the main 
+ *	loop may be replaced by a digit shift and mask.
+ *
+ * Results:
+ *	Returns the string of significant decimal digits, in newly
+ *	allocated memory.
+ *
+ * Side effects:
+ *	Stores the location of the decimal point in '*decpt' and the
+ *      location of the terminal null byte in '*endPtr'.
+ *
+ *-----------------------------------------------------------------------------
+ */
+
+inline static char*
+StrictBignumConversionPowD(Double* dPtr,
+				/* Original number to convert */
+			   int convType,
+				/* Type of conversion (shortest, Steele,
+				   E format, F format) */
+			   Tcl_WideUInt bw,
+				/* Integer significand */
+			   int b2, int b5,
+				/* Scale factor for the significand
+				 * in the numerator */
+			   int sd,
+				/* Scale factor for the denominator */
+			   int k,
+				/* Number of output digits before the decimal
+				 * point */
+			   int len,
+				/* Number of digits to allocate */
+			   int ilim,
+				/* Number of digits to convert if b >= s */
+			   int ilim1,
+				/* Number of digits to convert if b < s */
+			   int* decpt,
+				/* OUTPUT: Position of the decimal point */
+			   char** endPtr)
+				/* OUTPUT: Position of the terminal '\0'
+				 *         at the end of the returned string */
+{
+    
+    char* retval = ckalloc(len + 1);
+				/* Output buffer */
+    mp_int b;			/* Numerator of the fraction being converted */
+    mp_digit digit;		/* Current output digit */
+    char* s = retval;		/* Cursor in the output buffer */
+    int i;			/* Index in the output buffer */
+    mp_int temp;
+
+    /* 
+     * b = bw * 2**b2 * 5**b5
      */
 
-    k = (int) ceil(log(v) / log(10.));
-    if (k >= 0) {
-	sfac2 += k;
-	sfac5 = k;
-    } else {
-	rfac2 -= k;
-	mplusfac2 -= k;
-	mminusfac2 -= k;
-	rfac5 = -k;
+    TclBNInitBignumFromWideUInt(&b, bw);
+    MulPow5(&b, b5, &b);
+    mp_mul_2d(&b, b2, &b);
+
+    /* Adjust if the logarithm was guessed wrong */
+
+    if (b.used <= sd) {
+	mp_mul_d(&b, 10, &b);
+	ilim = ilim1;
+	--k;
     }
+    mp_init(&temp);
 
-    /*
-     * Scale r, s, mplus, mminus by the appropriate powers of 2 and 5.
+    /* 
+     * Loop through the digits. Do division and mod by s == 2**(sd*DIGIT_BIT)
+     * by mp_digit extraction 
      */
 
-    mp_init_set(&mplus, 1);
-    for (i=0 ; i<=8 ; ++i) {
-	if (rfac5 & (1 << i)) {
-	    mp_mul(&mplus, pow5+i, &mplus);
+    i = 1;
+    for (;;) {
+	if (b.used <= sd) {
+	    digit = 0;
+	} else {
+	    digit = b.dp[sd];
+	    if (b.used > sd+1 || digit >= 10) {
+		Tcl_Panic("wrong digit!");
+	    }
+	    --b.used; mp_clamp(&b);
 	}
-    }
-    mp_mul(&r, &mplus, &r);
-    mp_mul_2d(&r, rfac2, &r);
-    mp_init_copy(&mminus, &mplus);
-    mp_mul_2d(&mplus, mplusfac2, &mplus);
-    mp_mul_2d(&mminus, mminusfac2, &mminus);
-    mp_init_set(&s, 1);
-    for (i=0 ; i<=8 ; ++i) {
-	if (sfac5 & (1 << i)) {
-	    mp_mul(&s, pow5+i, &s);
+
+	/*
+	 * Have we converted all the requested digits?
+	 */
+	*s++ = '0' + digit;
+	if (i == ilim) {
+	    if (ShouldBankerRoundUpPowD(&b, sd, digit&1)) {
+		s = BumpUp(s, retval, &k);
+	    }
+	    break;
 	}
+	
+	/* Advance to the next digit */
+	
+	mp_mul_d(&b, 10, &b);
+	++i;
     }
-    mp_mul_2d(&s, sfac2, &s);
 
-    /*
-     * It is possible for k to be off by one because we used an inexact
-     * logarithm.
+    /* 
+     * Endgame - store the location of the decimal point and the end of the
+     * string.
      */
+    mp_clear_multi(&b, &temp, NULL);
+    *s = '\0';
+    *decpt = k;
+    if (endPtr) {
+	*endPtr = s;
+    }
+    return retval;
+}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * ShouldBankerRoundUp --
+ *
+ *	Tests whether a digit should be rounded up or down when finishing
+ *	bignum-based floating point conversion.
+ *
+ * Results:
+ *	Returns 1 if the number needs to be rounded up, 0 otherwise.
+ *
+ *-----------------------------------------------------------------------------
+ */
 
-    mp_init(&temp);
-    mp_add(&r, &mplus, &temp);
-    i = mp_cmp_mag(&temp, &s);
-    if (i>0 || (highOK && i==0)) {
-	mp_mul_d(&s, 10, &s);
-	k++;
-    } else {
-	mp_mul_d(&temp, 10, &temp);
-	i = mp_cmp_mag(&temp, &s);
-	if (i<0 || (highOK && i==0)) {
-	    mp_mul_d(&r, 10, &r);
-	    mp_mul_d(&mplus, 10, &mplus);
-	    mp_mul_d(&mminus, 10, &mminus);
-	    k--;
+inline static int
+ShouldBankerRoundUp(mp_int* twor,
+				/* 2x the remainder from thd division that
+				 * produced the last digit */
+		    mp_int* S,	/* Denominator */
+		    int isodd)	/* Flag == 1 if the last digit is odd */
+{
+    int r = mp_cmp_mag(twor, S);
+    switch (r) {
+    case MP_LT:
+	return 0;
+    case MP_EQ:
+	return isodd;
+    case MP_GT:
+	return 1;
+    }
+    Tcl_Panic("in ShouldBankerRoundUp, trichotomy fails!");
+    return 0;
+}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * ShouldBankerRoundUpToNext --
+ *
+ *	Tests whether the remainder is great enough to force rounding
+ *	to the next higher digit.
+ *
+ * Results:
+ *	Returns 1 if the number should be rounded up, 0 otherwise.
+ *
+ *-----------------------------------------------------------------------------
+ */
+
+inline static int
+ShouldBankerRoundUpToNext(mp_int* b,
+				/* Remainder from the division that produced
+				 * the last digit. */
+			  mp_int* m,
+				/* Numerator of the rounding tolerance */
+			  mp_int* S,
+				/* Denominator */
+			  int convType,
+				/* Conversion type: STEELE0 defeats
+				 * round-to-even. (Not sure why one would
+				 * want this; I coped it from Gay. FIXME */
+			  int isodd,
+				/* 1 if the integer significand is odd */
+			  mp_int* temp)
+				/* Work area needed for the calculation */
+{
+    int r;
+    /* Compare b and S-m: this is the same as comparing B+m and S. */
+    mp_add(b, m, temp);
+    r = mp_cmp_mag(temp, S);
+    switch(r) {
+    case MP_LT:
+	return 0;
+    case MP_EQ:
+	if (convType == TCL_DD_STEELE0) {
+	    return 0;
+	} else {
+	    return isodd;
 	}
+    case MP_GT:
+	return 1;
     }
+    Tcl_Panic("in ShouldBankerRoundUpToNext, trichotomy fails!");
+    return 0;
+}
+		 
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * ShorteningBignumConversion --
+ *
+ *	Convert a floating point number to a variable-length digit string
+ *	using the multiprecision method.
+ *
+ * Results:
+ *	Returns the string of digits.
+ *
+ * Side effects:
+ *	Stores the position of the decimal point in *decpt.
+ *	Stores a pointer to the end of the number in *endPtr.
+ *
+ *-----------------------------------------------------------------------------
+ */
+
+inline static char*
+ShorteningBignumConversion(Double* dPtr,
+				/* Original number being converted */
+			   int convType,
+				/* Conversion type */
+			   Tcl_WideUInt bw,
+				/* Integer significand and exponent */
+			   int b2,
+				/* Scale factor for the significand */
+			   int m2plus, int m2minus,
+				/* Scale factors for 1/2 ulp in numerator */
+			   int s2, int s5,
+				/* Scale factors for denominator */
+			   int k,
+				/* Guessed position of the decimal point */
+			   int len,
+				/* Size of the digit buffer to allocate */
+			   int ilim,
+				/* Number of digits to convert if b >= s */
+			   int ilim1,
+				/* Number of digits to convert if b < s */
+			   int* decpt,
+				/* OUTPUT: Position of the decimal point */
+			   char** endPtr)
+				/* OUTPUT: Pointer to the end of the number */
+{
+    char* retval = ckalloc(len+1);
+				/* Buffer of digits to return */
+    char* s = retval;		/* Cursor in the return value */
+    mp_int b;			/* Numerator of the result */
+    mp_int mminus;		/* 1/2 ulp below the result */
+    mp_int mplus;		/* 1/2 ulp above the result */
+    mp_int S;			/* Denominator of the result */
+    mp_int dig;			/* Current digit of the result */
+    int digit;			/* Current digit of the result */
+    mp_int temp;		/* Work area */
+    int minit = 1;		/* Fudge factor for when we misguess k */
+    int i;
+    int r1;
 
     /*
-     * At this point, k contains the power of ten by which we're scaling the
-     * result. r/s is at least 1/10 and strictly less than ten, and v = r/s *
-     * 10**k. mplus and mminus give the rounding limits.
+     * b = bw * 2**b2 * 5**b5
+     * S = 2**s2 * 5*s5
      */
 
-    for (;;) {
-	int tc1, tc2;
+    TclBNInitBignumFromWideUInt(&b, bw);
+    mp_mul_2d(&b, b2, &b);
+    mp_init_set_int(&S, 1);
+    MulPow5(&S, s5, &S); mp_mul_2d(&S, s2, &S);
 
-	mp_mul_d(&r, 10, &r);
-	mp_div(&r, &s, &temp, &r);	/* temp = 10r / s; r = 10r mod s */
-	i = temp.dp[0];
-	mp_mul_d(&mplus, 10, &mplus);
-	mp_mul_d(&mminus, 10, &mminus);
-	tc1 = mp_cmp_mag(&r, &mminus);
-	if (lowOK) {
-	    tc1 = (tc1 <= 0);
-	} else {
-	    tc1 = (tc1 < 0);
+    /*
+     * Handle the case where we guess the position of the decimal point 
+     * wrong. 
+     */
+    
+    if (mp_cmp_mag(&b, &S) == MP_LT) {
+	mp_mul_d(&b, 10, &b);
+	minit = 10;
+	ilim =ilim1;
+	--k;
+    }
+
+    /* mminus = 2**m2minus * 5**m5 */
+
+    mp_init_set_int(&mminus, minit);
+    mp_mul_2d(&mminus, m2minus, &mminus);
+    if (m2plus > m2minus) {
+	mp_init_copy(&mplus, &mminus);
+	mp_mul_2d(&mplus, m2plus-m2minus, &mplus);
+    }
+    mp_init(&temp);
+
+    /* Loop through the digits */
+
+    mp_init(&dig);
+    i = 1;
+    for (;;) {
+	mp_div(&b, &S, &dig, &b);
+	if (dig.used > 1 || dig.dp[0] >= 10) {
+	    Tcl_Panic("wrong digit!");
 	}
-	mp_add(&r, &mplus, &temp);
-	tc2 = mp_cmp_mag(&temp, &s);
-	if (highOK) {
-	    tc2 = (tc2 >= 0);
-	} else {
-	    tc2 = (tc2 > 0);
+	digit = dig.dp[0];
+
+	/* 
+	 * Does the current digit leave us with a remainder small enough to
+	 * round to it?
+	 */
+
+	r1 = mp_cmp_mag(&b, (m2plus > m2minus)? &mplus : &mminus);
+	if (r1 == MP_LT
+	    || (r1 == MP_EQ
+		&& convType != TCL_DD_STEELE0
+		&& (dPtr->w.word1 & 1) == 0)) {
+		mp_mul_2d(&b, 1, &b);
+	    if (ShouldBankerRoundUp(&b, &S, digit&1)) {
+		++digit;
+		if (digit == 10) {
+		    *s++ = '9';
+		    s = BumpUp(s, retval, &k);
+		    break;
+		}
+	    }
+	    *s++ = '0' + digit;
+	    break;
 	}
-	if (!tc1) {
-	    if (!tc2) {
-		*buffer++ = '0' + i;
-	    } else {
-		c = (char) (i + '1');
+
+	/*
+	 * Does the current digit leave us with a remainder large enough
+	 * to commit to rounding up to the next higher digit?
+	 */
+
+	if (ShouldBankerRoundUpToNext(&b, &mminus, &S, convType,
+				      dPtr->w.word1 & 1, &temp)) {
+	    ++digit;
+	    if (digit == 10) {
+		*s++ = '9';
+		s = BumpUp(s, retval, &k);
 		break;
 	    }
-	} else {
-	    if (!tc2) {
-		c = (char) (i + '0');
-	    } else {
-		mp_mul_2d(&r, 1, &r);
-		n = mp_cmp_mag(&r, &s);
-		if (n < 0) {
-		    c = (char) (i + '0');
-		} else {
-		    c = (char) (i + '1');
-		}
+	    *s++ = '0' + digit;
+	    break;
+	}
+
+	/* Have we converted all the requested digits? */
+
+	*s++ = '0' + digit;
+	if (i == ilim) {
+	    mp_mul_2d(&b, 1, &b);
+	    if (ShouldBankerRoundUp(&b, &S, digit&1)) {
+		s = BumpUp(s, retval, &k);  
 	    }
 	    break;
 	}
-    };
-    *buffer++ = c;
-    *buffer++ = '\0';
 
-    /*
-     * Free memory, and return.
+	/* Advance to the next digit */
+
+	if (s5 > 0) {
+
+	    /* Can possibly shorten the denominator */
+	    mp_mul_2d(&b, 1, &b);
+	    mp_mul_2d(&mminus, 1, &mminus);
+	    if (m2plus > m2minus) {
+		mp_mul_2d(&mplus, 1, &mplus);
+	    }
+	    mp_div_d(&S, 5, &S, NULL);
+	    --s5;
+	    /* 
+	     * TODO: It might possibly be a win to fall back to
+	     *       int64 arithmetic here if S < 2**64/10. But it's
+	     *       a win only for a fairly narrow range of magnitudes
+	     *       so perhaps not worth bothering. We already know that
+	     *       we shorten the denominator by at least 1 mp_digit, perhaps
+	     *       2. as we do the conversion for 17 digits of significance.
+	     * Possible savings:
+	     * 10**26   1 trip through loop before fallback possible
+	     * 10**27   1 trip
+	     * 10**28   2 trips     
+	     * 10**29   3 trips
+	     * 10**30   4 trips
+	     * 10**31   5 trips
+	     * 10**32   6 trips
+	     * 10**33   7 trips
+	     * 10**34   8 trips
+	     * 10**35   9 trips
+	     * 10**36  10 trips
+	     * 10**37  11 trips
+	     * 10**38  12 trips
+	     * 10**39  13 trips
+	     * 10**40  14 trips
+	     * 10**41  15 trips
+	     * 10**42  16 trips
+	     * thereafter  no gain.
+	     */
+	} else {
+	    mp_mul_d(&b, 10, &b);
+	    mp_mul_d(&mminus, 10, &mminus);
+	    if (m2plus > m2minus) {
+		mp_mul_2d(&mplus, 10, &mplus);
+	    }
+	}
+
+	++i;
+    }
+
+
+    /* 
+     * Endgame - store the location of the decimal point and the end of the
+     * string.
      */
+    if (m2plus > m2minus) {
+	mp_clear(&mplus);
+    }
+    mp_clear_multi(&b, &mminus, &temp, NULL);
+    *s = '\0';
+    *decpt = k;
+    if (endPtr) {
+	*endPtr = s;
+    }
+    return retval;
 
-    mp_clear_multi(&r, &s, &mplus, &mminus, &temp, NULL);
-    return k;
 }
-
+		 
 /*
- *----------------------------------------------------------------------
+ *-----------------------------------------------------------------------------
  *
- * AbsoluteValue --
+ * StrictBignumConversion --
  *
- *	Splits a 'double' into its absolute value and sign.
+ *	Convert a floating point number to a fixed-length digit string
+ *	using the multiprecision method.
  *
  * Results:
- *	Returns the absolute value.
+ *	Returns the string of digits.
  *
  * Side effects:
- *	Stores the signum in '*signum'.
+ *	Stores the position of the decimal point in *decpt.
+ *	Stores a pointer to the end of the number in *endPtr.
  *
- *----------------------------------------------------------------------
+ *-----------------------------------------------------------------------------
  */
 
-static double
-AbsoluteValue(
-    double v,			/* Number to split */
-    int *signum)		/* (Output) Sign of the number 1=-, 0=+ */
+inline static char*
+StrictBignumConversion(Double* dPtr,
+				/* Original number being converted */
+		       int convType,
+				/* Conversion type */
+		       Tcl_WideUInt bw,
+				/* Integer significand and exponent */
+		       int b2,	/* Scale factor for the significand */
+		       int s2, int s5,
+				/* Scale factors for denominator */
+		       int k,	/* Guessed position of the decimal point */
+		       int len,	/* Size of the digit buffer to allocate */
+		       int ilim,
+				/* Number of digits to convert if b >= s */
+		       int ilim1,
+				/* Number of digits to convert if b < s */
+		       int* decpt,
+				/* OUTPUT: Position of the decimal point */
+		       char** endPtr)
+				/* OUTPUT: Pointer to the end of the number */
 {
+    char* retval = ckalloc(len+1);
+				/* Buffer of digits to return */
+    char* s = retval;		/* Cursor in the return value */
+    mp_int b;			/* Numerator of the result */
+    mp_int S;			/* Denominator of the result */
+    mp_int dig;			/* Current digit of the result */
+    int digit;			/* Current digit of the result */
+    mp_int temp;		/* Work area */
+    int g;			/* Size of the current digit groun */
+    int i, j;
+    
     /*
-     * 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.)
+     * b = bw * 2**b2 * 5**b5
+     * S = 2**s2 * 5*s5
      */
 
-#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 (n770_fp) {
-	bitwhack.iv = Nokia770Twiddle(bitwhack.iv);
+    TclBNInitBignumFromWideUInt(&b, bw);
+    mp_mul_2d(&b, b2, &b);
+    mp_init_set_int(&S, 1);
+    MulPow5(&S, s5, &S); mp_mul_2d(&S, s2, &S);
+
+    /*
+     * Handle the case where we guess the position of the decimal point 
+     * wrong. 
+     */
+    
+    if (mp_cmp_mag(&b, &S) == MP_LT) {
+	mp_mul_d(&b, 10, &b);
+	ilim =ilim1;
+	--k;
     }
-    if (bitwhack.iv & ((Tcl_WideUInt) 1 << 63)) {
-	*signum = 1;
-	bitwhack.iv &= ~((Tcl_WideUInt) 1 << 63);
-	if (n770_fp) {
-	    bitwhack.iv = Nokia770Twiddle(bitwhack.iv);
+    mp_init(&temp);
+
+    /* Convert the leading digit */
+
+    mp_init(&dig);
+    i = 0;
+    mp_div(&b, &S, &dig, &b);
+    if (dig.used > 1 || dig.dp[0] >= 10) {
+	    Tcl_Panic("wrong digit!");
+	}
+    digit = dig.dp[0];
+
+    /* Is a single digit all that was requested? */
+
+    *s++ = '0' + digit;
+    if (++i >= ilim) {
+	mp_mul_2d(&b, 1, &b);
+	if (ShouldBankerRoundUp(&b, &S, digit&1)) {
+	    s = BumpUp(s, retval, &k);  
 	}
-	v = bitwhack.dv;
     } else {
-	*signum = 0;
+
+	for (;;) {
+
+	    /* Shift by a group of digits. */
+
+	    g = ilim - i;
+	    if (g > DIGIT_GROUP) {
+		g = DIGIT_GROUP;
+	    }
+	    if (s5 >= g) {
+		mp_div_d(&S, dpow5[g], &S, NULL);
+		s5 -= g;
+	    } else if (s5 > 0) {
+		mp_div_d(&S, dpow5[s5], &S, NULL);
+		mp_mul_d(&b, dpow5[g - s5], &b);
+		s5 = 0;
+	    } else {
+		mp_mul_d(&b, dpow5[g], &b);
+	    }
+	    mp_mul_2d(&b, g, &b);
+	    
+	    /*
+	     * As with the shortening bignum conversion, it's possible at
+	     * this point that we will have reduced the denominator to
+	     * less than 2**64/10, at which point it would be possible to
+	     * fall back to to int64 arithmetic. But the potential payoff
+	     * is tremendously less - unless we're working in F format -
+	     * because we know that three groups of digits will always
+	     * suffice for %#.17e, the longest format that doesn't introduce
+	     * empty precision.
+	     */
+
+	    /* Extract the next digit */
+	    
+	    mp_div(&b, &S, &dig, &b);
+	    if (dig.used > 1) {
+		Tcl_Panic("wrong digit!");
+	    }
+	    digit = dig.dp[0];
+	    for (j = g-1; j >= 0; --j) {
+		int t = itens[j];
+		*s++ = digit / t + '0';
+		digit %= t;
+	    }
+	    i += g;
+	    
+	    /* Have we converted all the requested digits? */
+	    
+	    if (i == ilim) {
+		mp_mul_2d(&b, 1, &b);
+		if (ShouldBankerRoundUp(&b, &S, digit&1)) {
+		    s = BumpUp(s, retval, &k);  
+		}
+		break;
+	    }
+	}
     }
-#endif
-    return v;
+    /* 
+     * Endgame - store the location of the decimal point and the end of the
+     * string.
+     */
+    mp_clear_multi(&b, &temp, NULL);
+    *s = '\0';
+    *decpt = k;
+    if (endPtr) {
+	*endPtr = s;
+    }
+    return retval;
+
 }
 
 /*
- *----------------------------------------------------------------------
+ *-----------------------------------------------------------------------------
  *
- * GetIntegerTimesPower --
+ * TclDoubleDigits --
  *
- *	Converts a floating point number to an exact integer times a power of
- *	the floating point radix.
+ *	Core of Tcl's conversion of double-precision floating point numbers
+ *	to decimal.
  *
  * Results:
- *	Returns 1 if it converted the smallest significand, 0 otherwise.
+ *	Returns a newly-allocated string of digits.
  *
  * Side effects:
- *	Initializes the integer value (does not just assign it), and stores
- *	the exponent.
+ *	Sets *decpt to the index of the character in the string before the
+ *	place that the decimal point should go. If 'endPtr' is not NULL,
+ *	sets endPtr to point to the terminating '\0' byte of the string.
+ *	Sets *sign to 1 if a minus sign should be printed with the number,
+ *	or 0 if a plus sign (or no sign) should appear.
  *
- *----------------------------------------------------------------------
+ * This function is a service routine that produces the string of digits
+ * for floating-point-to-decimal conversion. It can do a number of things
+ * according to the 'flags' argument. Valid values for 'flags' include:
+ *	TCL_DD_SHORTEST - This is the default for floating point conversion
+ *		if ::tcl_precision is 0. It constructs the shortest string
+ *		of digits that will reconvert to the given number when scanned.
+ *		For floating point numbers that are exactly between two
+ *		decimal numbers, it resolves using the 'round to even' rule.
+ *		With this value, the 'ndigits' parameter is ignored.
+ *	TCL_DD_STEELE - This value is not recommended and may be removed
+ *		in the future. It follows the conversion algorithm outlined
+ *		in "How to Print Floating-Point Numbers Accurately" by
+ *		Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, 
+ *		pp. 112-126]. This rule has the effect of rendering 1e23
+ *		as 9.9999999999999999e22 - which is a 'better' approximation
+ *		in the sense that it will reconvert correctly even if
+ *		a subsequent input conversion is 'round up' or 'round down'
+ *		rather than 'round to nearest', but is surprising otherwise.
+ *	TCL_DD_E_FORMAT - This value is used to prepare numbers for %e
+ *		format conversion (or for default floating->string if
+ *		tcl_precision is not 0). It constructs a string of at most
+ *		'ndigits' digits, choosing the one that is closest to the
+ *		given number (and resolving ties with 'round to even').
+ *		It is allowed to return fewer than 'ndigits' if the number
+ *		converts exactly; if the TCL_DD_E_FORMAT|TCL_DD_SHORTEN_FLAG 
+ *		is supplied instead, it is also allowed to return fewer digits
+ *		if the shorter string will still reconvert to the given
+ *		input number.
+ *	TCL_DD_F_FORMAT - This value is used to prepare numbers for %f
+ *		format conversion. It requests that conversion proceed until
+ *		'ndigits' digits after the decimal point have been converted.
+ *		It is possible for this format to result in a zero-length 
+ *		string if the number is sufficiently small. Again, it
+ *		is permissible for TCL_DD_F_FORMAT to return fewer digits
+ *		for a number that converts exactly, and changing the
+ *		argument to TCL_DD_F_FORMAT|TCL_DD_SHORTEN_FLAG will allow
+ *		the routine also to return fewer digits if the shorter string
+ *		will still reconvert without loss to the given input number.
+ *
+ *	To any of these flags may be OR'ed TCL_DD_NO_QUICK; this flag
+ *	requires all calculations to be done in exact arithmetic. Normally,
+ *	E and F format with fewer than about 14 digits will be done with
+ *	a quick floating point approximation and fall back on the exact
+ *	arithmetic only if the input number is close enough to the
+ *	midpoint between two decimal strings that more precision is needed
+ *	to resolve which string is correct.
+ *
+ * The value stored in the 'decpt' argument on return may be negative 
+ * (indicating that the decimal point falls to the left of the string) 
+ * or greater than the length of the string.  In addition, the value -9999
+ * is used as a sentinel to indicate that the string is one of the special
+ * values "Infinity" and "NaN", and that no decimal point should be inserted.
+ * 
+ *-----------------------------------------------------------------------------
  */
-
-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*/
+char*
+TclDoubleDigits(double dv,	/* Number to convert */
+		int ndigits,	/* Number of digits requested */
+		int flags,	/* Conversion flags */
+		int* decpt,	/* OUTPUT: Position of the decimal point */
+		int* sign,	/* OUTPUT: 1 if the result is negative */
+		char** endPtr)	/* OUTPUT: If not NULL, receives a pointer
+				 *         to one character beyond the end
+				 *         of the returned string */
 {
-    double a, f;
-    int e, i, n;
+    int convType = (flags & TCL_DD_CONVERSION_TYPE_MASK);
+				/* Type of conversion being performed
+				 * TCL_DD_SHORTEST0
+				 * TCL_DD_STEELE0
+				 * TCL_DD_E_FORMAT
+				 * TCL_DD_F_FORMAT */
+    Double d;			/* Union for deconstructing doubles */
+    Tcl_WideUInt bw;		/* Integer significand */
+    int be;			/* Power of 2 by which b must be multiplied */
+    int bbits;			/* Number of bits needed to represent b */
+    int denorm;			/* Flag == 1 iff the input number was
+				 * denormalized */
+    int k;			/* Estimate of floor(log10(d)) */
+    int k_check;		/* Flag == 1 if d is near enough to a 
+				 * power of ten that k must be checked */
+    int b2, b5, s2, s5;		/* Powers of 2 and 5 in the numerator and
+				 * denominator of intermediate results */
+    int ilim, ilim1;
+    char* retval;		/* Return value from this function */
+    int i;
 
-    /*
-     * Develop f and e such that v = f * FLT_RADIX**e, with
-     * 1.0/FLT_RADIX <= f < 1.
+    /* Put the input number into a union for bit-whacking */
+
+    d.d = dv;
+
+    /* 
+     * Handle the cases of negative numbers (by taking the absolute value:
+     * this includes -Inf and -NaN!), infinity, Not a Number, and zero.
      */
 
-    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);
+    TakeAbsoluteValue(&d, sign);
+    if ((d.w.word0 & EXP_MASK) == EXP_MASK) {
+	return FormatInfAndNaN(&d, decpt, endPtr);
     }
-    e = (e - n) / log2FLT_RADIX;
-#endif
-    if (f == 1.0) {
-	f = 1.0 / FLT_RADIX;
-	e += 1;
+    if (d.d == 0.0) {
+	return FormatZero(decpt, endPtr);
     }
 
+    /* 
+     * Unpack the floating point into a wide integer and an exponent.
+     * Determine the number of bits that the big integer requires, and
+     * compute a quick approximation (which may be one too high) of
+     * ceil(log10(d.d)).
+     */
+    denorm = ((d.w.word0 & EXP_MASK) == 0);
+    DoubleToExpAndSig(d.d, &bw, &be, &bbits);
+    k = ApproximateLog10(bw, be, bbits);
+    k = BetterLog10(d.d, k, &k_check);
+
+    /* At this point, we have:
+     *	  d is the number to convert.
+     *    bw are significand and exponent: d == bw*2**be, 
+     *    bbits is the length of bw: 2**bbits-1 <= bw < 2**bbits
+     *	  k is either ceil(log10(d)) or ceil(log10(d))+1. k_check is 0
+     *      if we know that k is exactly ceil(log10(d)) and 1 if we need to
+     *	    check.
+     *    We want a rational number 
+     *      r = b * 10**(1-k) = bw * 2**b2 * 5**b5 / (2**s2 / 5**s5),
+     *    with b2, b5, s2, s5 >= 0.  Note that the most significant decimal
+     *    digit is floor(r) and that successive digits can be obtained
+     *    by setting r <- 10*floor(r) (or b <= 10 * (b % S)).
+     *    Find appropriate b2, b5, s2, s5.
+     */
+
+    ComputeScale(be, k, &b2, &b5, &s2, &s5);
+
     /*
-     * If the original number was denormalized, adjust e and f to be denormal
-     * as well.
+     * Correct an incorrect caller-supplied 'ndigits'.
+     * Also determine:
+     *	i = The maximum number of decimal digits that will be returned in the
+     *      formatted string.  This is k + 1 + ndigits for F format, 18 for
+     *	    shortest and Steele, and ndigits for E format.
+     *  ilim = The number of significant digits to convert if
+     *         k has been guessed correctly. This is -1 for shortest and Steele
+     *         (which stop when all significance has been lost), 'ndigits'
+     *	       for E format, and 'k + 1 + ndigits' for F format.
+     *  ilim1 = The minimum number of significant digits to convert if
+     *	        k has been guessed 1 too high. This, too, is -1 for shortest
+     *	        and Steele, and 'ndigits' for E format, but it's 'ndigits-1'
+     *	        for F format.
      */
 
-    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;
+    SetPrecisionLimits(convType, k, &ndigits, &i, &ilim, &ilim1);
+
+    /* 
+     * Try to do low-precision conversion in floating point rather
+     * than resorting to expensive multiprecision arithmetic
+     */
+    if (ilim >= 0 && ilim <= QUICK_MAX && !(flags & TCL_DD_NO_QUICK)) {
+	if ((retval = QuickConversion(d.d, k, k_check, flags,
+				      i, ilim, ilim1,
+				      decpt, endPtr)) != NULL) {
+	    return retval;
+	}
     }
 
-    /*
-     * 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);
+    /* 
+     * For shortening conversions, determine the upper and lower bounds
+     * for the remainder at which we can stop.
+     *   m+ = (2**m2plus * 5**m5) / (2**s2 * 5**s5) is the limit on the
+     *        high side, and
+     *   m- = (2**m2minus * 5**m5) / (2**s2 * 5**s5) is the limit on the
+     *        low side.
+     *   We may need to increase s2 to put m2plus, m2minus, b2 over a
+     *   common denominator.
+     */
+
+    if (flags & TCL_DD_SHORTEN_FLAG) {
+	int m2minus = b2;
+	int m2plus;
+	int m5 = b5;
+	int len = i;
+
+	/* 
+	 * Find the quantity i so that (2**i*5**b5)/(2**s2*5**s5)
+	 * is 1/2 unit in the least significant place of the floating 
+	 * point number.
+	 */
+	if (denorm) {
+	    i = be + EXPONENT_BIAS + (FP_PRECISION-1);
+	} else {
+	    i = 1 + FP_PRECISION - bbits;
+	}
+	b2 += i;
+	s2 += i;
+
+	/* 
+	 * Reduce the fractions to lowest terms, since the above calculation
+	 * may have left excess powers of 2 in numerator and denominator
+	 */
+	CastOutPowersOf2(&b2, &m2minus, &s2);
+
+	/*
+	 * In the special case where bw==1, the nearest floating point number
+	 * to it on the low side is 1/4 ulp below it. Adjust accordingly.
+	 */
+	m2plus = m2minus;
+	if (!denorm && bw == 1) {
+	    ++b2;
+	    ++s2;
+	    ++m2plus;
+	}
+
+	if (s5+1 < N_LOG2POW5
+	    && s2+1 + log2pow5[s5+1] <= 64) {
+	    /*
+	     * If 10*2**s2*5**s5 == 2**(s2+1)+5**(s5+1) fits in a 64-bit
+	     * word, then all our intermediate calculations can be done
+	     * using exact 64-bit arithmetic with no need for expensive
+	     * multiprecision operations. (This will be true for all numbers
+	     * in the range [1.0e-3 .. 1.0e+24]).
+	     */
+
+	    return ShorteningInt64Conversion(&d, convType, bw, b2, b5,
+					     m2plus, m2minus, m5,
+					     s2, s5, k, len, ilim, ilim1, 
+					     decpt, endPtr);
+	} else if (s5 == 0) {
+	    /*
+	     * The denominator is a power of 2, so we can replace division
+	     * by digit shifts. First we round up s2 to a multiple of
+	     * DIGIT_BIT, and adjust m2 and b2 accordingly. Then we launch
+	     * into a version of the comparison that's specialized for
+	     * the 'power of mp_digit in the denominator' case.
+	     */
+	    if (s2 % DIGIT_BIT != 0) {
+		int delta = DIGIT_BIT - (s2 % DIGIT_BIT);
+		b2 += delta;
+		m2plus += delta;
+		m2minus += delta;
+		s2 += delta;
+	    }
+	    return ShorteningBignumConversionPowD(&d, convType, bw, b2, b5,
+						  m2plus, m2minus, m5,
+						  s2/DIGIT_BIT, k, len, 
+						  ilim, ilim1, decpt, endPtr);
+	} else {
+
+	    /* 
+	     * Alas, there's no helpful special case; use full-up
+	     * bignum arithmetic for the conversion
+	     */
+
+	    return ShorteningBignumConversion(&d, convType, bw,
+					      b2, m2plus, m2minus,
+					      s2, s5, k, len,
+					      ilim, ilim1, decpt, endPtr);
+
+	}
+
+    } else {
+
+	/* Non-shortening conversion */
+
+	int len = i;
+
+	/* Reduce numerator and denominator to lowest terms */
+
+	if (b2 >= s2 && s2 > 0) {
+	    b2 -= s2; s2 = 0;
+	} else if (s2 >= b2 && b2 > 0) {
+	    s2 -= b2; b2 = 0;
+	}
+
+	if (s5+1 < N_LOG2POW5
+	    && s2+1 + log2pow5[s5+1] <= 64) {
+	    /*
+	     * If 10*2**s2*5**s5 == 2**(s2+1)+5**(s5+1) fits in a 64-bit
+	     * word, then all our intermediate calculations can be done
+	     * using exact 64-bit arithmetic with no need for expensive
+	     * multiprecision operations.
+	     */
+
+	    return StrictInt64Conversion(&d, convType, bw, b2, b5,
+					 s2, s5, k, len, ilim, ilim1, 
+					 decpt, endPtr);
+
+	} else if (s5 == 0) {
+	    /*
+	     * The denominator is a power of 2, so we can replace division
+	     * by digit shifts. First we round up s2 to a multiple of
+	     * DIGIT_BIT, and adjust m2 and b2 accordingly. Then we launch
+	     * into a version of the comparison that's specialized for
+	     * the 'power of mp_digit in the denominator' case.
+	     */
+	    if (s2 % DIGIT_BIT != 0) {
+		int delta = DIGIT_BIT - (s2 % DIGIT_BIT);
+		b2 += delta;
+		s2 += delta;
+	    }
+	    return StrictBignumConversionPowD(&d, convType, bw, b2, b5,
+					      s2/DIGIT_BIT, k, len, 
+					      ilim, ilim1, decpt, endPtr);
+	} else {
+	    /*
+	     * There are no helpful special cases, but at least we know
+	     * in advance how many digits we will convert. We can run the
+	     * conversion in steps of DIGIT_GROUP digits, so as to
+	     * have many fewer mp_int divisions.
+	     */
+	    return StrictBignumConversion(&d, convType, bw, b2, s2, s5,
+					  k, len, ilim, ilim1, decpt, endPtr);
+	}
+    } 
 }
+
 
 /*
  *----------------------------------------------------------------------
@@ -2237,6 +4380,11 @@ TclInitDoubleConversion(void)
     for (i=0; i<8; ++i) {
 	mp_sqr(pow5+i, pow5+i+1);
     }
+    mp_init_set_int(pow5_13, 1220703125);
+    for (i = 1; i < 5; ++i) {
+	mp_init(pow5_13 + i);
+	mp_sqr(pow5_13 + i - 1, pow5_13 + i);
+    }
 
     /*
      * Determine the number of decimal digits to the left and right of the
@@ -2293,7 +4441,7 @@ TclFinalizeDoubleConversion(void)
 {
     int i;
 
-    Tcl_Free((char *) pow10_wide);
+    ckfree((char *) pow10_wide);
     for (i=0; i<9; ++i) {
 	mp_clear(pow5 + i);
     }
@@ -2433,6 +4581,20 @@ TclBignumToDouble(
 	return -r;
     }
 }
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * TclCeil --
+ *
+ *	Computes the smallest floating point number that is at least the
+ *	mp_int argument.
+ *
+ * Results:
+ *	Returns the floating point number.
+ *
+ *-----------------------------------------------------------------------------
+ */
 
 double
 TclCeil(
@@ -2476,6 +4638,20 @@ TclCeil(
     mp_clear(&b);
     return r;
 }
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * TclFloor --
+ *
+ *	Computes the largest floating point number less than or equal to
+ *	the mp_int argument.
+ *
+ * Results:
+ *	Returns the floating point value.
+ *
+ *-----------------------------------------------------------------------------
+ */
 
 double
 TclFloor(
diff --git a/generic/tclStubInit.c b/generic/tclStubInit.c
index b804b9a..d308f34 100644
--- a/generic/tclStubInit.c
+++ b/generic/tclStubInit.c
@@ -8,7 +8,7 @@
  * See the file "license.terms" for information on usage and redistribution
  * of this file, and for a DISCLAIMER OF ALL WARRANTIES.
  *
- * RCS: @(#) $Id: tclStubInit.c,v 1.197 2010/09/16 14:49:37 nijtmans Exp $
+ * RCS: @(#) $Id: tclStubInit.c,v 1.198 2010/11/28 23:20:11 kennykb Exp $
  */
 
 #include "tclInt.h"
@@ -305,6 +305,7 @@ static const TclIntStubs tclIntStubs = {
     TclInitRewriteEnsemble, /* 246 */
     TclResetRewriteEnsemble, /* 247 */
     TclCopyChannel, /* 248 */
+    TclDoubleDigits, /* 249 */
 };
 
 static const TclIntPlatStubs tclIntPlatStubs = {
@@ -460,6 +461,8 @@ const TclTomMathStubs tclTomMathStubs = {
     TclBN_s_mp_mul_digs, /* 58 */
     TclBN_s_mp_sqr, /* 59 */
     TclBN_s_mp_sub, /* 60 */
+    TclBN_mp_init_set_int, /* 61 */
+    TclBN_mp_set_int, /* 62 */
 };
 
 static const TclStubHooks tclStubHooks = {
diff --git a/generic/tclTest.c b/generic/tclTest.c
index b845d72..80dbbae 100644
--- a/generic/tclTest.c
+++ b/generic/tclTest.c
@@ -14,9 +14,11 @@
  * See the file "license.terms" for information on usage and redistribution of
  * this file, and for a DISCLAIMER OF ALL WARRANTIES.
  *
- * RCS: @(#) $Id: tclTest.c,v 1.154 2010/09/27 19:42:38 msofer Exp $
+ * RCS: @(#) $Id: tclTest.c,v 1.155 2010/11/28 23:20:11 kennykb Exp $
  */
 
+#include <math.h>
+
 #undef STATIC_BUILD
 #ifndef USE_TCL_STUBS
 #   define USE_TCL_STUBS
@@ -233,6 +235,9 @@ static int		TestdelCmd(ClientData dummy,
 			    Tcl_Interp *interp, int argc, const char **argv);
 static int		TestdelassocdataCmd(ClientData dummy,
 			    Tcl_Interp *interp, int argc, const char **argv);
+static int		TestdoubledigitsObjCmd(ClientData dummy,
+					       Tcl_Interp* interp,
+					       int objc, Tcl_Obj* const objv[]);
 static int		TestdstringCmd(ClientData dummy,
 			    Tcl_Interp *interp, int argc, const char **argv);
 static int		TestencodingObjCmd(ClientData dummy,
@@ -569,6 +574,8 @@ Tcltest_Init(
     Tcl_CreateCommand(interp, "testdel", TestdelCmd, NULL, NULL);
     Tcl_CreateCommand(interp, "testdelassocdata", TestdelassocdataCmd,
 	    NULL, NULL);
+    Tcl_CreateObjCommand(interp, "testdoubledigits", TestdoubledigitsObjCmd,
+			 NULL, NULL);
     Tcl_DStringInit(&dstring);
     Tcl_CreateCommand(interp, "testdstring", TestdstringCmd, NULL,
 	    NULL);
@@ -1594,6 +1601,102 @@ TestdelassocdataCmd(
 }
 
 /*
+ *-----------------------------------------------------------------------------
+ *
+ * TestdoubledigitsCmd --
+ *
+ *	This procedure implements the 'testdoubledigits' command. It is
+ *	used to test the low-level floating-point formatting primitives
+ *	in Tcl.
+ *
+ * Usage:
+ *	testdoubledigits fpval ndigits type ?shorten"
+ *
+ * Parameters:
+ *	fpval - Floating-point value to format.
+ *	ndigits - Digit count to request from Tcl_DoubleDigits
+ *	type - One of 'shortest', 'Steele', 'e', 'f'
+ *	shorten - Indicates that the 'shorten' flag should be passed in.
+ *
+ *-----------------------------------------------------------------------------
+ */
+
+static int
+TestdoubledigitsObjCmd(ClientData unused,
+				/* NULL */
+		       Tcl_Interp* interp,
+				/* Tcl interpreter */
+		       int objc,
+				/* Parameter count */
+		       Tcl_Obj* const objv[])
+				/* Parameter vector */
+{
+    static const char* options[] = {
+	"shortest",
+	"Steele",
+	"e",
+	"f",
+	NULL
+    };
+    static const int types[] = {
+	TCL_DD_SHORTEST,
+	TCL_DD_STEELE,
+	TCL_DD_E_FORMAT,
+	TCL_DD_F_FORMAT
+    };
+
+    const Tcl_ObjType* doubleType;
+    double d;
+    int status;
+    int ndigits;
+    int type;
+    int decpt;
+    int signum;
+    char* str;
+    char* endPtr;
+    Tcl_Obj* strObj;
+    Tcl_Obj* retval;
+
+    if (objc < 4 || objc > 5) {
+	Tcl_WrongNumArgs(interp, 1, objv, "fpval ndigits type ?shorten?");
+	return TCL_ERROR;
+    }
+    status = Tcl_GetDoubleFromObj(interp, objv[1], &d);
+    if (status != TCL_OK) {
+	doubleType = Tcl_GetObjType("double");
+	if (objv[1]->typePtr == doubleType
+	    || TclIsNaN(objv[1]->internalRep.doubleValue)) {
+	    status = TCL_OK;
+	    memcpy(&d, &(objv[1]->internalRep.doubleValue), sizeof(double));
+	}
+    }
+    if (status != TCL_OK
+	|| Tcl_GetIntFromObj(interp, objv[2], &ndigits) != TCL_OK
+	|| Tcl_GetIndexFromObj(interp, objv[3], options, "conversion type",
+			       TCL_EXACT, &type) != TCL_OK) {
+	fprintf(stderr, "bad value? %g\n", d);
+	return TCL_ERROR;
+    }
+    type = types[type];
+    if (objc > 4) {
+	if (strcmp(Tcl_GetString(objv[4]), "shorten")) {
+	    Tcl_SetObjResult(interp, Tcl_NewStringObj("bad flag", -1));
+	    return TCL_ERROR;
+	}
+	type |= TCL_DD_SHORTEN_FLAG;
+    }
+    str = TclDoubleDigits(d, ndigits, type, &decpt, &signum, &endPtr);
+    strObj = Tcl_NewStringObj(str, endPtr-str);
+    ckfree(str);
+    retval = Tcl_NewListObj(1, &strObj);
+    Tcl_ListObjAppendElement(NULL, retval, Tcl_NewIntObj(decpt));
+    strObj = Tcl_NewStringObj(signum ? "-" : "+", 1);
+    Tcl_ListObjAppendElement(NULL, retval, strObj);
+    Tcl_SetObjResult(interp, retval);
+    return TCL_OK;
+}
+
+/*
  *----------------------------------------------------------------------
  *
  * TestdstringCmd --
diff --git a/generic/tclTomMath.decls b/generic/tclTomMath.decls
index 12fa7cd..0bd5ca5 100644
--- a/generic/tclTomMath.decls
+++ b/generic/tclTomMath.decls
@@ -13,7 +13,7 @@
 # See the file "license.terms" for information on usage and redistribution
 # of this file, and for a DISCLAIMER OF ALL WARRANTIES.
 #
-# RCS: @(#) $Id: tclTomMath.decls,v 1.8 2010/09/15 07:33:54 nijtmans Exp $
+# RCS: @(#) $Id: tclTomMath.decls,v 1.9 2010/11/28 23:20:11 kennykb Exp $
 
 library tcl
 
@@ -214,3 +214,9 @@ declare 59 {
 declare 60 {
     int TclBN_s_mp_sub(mp_int *a, mp_int *b, mp_int *c)
 }
+declare 61 {
+    int TclBN_mp_init_set_int(mp_int* a, unsigned long i)
+}
+declare 62 {
+    int TclBN_mp_set_int(mp_int* a, unsigned long i)
+}
\ No newline at end of file
diff --git a/generic/tclTomMathDecls.h b/generic/tclTomMathDecls.h
index 61e9c07..4e55ce1 100644
--- a/generic/tclTomMathDecls.h
+++ b/generic/tclTomMathDecls.h
@@ -11,7 +11,7 @@
  * See the file "license.terms" for information on usage and redistribution
  * of this file, and for a DISCLAIMER OF ALL WARRANTIES.
  *
- * RCS: @(#) $Id: tclTomMathDecls.h,v 1.13 2010/08/19 04:26:04 nijtmans Exp $
+ * RCS: @(#) $Id: tclTomMathDecls.h,v 1.14 2010/11/28 23:20:11 kennykb Exp $
  */
 
 #ifndef _TCLTOMMATHDECLS
@@ -79,6 +79,7 @@
 #define mp_init_copy TclBN_mp_init_copy
 #define mp_init_multi TclBN_mp_init_multi
 #define mp_init_set TclBN_mp_init_set
+#define mp_init_set_int TclBN_mp_init_set_int
 #define mp_init_size TclBN_mp_init_size
 #define mp_karatsuba_mul TclBN_mp_karatsuba_mul
 #define mp_karatsuba_sqr TclBN_mp_karatsuba_sqr
@@ -96,6 +97,7 @@
 #define mp_rshd TclBN_mp_rshd
 #define mp_s_rmap TclBNMpSRmap
 #define mp_set TclBN_mp_set
+#define mp_set_int TclBN_mp_set_int
 #define mp_shrink TclBN_mp_shrink
 #define mp_sqr TclBN_mp_sqr
 #define mp_sqrt TclBN_mp_sqrt
@@ -268,6 +270,10 @@ EXTERN int		TclBN_s_mp_mul_digs(mp_int *a, mp_int *b, mp_int *c,
 EXTERN int		TclBN_s_mp_sqr(mp_int *a, mp_int *b);
 /* 60 */
 EXTERN int		TclBN_s_mp_sub(mp_int *a, mp_int *b, mp_int *c);
+/* 61 */
+EXTERN int		TclBN_mp_init_set_int(mp_int*a, unsigned long i);
+/* 62 */
+EXTERN int		TclBN_mp_set_int(mp_int*a, unsigned long i);
 
 typedef struct TclTomMathStubs {
     int magic;
@@ -334,6 +340,8 @@ typedef struct TclTomMathStubs {
     int (*tclBN_s_mp_mul_digs) (mp_int *a, mp_int *b, mp_int *c, int digs); /* 58 */
     int (*tclBN_s_mp_sqr) (mp_int *a, mp_int *b); /* 59 */
     int (*tclBN_s_mp_sub) (mp_int *a, mp_int *b, mp_int *c); /* 60 */
+    int (*tclBN_mp_init_set_int) (mp_int*a, unsigned long i); /* 61 */
+    int (*tclBN_mp_set_int) (mp_int*a, unsigned long i); /* 62 */
 } TclTomMathStubs;
 
 #ifdef __cplusplus
@@ -472,6 +480,10 @@ extern const TclTomMathStubs *tclTomMathStubsPtr;
 	(tclTomMathStubsPtr->tclBN_s_mp_sqr) /* 59 */
 #define TclBN_s_mp_sub \
 	(tclTomMathStubsPtr->tclBN_s_mp_sub) /* 60 */
+#define TclBN_mp_init_set_int \
+	(tclTomMathStubsPtr->tclBN_mp_init_set_int) /* 61 */
+#define TclBN_mp_set_int \
+	(tclTomMathStubsPtr->tclBN_mp_set_int) /* 62 */
 
 #endif /* defined(USE_TCL_STUBS) */
 
diff --git a/generic/tclUtil.c b/generic/tclUtil.c
index 45cc8a1..a04c29c 100644
--- a/generic/tclUtil.c
+++ b/generic/tclUtil.c
@@ -11,7 +11,7 @@
  * See the file "license.terms" for information on usage and redistribution of
  * this file, and for a DISCLAIMER OF ALL WARRANTIES.
  *
- * RCS: @(#) $Id: tclUtil.c,v 1.118 2010/10/01 12:52:50 dkf Exp $
+ * RCS: @(#) $Id: tclUtil.c,v 1.119 2010/11/28 23:20:11 kennykb Exp $
  */
 
 #include "tclInt.h"
@@ -2234,129 +2234,100 @@ Tcl_PrintDouble(
     char *p, c;
     int exponent;
     int signum;
-    char buffer[TCL_DOUBLE_SPACE];
+    char* digits;
+    char* end;
     Tcl_UniChar ch;
 
     int *precisionPtr = Tcl_GetThreadData(&precisionKey, (int)sizeof(int));
 
     /*
-     * If *precisionPtr == 0, then use TclDoubleDigits to develop a decimal
-     * significand and exponent, then format it in E or F format as
-     * appropriate. If *precisionPtr != 0, use the native sprintf and then add
-     * a trailing ".0" if there is no decimal point in the rep.
+     * Handle NaN.
      */
+    
+    if (TclIsNaN(value)) {
+	TclFormatNaN(value, dst);
+	return;
+    }
 
-    if (*precisionPtr == 0) {
+    /*
+     * Handle infinities.
+     */
+    
+    if (TclIsInfinite(value)) {
 	/*
-	 * Handle NaN.
+	 * Remember to copy the terminating NUL too.
 	 */
-
-	if (TclIsNaN(value)) {
-	    TclFormatNaN(value, dst);
-	    return;
+	
+	if (value < 0) {
+	    memcpy(dst, "-Inf", 5);
+	} else {
+	    memcpy(dst, "Inf", 4);
 	}
+	return;
+    }
 
+    /*
+     * Ordinary (normal and denormal) values.
+     */
+    
+    if (*precisionPtr == 0) {
+	digits = TclDoubleDigits(value, -1, TCL_DD_SHORTEST,
+				 &exponent, &signum, &end);
+    } else {
+	digits = TclDoubleDigits(value, *precisionPtr, TCL_DD_E_FORMAT, 
+				 &exponent, &signum, &end);
+    }
+    if (signum) {
+	*dst++ = '-';
+    }
+    p = digits;
+    if (exponent < -4 || exponent > 16) {
 	/*
-	 * Handle infinities.
+	 * E format for numbers < 1e-3 or >= 1e17.
 	 */
-
-	if (TclIsInfinite(value)) {
-	    /*
-	     * Remember to copy the terminating NUL too.
-	     */
-
-	    if (value < 0) {
-		memcpy(dst, "-Inf", 5);
-	    } else {
-		memcpy(dst, "Inf", 4);
+	
+	*dst++ = *p++;
+	c = *p;
+	if (c != '\0') {
+	    *dst++ = '.';
+	    while (c != '\0') {
+		*dst++ = c;
+		c = *++p;
 	    }
-	    return;
 	}
-
+	sprintf(dst, "e%+d", exponent);
+    } else {
 	/*
-	 * Ordinary (normal and denormal) values.
+	 * F format for others.
 	 */
-
-	exponent = TclDoubleDigits(buffer, value, &signum);
-	if (signum) {
-	    *dst++ = '-';
+	
+	if (exponent < 0) {
+	    *dst++ = '0';
 	}
-	p = buffer;
-	if (exponent < -3 || exponent > 17) {
-	    /*
-	     * E format for numbers < 1e-3 or >= 1e17.
-	     */
-
-	    *dst++ = *p++;
-	    c = *p;
+	c = *p;
+	while (exponent-- >= 0) {
 	    if (c != '\0') {
-		*dst++ = '.';
-		while (c != '\0') {
-		    *dst++ = c;
-		    c = *++p;
-		}
-	    }
-	    sprintf(dst, "e%+d", exponent-1);
-	} else {
-	    /*
-	     * F format for others.
-	     */
-
-	    if (exponent <= 0) {
-		*dst++ = '0';
-	    }
-	    c = *p;
-	    while (exponent-- > 0) {
-		if (c != '\0') {
-		    *dst++ = c;
-		    c = *++p;
-		} else {
-		    *dst++ = '0';
-		}
-	    }
-	    *dst++ = '.';
-	    if (c == '\0') {
-		*dst++ = '0';
+		*dst++ = c;
+		c = *++p;
 	    } else {
-		while (++exponent < 0) {
-		    *dst++ = '0';
-		}
-		while (c != '\0') {
-		    *dst++ = c;
-		    c = *++p;
-		}
+		*dst++ = '0';
 	    }
-	    *dst++ = '\0';
 	}
-    } else {
-	/*
-	 * tcl_precision is supplied, pass it to the native sprintf.
-	 */
-
-	sprintf(dst, "%.*g", *precisionPtr, value);
-
-	/*
-	 * If the ASCII result looks like an integer, add ".0" so that it
-	 * doesn't look like an integer anymore. This prevents floating-point
-	 * values from being converted to integers unintentionally. Check for
-	 * ASCII specifically to speed up the function.
-	 */
-
-	for (p = dst; *p != 0;) {
-	    if (UCHAR(*p) < 0x80) {
-		c = *p++;
-	    } else {
-		p += Tcl_UtfToUniChar(p, &ch);
-		c = UCHAR(ch);
+	*dst++ = '.';
+	if (c == '\0') {
+	    *dst++ = '0';
+	} else {
+	    while (++exponent < -1) {
+		*dst++ = '0';
 	    }
-	    if ((c == '.') || isalpha(UCHAR(c))) {	/* INTL: ISO only. */
-		return;
+	    while (c != '\0') {
+		*dst++ = c;
+		c = *++p;
 	    }
 	}
-	p[0] = '.';
-	p[1] = '0';
-	p[2] = 0;
+	*dst++ = '\0';
     }
+    ckfree(digits);
 }
 
 /*
diff --git a/tests/util.test b/tests/util.test
index 994fc0f..1d31609 100644
--- a/tests/util.test
+++ b/tests/util.test
@@ -7,13 +7,14 @@
 # See the file "license.terms" for information on usage and redistribution
 # of this file, and for a DISCLAIMER OF ALL WARRANTIES.
 #
-# RCS: @(#) $Id: util.test,v 1.20 2008/10/14 16:35:44 dgp Exp $
+# RCS: @(#) $Id: util.test,v 1.21 2010/11/28 23:20:11 kennykb Exp $
 
 if {[lsearch [namespace children] ::tcltest] == -1} {
     package require tcltest
     namespace import -force ::tcltest::*
 }
 
+testConstraint controversialNaN 1
 testConstraint testdstring [llength [info commands testdstring]]
 testConstraint testconcatobj [llength [info commands testconcatobj]]
 
@@ -43,6 +44,10 @@ proc testIEEE {} {
 		ieeeValues(+Infinity)
 	    binary scan \x00\x00\x00\x00\x00\x00\xf8\x7f d \
 		ieeeValues(NaN)
+	    binary scan \x00\x00\x00\x00\x00\x00\xf8\xff d \
+		ieeeValues(-NaN)
+	    binary scan \xef\xcd\xab\x89\x67\x45\xfb\xff d \
+		ieeeValues(-NaN(3456789abcdef))
 	    set ieeeValues(littleEndian) 1
 	    return 1
 	}
@@ -65,6 +70,10 @@ proc testIEEE {} {
 		ieeeValues(+Infinity)
 	    binary scan \x7f\xf8\x00\x00\x00\x00\x00\x00 d \
 		ieeeValues(NaN)
+	    binary scan \xff\xf8\x00\x00\x00\x00\x00\x00 d \
+		ieeeValues(-NaN)
+	    binary scan \xff\xfb\x45\x67\x89\xab\xcd\xef d \
+		ieeeValues(-NaN(3456789abcdef))
 	    set ieeeValues(littleEndian) 0
 	    return 1
 	}
@@ -85,6 +94,30 @@ proc convertDouble { x } {
     return $result
 }
 
+proc verdonk_test {sig binexp shouldbe exp} {
+    regexp {([-+]?)([0-9a-f]+)} $sig -> signum sig
+    scan $sig %llx sig
+    if {$signum eq {-}} {
+	set signum [expr 1<<63]
+    } else {
+	set signum 0
+    }
+    regexp {E([-+]?[0-9]+)} $binexp -> binexp
+    set word [expr {$signum | (($binexp + 0x3ff)<<52)|($sig & ~(1<<52))}]
+    binary scan [binary format w $word] q double
+    regexp {([-+])(\d+)_(\d+)\&} $shouldbe -> signum digits1 digits2
+    regexp {E([-+]\d+)} $exp -> decexp
+    incr decexp [expr {[string length $digits1] - 1}]
+    lassign [testdoubledigits $double [string length $digits1] e] \
+	outdigits decpt outsign
+    if {[string index $digits2 0] >= 5} {
+	incr digits1
+    }
+    if {$outsign != $signum || $outdigits != $digits1 || $decpt != $decexp} {
+	return -code error "result is ${outsign}0.${outdigits}E$decpt\
+                            should be ${signum}0.${digits1}E$decexp"
+    }
+}
 
 test util-1.1 {TclFindElement procedure - binary element in middle of list} {
     lindex {0 foo\x00help 1} 1
@@ -1106,6 +1139,774 @@ test util-11.23 {Tcl_PrintDouble - scaling} {
     expr 1.1e17
 } {1.1e+17}
 
+test util-12.1 {TclDoubleDigits - Inf} ieeeFloatingPoint {
+     testdoubledigits Inf -1 shortest
+} {Infinity 9999 +}
+test util-12.2 {TclDoubleDigits - -Inf} ieeeFloatingPoint {
+     testdoubledigits -Inf -1 shortest
+} {Infinity 9999 -}
+test util-12.3 {TclDoubleDigits - NaN} ieeeFloatingPoint {
+     testdoubledigits $ieeeValues(NaN) -1 shortest
+} {NaN 9999 +}
+test util-12.4 {TclDoubleDigits - NaN} {*}{
+     -constraints {ieeeFloatingPoint && controversialNaN}
+     -body {
+	 testdoubledigits -NaN -1 shortest
+     }
+    -result {NaN 9999 -}
+}
+test util-12.5 {TclDoubleDigits - 0} {
+     testdoubledigits 0.0 -1 shortest
+} {0 0 +}
+test util-12.6 {TclDoubleDigits - -0} {
+     testdoubledigits -0.0 -1 shortest
+} {0 0 -}
+
+# Verdonk test vectors
+
+test util-13.1 {just over exact - 1 digits} {*}{
+    -body {
+        verdonk_test 1754e31cd072da E+1008 +4_000000000000000000& E+303
+    }
+    -result {}
+}
+test util-13.2 {just over exact - 1 digits} {*}{
+    -body {
+        verdonk_test -1afcef51f0fb5f E+265 -1_000000000000000000& E+80
+    }
+    -result {}
+}
+test util-13.3 {just over exact - 1 digits} {*}{
+    -body {
+        verdonk_test 1754e31cd072da E+1006 +1_000000000000000000& E+303
+    }
+    -result {}
+}
+test util-13.4 {just over exact - 1 digits} {*}{
+    -body {
+        verdonk_test -1754e31cd072da E+1007 -2_000000000000000000& E+303
+    }
+    -result {}
+}
+test util-13.5 {just over exact - 1 digits} {*}{
+    -body {
+        verdonk_test 1e07b27dd78b14 E-848 +1_00000000000000000& E-255
+    }
+    -result {}
+}
+test util-13.6 {just over exact - 1 digits} {*}{
+    -body {
+        verdonk_test -1e29e9c56687fe E-709 -7_00000000000000000& E-214
+    }
+    -result {}
+}
+test util-13.7 {just over exact - 1 digits} {*}{
+    -body {
+        verdonk_test 1be03d0bf225c7 E-137 +1_00000000000000000& E-41
+    }
+    -result {}
+}
+test util-13.8 {just over exact - 1 digits} {*}{
+    -body {
+        verdonk_test -1a2fe76a3f9475 E-499 -1_00000000000000000& E-150
+    }
+    -result {}
+}
+test util-13.9 {just under exact - 1 digits} {*}{
+    -body {
+        verdonk_test 19a2028368022e E+1019 +8_999999999999999999& E+306
+    }
+    -result {}
+}
+test util-13.10 {just under exact - 1 digits} {*}{
+    -body {
+        verdonk_test -1317e5ef3ab327 E+509 -1_999999999999999999& E+153
+    }
+    -result {}
+}
+test util-13.11 {just under exact - 1 digits} {*}{
+    -body {
+        verdonk_test 1317e5ef3ab327 E+510 +3_99999999999999999& E+153
+    }
+    -result {}
+}
+test util-13.12 {just under exact - 1 digits} {*}{
+    -body {
+        verdonk_test -1317e5ef3ab327 E+511 -7_99999999999999999& E+153
+    }
+    -result {}
+}
+test util-13.13 {just under exact - 1 digits} {*}{
+    -body {
+        verdonk_test 1eb8e84fa0b278 E-1008 +6_999999999999999999& E-304
+    }
+    -result {}
+}
+test util-13.14 {just under exact - 1 digits} {*}{
+    -body {
+        verdonk_test -13339131c46f8b E-1004 -6_999999999999999999& E-303
+    }
+    -result {}
+}
+test util-13.15 {just under exact - 1 digits} {*}{
+    -body {
+        verdonk_test 1c0f92a6276c9d E-162 +2_999999999999999999& E-49
+    }
+    -result {}
+}
+test util-13.16 {just under exact - 1 digits} {*}{
+    -body {
+        verdonk_test -15ce1f143d7ad2 E-443 -5_99999999999999999& E-134
+    }
+    -result {}
+}
+test util-13.17 {just over exact - 2 digits} {*}{
+    -body {
+        verdonk_test 1c0794d9d40e96 E-301 +43_000000000000000000& E-92
+    }
+    -result {}
+}
+test util-13.18 {just over exact - 2 digits} {*}{
+    -body {
+        verdonk_test -1c0794d9d40e96 E-300 -86_000000000000000000& E-92
+    }
+    -result {}
+}
+test util-13.19 {just over exact - 2 digits} {*}{
+    -body {
+        verdonk_test 1cd5bee57763e6 E-241 +51_000000000000000000& E-74
+    }
+    -result {}
+}
+test util-13.20 {just under exact - 2 digits} {*}{
+    -body {
+        verdonk_test 1d1c26db7d0dae E+651 +16_999999999999999999& E+195
+    }
+    -result {}
+}
+test util-13.21 {just under exact - 2 digits} {*}{
+    -body {
+        verdonk_test -13f7ced916872b E-5 -38_999999999999999999& E-3
+    }
+    -result {}
+}
+test util-13.22 {just over exact - 3 digits} {*}{
+    -body {
+        verdonk_test 17d93193f78fc6 E+588 +151_0000000000000000000& E+175
+    }
+    -result {}
+}
+test util-13.23 {just over exact - 3 digits} {*}{
+    -body {
+        verdonk_test -1a82a1631eeb30 E-625 -119_000000000000000000& E-190
+    }
+    -result {}
+}
+test util-13.24 {just under exact - 3 digits} {*}{
+    -body {
+        verdonk_test -16c309024bab4b E+290 -282_999999999999999999& E+85
+    }
+    -result {}
+}
+test util-13.25 {just over exact - 8 digits} {*}{
+    -body {
+        verdonk_test 1dbbac6f83a821 E-800 +27869147_0000000000000000000& E-248
+    }
+    -result {}
+}
+test util-13.26 {just under exact - 9 digits} {*}{
+    -body {
+        verdonk_test -1c569e968e0944 E+430 -491080653_9999999999999999999& E+121
+    }
+    -result {}
+}
+test util-13.27 {just under exact - 9 digits} {*}{
+    -body {
+        verdonk_test 1c569e968e0944 E+429 +245540326_9999999999999999999& E+121
+    }
+    -result {}
+}
+test util-13.28 {just over exact - 10 digits} {*}{
+    -body {
+        verdonk_test -1fc575867314ee E-330 -9078555839_0000000000000000000& E-109
+    }
+    -result {}
+}
+test util-13.29 {just under exact - 10 digits} {*}{
+    -body {
+        verdonk_test -1c569e968e0944 E+428 -1227701634_9999999999999999999& E+120
+    }
+    -result {}
+}
+test util-13.30 {just over exact - 11 digits} {*}{
+    -body {
+        verdonk_test 1fc575867314ee E-329 +18157111678_0000000000000000000& E-109
+    }
+    -result {}
+}
+test util-13.31 {just over exact - 14 digits} {*}{
+    -body {
+        verdonk_test -18bf7e7fa6f02a E-196 -15400733123779_0000000000000000000& E-72
+    }
+    -result {}
+}
+test util-13.32 {just over exact - 17 digits} {*}{
+    -body {
+        verdonk_test -13de005bd620df E+217 -26153245263757307_0000000000000000000& E+49
+    }
+    -result {}
+}
+test util-13.33 {just over exact - 18 digits} {*}{
+    -body {
+        verdonk_test 1f92bacb3cb40c E+718 +272104041512242479_0000000000000000000& E+199
+    }
+    -result {}
+}
+test util-13.34 {just over exact - 18 digits} {*}{
+    -body {
+        verdonk_test -1f92bacb3cb40c E+719 -544208083024484958_0000000000000000000& E+199
+    }
+    -result {}
+}
+test util-13.35 {just over half ulp - 1 digits} {*}{
+    -body {
+        verdonk_test 142dbf25096cf5 E+148 +4_500000000000000000& E+44
+    }
+    -result {}
+}
+test util-13.36 {just over half ulp - 1 digits} {*}{
+    -body {
+        verdonk_test -1afcef51f0fb5f E+263 -2_500000000000000000& E+79
+    }
+    -result {}
+}
+test util-13.37 {just over half ulp - 1 digits} {*}{
+    -body {
+        verdonk_test 102498ea6df0c4 E+145 +4_500000000000000000& E+43
+    }
+    -result {}
+}
+test util-13.38 {just over half ulp - 1 digits} {*}{
+    -body {
+        verdonk_test -1754e31cd072da E+1004 -2_500000000000000000& E+302
+    }
+    -result {}
+}
+test util-13.39 {just over half ulp - 1 digits} {*}{
+    -body {
+        verdonk_test 12deac01e2b4f7 E-557 +2_50000000000000000& E-168
+    }
+    -result {}
+}
+test util-13.40 {just over half ulp - 1 digits} {*}{
+    -body {
+        verdonk_test -1b1df536c13eee E-307 -6_50000000000000000& E-93
+    }
+    -result {}
+}
+test util-13.41 {just over half ulp - 1 digits} {*}{
+    -body {
+        verdonk_test 10711fed5b19a4 E-154 +4_50000000000000000& E-47
+    }
+    -result {}
+}
+test util-13.42 {just over half ulp - 1 digits} {*}{
+    -body {
+        verdonk_test -148d67e8b1e00d E-151 -4_50000000000000000& E-46
+    }
+    -result {}
+}
+test util-13.43 {just under half ulp - 1 digits} {*}{
+    -body {
+        verdonk_test 1c8c574c0c6be7 E+187 +3_49999999999999999& E+56
+    }
+    -result {}
+}
+test util-13.44 {just under half ulp - 1 digits} {*}{
+    -body {
+        verdonk_test -1756183c147514 E+206 -1_49999999999999999& E+62
+    }
+    -result {}
+}
+test util-13.45 {just under half ulp - 1 digits} {*}{
+    -body {
+        verdonk_test 12ab469676c410 E+203 +1_49999999999999999& E+61
+    }
+    -result {}
+}
+test util-13.46 {just under half ulp - 1 digits} {*}{
+    -body {
+        verdonk_test -1539684e774b48 E+246 -1_49999999999999999& E+74
+    }
+    -result {}
+}
+test util-13.47 {just under half ulp - 1 digits} {*}{
+    -body {
+        verdonk_test 12e5f5dfa4fe9d E-286 +9_499999999999999999& E-87
+    }
+    -result {}
+}
+test util-13.48 {just under half ulp - 1 digits} {*}{
+    -body {
+        verdonk_test -1bdc2417bf7787 E-838 -9_499999999999999999& E-253
+    }
+    -result {}
+}
+test util-13.49 {just under half ulp - 1 digits} {*}{
+    -body {
+        verdonk_test 1eb8e84fa0b278 E-1009 +3_499999999999999999& E-304
+    }
+    -result {}
+}
+test util-13.50 {just under half ulp - 1 digits} {*}{
+    -body {
+        verdonk_test -1e3cbc9907fdc8 E-290 -9_499999999999999999& E-88
+    }
+    -result {}
+}
+test util-13.51 {just over half ulp - 2 digits} {*}{
+    -body {
+        verdonk_test 10ad836f269a17 E-324 +30_500000000000000000& E-99
+    }
+    -result {}
+}
+test util-13.52 {just over half ulp - 2 digits} {*}{
+    -body {
+        verdonk_test -1b39ae1909c31b E-687 -26_500000000000000000& E-208
+    }
+    -result {}
+}
+test util-13.53 {just over half ulp - 3 digits} {*}{
+    -body {
+        verdonk_test 1b2ab18615fcc6 E-576 +686_500000000000000000& E-176
+    }
+    -result {}
+}
+test util-13.54 {just over half ulp - 3 digits} {*}{
+    -body {
+        verdonk_test -13e1f90a573064 E-624 -178_500000000000000000& E-190
+    }
+    -result {}
+}
+test util-13.55 {just under half ulp - 3 digits} {*}{
+    -body {
+        verdonk_test 16c309024bab4b E+289 +141_499999999999999999& E+85
+    }
+    -result {}
+}
+test util-13.56 {just under half ulp - 4 digits} {*}{
+    -body {
+        verdonk_test -159bd3ad46e346 E+193 -1695_499999999999999999& E+55
+    }
+    -result {}
+}
+test util-13.57 {just under half ulp - 4 digits} {*}{
+    -body {
+        verdonk_test 1df4170f0fdecc E+124 +3981_499999999999999999& E+34
+    }
+    -result {}
+}
+test util-13.58 {just over half ulp - 6 digits} {*}{
+    -body {
+        verdonk_test 17e1e0f1c7a4ac E+415 +126300_5000000000000000000& E+120
+    }
+    -result {}
+}
+test util-13.59 {just over half ulp - 6 digits} {*}{
+    -body {
+        verdonk_test -1dda592e398dd7 E+418 -126300_5000000000000000000& E+121
+    }
+    -result {}
+}
+test util-13.60 {just under half ulp - 7 digits} {*}{
+    -body {
+        verdonk_test -1e597c0b94b7ae E+453 -4411845_499999999999999999& E+130
+    }
+    -result {}
+}
+test util-13.61 {just under half ulp - 9 digits} {*}{
+    -body {
+        verdonk_test 1c569e968e0944 E+427 +613850817_4999999999999999999& E+120
+    }
+    -result {}
+}
+test util-13.62 {just under half ulp - 9 digits} {*}{
+    -body {
+        verdonk_test -1c569e968e0944 E+428 -122770163_49999999999999999999& E+121
+    }
+    -result {}
+}
+test util-13.63 {just over half ulp - 18 digits} {*}{
+    -body {
+        verdonk_test 17ae0c186d8709 E+719 +408156062268363718_5000000000000000000& E+199
+    }
+    -result {}
+}
+test util-13.64 {just over exact - 1 digits} {*}{
+    -body {
+        verdonk_test 152d02c7e14af7 E+76 +1_0000000000000000& E+23
+    }
+    -result {}
+}
+test util-13.65 {just over exact - 1 digits} {*}{
+    -body {
+        verdonk_test -19d971e4fe8402 E+89 -1_0000000000000000& E+27
+    }
+    -result {}
+}
+test util-13.66 {just over exact - 1 digits} {*}{
+    -body {
+        verdonk_test 19d971e4fe8402 E+90 +2_0000000000000000& E+27
+    }
+    -result {}
+}
+test util-13.67 {just over exact - 1 digits} {*}{
+    -body {
+        verdonk_test -19d971e4fe8402 E+91 -4_0000000000000000& E+27
+    }
+    -result {}
+}
+test util-13.68 {just over exact - 1 digits} {*}{
+    -body {
+        verdonk_test 15798ee2308c3a E-27 +1_0000000000000000& E-8
+    }
+    -result {}
+}
+test util-13.69 {just over exact - 1 digits} {*}{
+    -body {
+        verdonk_test -15798ee2308c3a E-26 -2_0000000000000000& E-8
+    }
+    -result {}
+}
+test util-13.70 {just over exact - 1 digits} {*}{
+    -body {
+        verdonk_test 15798ee2308c3a E-25 +4_0000000000000000& E-8
+    }
+    -result {}
+}
+test util-13.71 {just over exact - 1 digits} {*}{
+    -body {
+        verdonk_test -1ef2d0f5da7dd9 E-84 -1_0000000000000000& E-25
+    }
+    -result {}
+}
+test util-13.72 {just under exact - 1 digits} {*}{
+    -body {
+        verdonk_test 1a784379d99db4 E+78 +4_9999999999999999& E+23
+    }
+    -result {}
+}
+test util-13.73 {just under exact - 1 digits} {*}{
+    -body {
+        verdonk_test -1a784379d99db4 E+80 -1_9999999999999999& E+24
+    }
+    -result {}
+}
+test util-13.74 {just under exact - 1 digits} {*}{
+    -body {
+        verdonk_test 13da329b633647 E+81 +2_9999999999999999& E+24
+    }
+    -result {}
+}
+test util-13.75 {just under exact - 1 digits} {*}{
+    -body {
+        verdonk_test -1cf389cd46047d E+85 -6_9999999999999999& E+25
+    }
+    -result {}
+}
+test util-13.76 {just under exact - 1 digits} {*}{
+    -body {
+        verdonk_test 19999999999999 E-3 +1_99999999999999999& E-1
+    }
+    -result {}
+}
+test util-13.77 {just under exact - 1 digits} {*}{
+    -body {
+        verdonk_test -13333333333333 E-2 -2_99999999999999999& E-1
+    }
+    -result {}
+}
+test util-13.78 {just under exact - 1 digits} {*}{
+    -body {
+        verdonk_test 16849b86a12b9b E-48 +4_99999999999999999& E-15
+    }
+    -result {}
+}
+test util-13.79 {just under exact - 1 digits} {*}{
+    -body {
+        verdonk_test -16849b86a12b9b E-46 -1_99999999999999999& E-14
+    }
+    -result {}
+}
+test util-13.80 {just over exact - 2 digits} {*}{
+    -body {
+        verdonk_test 17ccfc73126788 E-71 +63_00000000000000000& E-23
+    }
+    -result {}
+}
+test util-13.81 {just over exact - 2 digits} {*}{
+    -body {
+        verdonk_test -1dc03b8fd7016a E-68 -63_00000000000000000& E-22
+    }
+    -result {}
+}
+test util-13.82 {just under exact - 2 digits} {*}{
+    -body {
+        verdonk_test 13f7ced916872b E-5 +38_999999999999999999& E-3
+    }
+    -result {}
+}
+test util-13.83 {just over exact - 3 digits} {*}{
+    -body {
+        verdonk_test 1b297cad9f70b6 E+97 +269_000000000000000000& E+27
+    }
+    -result {}
+}
+test util-13.84 {just over exact - 3 digits} {*}{
+    -body {
+        verdonk_test -1b297cad9f70b6 E+98 -538_00000000000000000& E+27
+    }
+    -result {}
+}
+test util-13.85 {just over exact - 3 digits} {*}{
+    -body {
+        verdonk_test 1cdc06b20ef183 E-82 +373_00000000000000000& E-27
+    }
+    -result {}
+}
+test util-13.86 {just over exact - 4 digits} {*}{
+    -body {
+        verdonk_test 1b297cad9f70b6 E+96 +1345_00000000000000000& E+26
+    }
+    -result {}
+}
+# this one is not 4 digits, it is 3, and it is covered above.
+test util-13.87 {just over exact - 4 digits} {*}{
+    -constraints knownBadTest
+    -body {
+        verdonk_test -1b297cad9f70b6 E+97 -2690_00000000000000000& E+26
+    }
+    -result {}
+}
+test util-13.88 {just over exact - 5 digits} {*}{
+    -body {
+        verdonk_test -150a246ecd44f3 E-63 -14257_00000000000000000& E-23
+    }
+    -result {}
+}
+test util-13.89 {just under exact - 6 digits} {*}{
+    -body {
+        verdonk_test -119b96f36ec68b E-19 -209900_999999999999999999& E-11
+    }
+    -result {}
+}
+test util-13.90 {just over exact - 11 digits} {*}{
+    -body {
+        verdonk_test 1c06d366394441 E-35 +50980203373_000000000000000000& E-21
+    }
+    -result {}
+}
+test util-13.91 {just under exact - 12 digits} {*}{
+    -body {
+        verdonk_test -1f58ac4db68c90 E+122 -104166211810_99999999999999999& E+26
+    }
+    -result {}
+}
+test util-13.92 {just over half ulp - 1 digits} {*}{
+    -body {
+        verdonk_test 19d971e4fe8402 E+87 +2_5000000000000000& E+26
+    }
+    -result {}
+}
+test util-13.93 {just over half ulp - 1 digits} {*}{
+    -body {
+        verdonk_test -1dc74be914d16b E+81 -4_500000000000000& E+24
+    }
+    -result {}
+}
+test util-13.94 {just over half ulp - 1 digits} {*}{
+    -body {
+        verdonk_test 14adf4b7320335 E+84 +2_500000000000000& E+25
+    }
+    -result {}
+}
+test util-13.95 {just over half ulp - 1 digits} {*}{
+    -body {
+        verdonk_test -1ae22487c1042b E+85 -6_5000000000000000& E+25
+    }
+    -result {}
+}
+test util-13.96 {just over half ulp - 1 digits} {*}{
+    -body {
+        verdonk_test 187fe49aab41e0 E-54 +8_5000000000000000& E-17
+    }
+    -result {}
+}
+test util-13.97 {just over half ulp - 1 digits} {*}{
+    -body {
+        verdonk_test -1f5c05e4b23fd7 E-61 -8_5000000000000000& E-19
+    }
+    -result {}
+}
+test util-13.98 {just over half ulp - 1 digits} {*}{
+    -body {
+        verdonk_test 1faa7ab552a552 E-42 +4_5000000000000000& E-13
+    }
+    -result {}
+}
+test util-13.99 {just over half ulp - 1 digits} {*}{
+    -body {
+        verdonk_test -1b7cdfd9d7bdbb E-36 -2_5000000000000000& E-11
+    }
+    -result {}
+}
+test util-13.100 {just under half ulp - 1 digits} {*}{
+    -body {
+        verdonk_test 13da329b633647 E+80 +1_4999999999999999& E+24
+    }
+    -result {}
+}
+test util-13.101 {just under half ulp - 1 digits} {*}{
+    -body {
+        verdonk_test -1cf389cd46047d E+84 -3_49999999999999999& E+25
+    }
+    -result {}
+}
+test util-13.102 {just under half ulp - 1 digits} {*}{
+    -body {
+        verdonk_test 1f04ef12cb04cf E+85 +7_4999999999999999& E+25
+    }
+    -result {}
+}
+test util-13.103 {just under half ulp - 1 digits} {*}{
+    -body {
+        verdonk_test -1f04ef12cb04cf E+86 -1_4999999999999999& E+26
+    }
+    -result {}
+}
+test util-13.104 {just under half ulp - 1 digits} {*}{
+    -body {
+        verdonk_test 13333333333333 E-3 +1_49999999999999999& E-1
+    }
+    -result {}
+}
+test util-13.105 {just under half ulp - 1 digits} {*}{
+    -body {
+        verdonk_test -107e1fe91b0b70 E-36 -1_49999999999999999& E-11
+    }
+    -result {}
+}
+test util-13.106 {just under half ulp - 1 digits} {*}{
+    -body {
+        verdonk_test 149da7e361ce4c E-33 +1_49999999999999999& E-10
+    }
+    -result {}
+}
+test util-13.107 {just under half ulp - 1 digits} {*}{
+    -body {
+        verdonk_test -19c511dc3a41df E-30 -1_49999999999999999& E-9
+    }
+    -result {}
+}
+test util-13.108 {just over half ulp - 2 digits} {*}{
+    -body {
+        verdonk_test -1aa83d74267822 E+93 -16_5000000000000000& E+27
+    }
+    -result {}
+}
+test util-13.109 {just over half ulp - 2 digits} {*}{
+    -body {
+        verdonk_test 18f1d5969453de E+89 +96_5000000000000000& E+25
+    }
+    -result {}
+}
+test util-13.110 {just over half ulp - 2 digits} {*}{
+    -body {
+        verdonk_test 11d9bd564dcda6 E-70 +94_50000000000000000& E-23
+    }
+    -result {}
+}
+test util-13.111 {just over half ulp - 2 digits} {*}{
+    -body {
+        verdonk_test -1a58973ecbede6 E-48 -58_50000000000000000& E-16
+    }
+    -result {}
+}
+test util-13.112 {just over half ulp - 3 digits} {*}{
+    -body {
+        verdonk_test 1b297cad9f70b6 E+95 +672_50000000000000000& E+26
+    }
+    -result {}
+}
+test util-13.113 {just over half ulp - 3 digits} {*}{
+    -body {
+        verdonk_test -1b297cad9f70b6 E+96 -134_500000000000000000& E+27
+    }
+    -result {}
+}
+test util-13.114 {just over half ulp - 3 digits} {*}{
+    -body {
+        verdonk_test 1cdc06b20ef183 E-83 +186_50000000000000000& E-27
+    }
+    -result {}
+}
+test util-13.115 {just over half ulp - 3 digits} {*}{
+    -body {
+        verdonk_test -136071dcae4565 E-47 -860_50000000000000000& E-17
+    }
+    -result {}
+}
+test util-13.116 {just over half ulp - 6 digits} {*}{
+    -body {
+        verdonk_test 1cb968d297dde8 E+99 +113788_50000000000000000& E+25
+    }
+    -result {}
+}
+test util-13.117 {just over half ulp - 6 digits} {*}{
+    -body {
+        verdonk_test -11f3e1839eeab1 E+103 -113788_50000000000000000& E+26
+    }
+    -result {}
+}
+test util-13.118 {just under half ulp - 9 digits} {*}{
+    -body {
+        verdonk_test 1e9cec176c96f8 E+117 +317903333_49999999999999999& E+27
+    }
+    -result {}
+}
+test util-13.119 {just over half ulp - 11 digits} {*}{
+    -body {
+        verdonk_test 1c06d366394441 E-36 +25490101686_500000000000000000& E-21
+    }
+    -result {}
+}
+test util-13.120 {just under half ulp - 11 digits} {*}{
+    -body {
+        verdonk_test 1f58ac4db68c90 E+121 +52083105905_49999999999999999& E+26
+    }
+    -result {}
+}
+
+test util-14.1 {funky NaN} {*}{
+    -constraints {ieeeFloatingPoint && controversialNaN}
+    -body {
+	set ieeeValues(-NaN)
+    }
+    -result -NaN
+}
+
+test util-14.2 {funky NaN} {*}{
+    -constraints {ieeeFloatingPoint && controversialNaN}
+    -body {
+	set ieeeValues(-NaN(3456789abcdef))
+    }
+    -result -NaN(3456789abcdef)
+}
+
 # cleanup
 ::tcltest::cleanupTests
 return
+
+# Local Variables:
+# mode: tcl
+# End:
\ No newline at end of file
diff --git a/unix/Makefile.in b/unix/Makefile.in
index fad07e6..8eebff2 100644
--- a/unix/Makefile.in
+++ b/unix/Makefile.in
@@ -4,7 +4,7 @@
 # "./configure", which is a configuration script generated by the "autoconf"
 # program (constructs like "@foo@" will get replaced in the actual Makefile.
 #
-# RCS: @(#) $Id: Makefile.in,v 1.310 2010/11/04 15:55:45 dkf Exp $
+# RCS: @(#) $Id: Makefile.in,v 1.311 2010/11/28 23:20:11 kennykb Exp $
 
 VERSION 		= @TCL_VERSION@
 MAJOR_VERSION		= @TCL_MAJOR_VERSION@
@@ -322,12 +322,13 @@ TOMMATH_OBJS = bncore.o bn_reverse.o bn_fast_s_mp_mul_digs.o \
 	bn_mp_div_2d.o bn_mp_div_3.o \
         bn_mp_exch.o bn_mp_expt_d.o bn_mp_grow.o bn_mp_init.o \
 	bn_mp_init_copy.o bn_mp_init_multi.o bn_mp_init_set.o \
-	bn_mp_init_size.o bn_mp_karatsuba_mul.o \
+	bn_mp_init_set_int.o bn_mp_init_size.o bn_mp_karatsuba_mul.o \
 	bn_mp_karatsuba_sqr.o \
         bn_mp_lshd.o bn_mp_mod.o bn_mp_mod_2d.o bn_mp_mul.o bn_mp_mul_2.o \
         bn_mp_mul_2d.o bn_mp_mul_d.o bn_mp_neg.o bn_mp_or.o \
 	bn_mp_radix_size.o bn_mp_radix_smap.o \
-        bn_mp_read_radix.o bn_mp_rshd.o bn_mp_set.o bn_mp_shrink.o \
+        bn_mp_read_radix.o bn_mp_rshd.o bn_mp_set.o bn_mp_set_int.o \
+	bn_mp_shrink.o \
 	bn_mp_sqr.o bn_mp_sqrt.o bn_mp_sub.o bn_mp_sub_d.o \
         bn_mp_to_unsigned_bin.o bn_mp_to_unsigned_bin_n.o \
 	bn_mp_toom_mul.o bn_mp_toom_sqr.o bn_mp_toradix_n.o \
@@ -496,6 +497,7 @@ TOMMATH_SRCS = \
 	$(TOMMATH_DIR)/bn_mp_init_copy.c \
 	$(TOMMATH_DIR)/bn_mp_init_multi.c \
 	$(TOMMATH_DIR)/bn_mp_init_set.c \
+	$(TOMMATH_DIR)/bn_mp_init_set_int.c \
 	$(TOMMATH_DIR)/bn_mp_init_size.c \
 	$(TOMMATH_DIR)/bn_mp_karatsuba_mul.c \
 	$(TOMMATH_DIR)/bn_mp_karatsuba_sqr.c \
@@ -513,6 +515,7 @@ TOMMATH_SRCS = \
 	$(TOMMATH_DIR)/bn_mp_read_radix.c \
 	$(TOMMATH_DIR)/bn_mp_rshd.c \
 	$(TOMMATH_DIR)/bn_mp_set.c \
+	$(TOMMATH_DIR)/bn_mp_set_int.c \
 	$(TOMMATH_DIR)/bn_mp_shrink.c \
 	$(TOMMATH_DIR)/bn_mp_sqr.c \
 	$(TOMMATH_DIR)/bn_mp_sqrt.c \
@@ -1380,6 +1383,9 @@ bn_mp_init_multi.o: $(TOMMATH_DIR)/bn_mp_init_multi.c $(MATHHDRS)
 bn_mp_init_set.o: $(TOMMATH_DIR)/bn_mp_init_set.c $(MATHHDRS)
 	$(CC) -c $(CC_SWITCHES) $(TOMMATH_DIR)/bn_mp_init_set.c
 
+bn_mp_init_set_int.o: $(TOMMATH_DIR)/bn_mp_init_set_int.c $(MATHHDRS)
+	$(CC) -c $(CC_SWITCHES) $(TOMMATH_DIR)/bn_mp_init_set_int.c
+
 bn_mp_init_size.o:$(TOMMATH_DIR)/bn_mp_init_size.c $(MATHHDRS)
 	$(CC) -c $(CC_SWITCHES) $(TOMMATH_DIR)/bn_mp_init_size.c
 
@@ -1431,6 +1437,9 @@ bn_mp_rshd.o: $(TOMMATH_DIR)/bn_mp_rshd.c $(MATHHDRS)
 bn_mp_set.o: $(TOMMATH_DIR)/bn_mp_set.c $(MATHHDRS)
 	$(CC) -c $(CC_SWITCHES) $(TOMMATH_DIR)/bn_mp_set.c
 
+bn_mp_set_int.o: $(TOMMATH_DIR)/bn_mp_set_int.c $(MATHHDRS)
+	$(CC) -c $(CC_SWITCHES) $(TOMMATH_DIR)/bn_mp_set_int.c
+
 bn_mp_shrink.o: $(TOMMATH_DIR)/bn_mp_shrink.c $(MATHHDRS)
 	$(CC) -c $(CC_SWITCHES) $(TOMMATH_DIR)/bn_mp_shrink.c
 
diff --git a/win/Makefile.in b/win/Makefile.in
index 3097849..74d006c 100644
--- a/win/Makefile.in
+++ b/win/Makefile.in
@@ -4,7 +4,7 @@
 # "./configure", which is a configuration script generated by the "autoconf"
 # program (constructs like "@foo@" will get replaced in the actual Makefile.
 #
-# RCS: @(#) $Id: Makefile.in,v 1.185 2010/11/04 21:48:23 nijtmans Exp $
+# RCS: @(#) $Id: Makefile.in,v 1.186 2010/11/28 23:20:11 kennykb Exp $
 
 VERSION = @TCL_VERSION@
 
@@ -318,6 +318,7 @@ TOMMATH_OBJS = \
 	bn_mp_init_copy.${OBJEXT} \
 	bn_mp_init_multi.${OBJEXT} \
 	bn_mp_init_set.${OBJEXT} \
+	bn_mp_init_set_int.${OBJEXT} \
 	bn_mp_init_size.${OBJEXT} \
 	bn_mp_karatsuba_mul.${OBJEXT} \
 	bn_mp_karatsuba_sqr.$(OBJEXT) \
@@ -335,6 +336,7 @@ TOMMATH_OBJS = \
 	bn_mp_read_radix.${OBJEXT} \
 	bn_mp_rshd.${OBJEXT} \
 	bn_mp_set.${OBJEXT} \
+	bn_mp_set_int.${OBJEXT} \
 	bn_mp_shrink.${OBJEXT} \
 	bn_mp_sqr.${OBJEXT} \
 	bn_mp_sqrt.${OBJEXT} \
diff --git a/win/makefile.vc b/win/makefile.vc
index 9cc1a71..988d823 100644
--- a/win/makefile.vc
+++ b/win/makefile.vc
@@ -13,7 +13,7 @@
 # Copyright (c) 2003-2008 Pat Thoyts.
 #
 #------------------------------------------------------------------------------
-# RCS: @(#) $Id: makefile.vc,v 1.216 2010/11/04 21:48:23 nijtmans Exp $
+# RCS: @(#) $Id: makefile.vc,v 1.217 2010/11/28 23:20:11 kennykb Exp $
 #------------------------------------------------------------------------------
 
 # Check to see we are configured to build with MSVC (MSDEVDIR or MSVCDIR)
@@ -369,6 +369,7 @@ TOMMATHOBJS = \
 	$(TMP_DIR)\bn_mp_init_copy.obj \
 	$(TMP_DIR)\bn_mp_init_multi.obj \
 	$(TMP_DIR)\bn_mp_init_set.obj \
+	$(TMP_DIR)\bn_mp_init_set_int.obj \
 	$(TMP_DIR)\bn_mp_init_size.obj \
 	$(TMP_DIR)\bn_mp_karatsuba_mul.obj \
 	$(TMP_DIR)\bn_mp_karatsuba_sqr.obj \
@@ -386,6 +387,7 @@ TOMMATHOBJS = \
 	$(TMP_DIR)\bn_mp_read_radix.obj \
 	$(TMP_DIR)\bn_mp_rshd.obj \
 	$(TMP_DIR)\bn_mp_set.obj \
+	$(TMP_DIR)\bn_mp_set_int.obj \
 	$(TMP_DIR)\bn_mp_shrink.obj \
 	$(TMP_DIR)\bn_mp_sqr.obj \
 	$(TMP_DIR)\bn_mp_sqrt.obj \