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.

1 files changed, 2474 insertions, 298 deletions

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(