path: root/generic/tclStrToD.c

/*
 *----------------------------------------------------------------------
 *
 * tclStrToD.c --
 *
 *	This file contains a collection of procedures for managing conversions
 *	to/from floating-point in Tcl. They include TclParseNumber, which
 *	parses numbers from strings; TclDoubleDigits, which formats numbers
 *	into strings of digits, and procedures for interconversion among
 *	'double' and 'mp_int' types.
 *
 * Copyright (c) 2005 by Kevin B. Kenny.  All rights reserved.
 *
 * See the file "license.terms" for information on usage and redistribution of
 * this file, and for a DISCLAIMER OF ALL WARRANTIES.
 *
 * RCS: @(#) $Id: tclStrToD.c,v 1.18 2005/12/19 19:03:17 dgp Exp $
 *
 *----------------------------------------------------------------------
 */

#include <tclInt.h>
#include <stdio.h>
#include <stdlib.h>
#include <float.h>
#include <limits.h>
#include <math.h>
#include <ctype.h>
#include <tommath.h>

/*
 * Define TIP_114_FORMATS to accept 0b and 0o for binary and octal strings.
 * Define KILL_OCTAL as well as TIP_114_FORMATS to suppress interpretation of
 * numbers with leading zero as octal. (Ceterum censeo: numeros octonarios
 * delendos esse.)
 */

#define	TIP_114_FORMATS
#undef	KILL_OCTAL

#ifndef TIP_114_FORMATS
#undef	KILL_OCTAL
#endif

/*
 * This code supports (at least hypothetically), IBM, Cray, VAX and IEEE-754
 * floating point; of these, only IEEE-754 can represent NaN. IEEE-754 can be
 * uniquely determined by radix and by the widths of significand and exponent.
 */

#if (FLT_RADIX == 2) && (DBL_MANT_DIG == 53) && (DBL_MAX_EXP == 1024)
#   define IEEE_FLOATING_POINT
#endif

/*
 * gcc on x86 needs access to rounding controls, because of a questionable
 * feature where it retains intermediate results as IEEE 'long double' values
 * somewhat unpredictably. It is tempting to include fpu_control.h, but that
 * file exists only on Linux; it is missing on Cygwin and MinGW. Most gcc-isms
 * and ix86-isms are factored out here.
 */

#if defined(__GNUC__) && defined(__i386)
typedef unsigned int fpu_control_t __attribute__ ((__mode__ (__HI__)));
#define _FPU_GETCW(cw) __asm__ __volatile__ ("fnstcw %0" : "=m" (*&cw))
#define _FPU_SETCW(cw) __asm__ __volatile__ ("fldcw %0" : : "m" (*&cw))
#   define FPU_IEEE_ROUNDING	0x027f
#   define ADJUST_FPU_CONTROL_WORD
#endif

/*
 * HP's PA_RISC architecture uses 7ff4000000000000 to represent a quiet NaN.
 * Everyone else uses 7ff8000000000000. (Why, HP, why?)
 */

#ifdef __hppa
#   define NAN_START 0x7ff4
#   define NAN_MASK (((Tcl_WideUInt) 1) << 50)
#else
#   define NAN_START 0x7ff8
#   define NAN_MASK (((Tcl_WideUInt) 1) << 51)
#endif

/*
 * Constants used by this file (most of which are only ever calculated at
 * runtime).
 */

static int maxpow10_wide;	/* The powers of ten that can be represented
				 * exactly as wide integers. */
static Tcl_WideUInt *pow10_wide;
#define MAXPOW	22
static double pow10[MAXPOW+1];	/* The powers of ten that can be represented
				 * exactly as IEEE754 doubles. */
static int mmaxpow;		/* Largest power of ten that can be
				 * represented exactly in a 'double'. */
static int log10_DIGIT_MAX;	/* The number of decimal digits that fit in an
				 * mp_digit. */
static int log2FLT_RADIX;	/* Logarithm of the floating point radix. */
static int mantBits;		/* Number of bits in a double's significand */
static mp_int pow5[9];		/* Table of powers of 5**(2**n), up to
				 * 5**256 */
static double tiny;		/* The smallest representable double */
static int maxDigits;		/* The maximum number of digits to the left of
				 * the decimal point of a double. */
static int minDigits;		/* The maximum number of digits to the right
				 * of the decimal point in a double. */
static int mantDIGIT;		/* Number of mp_digit's needed to hold the
				 * significand of a double. */
static CONST double pow_10_2_n[] = {	/* Inexact higher powers of ten. */
    1.0,
    100.0,
    10000.0,
    1.0e+8,
    1.0e+16,
    1.0e+32,
    1.0e+64,
    1.0e+128,
    1.0e+256
};

/*
 * 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(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 double		Pow10TimesFrExp(int exponent, double fraction,
			    int *machexp);
static double		RefineApproximation(double approx,
			    mp_int *exactSignificand, int exponent);
static double		SafeLdExp(double fraction, int exponent);

/*
 *----------------------------------------------------------------------
 *
 * TclParseNumber --
 *
 *	Scans bytes, interpreted as characters in Tcl's internal encoding,
 *	and parses the longest prefix that is the string representation of
 *	a number in a format recognized by Tcl.
 *
 *	The arguments bytes, numBytes, and objPtr are the inputs which
 *	determine the string to be parsed.  If bytes is non-NULL, it
 *	points to the first byte to be scanned.  If bytes is NULL, then objPtr
 *	must be non-NULL, and the string representation of objPtr will be
 *	scanned (generated first, if necessary).  The numBytes argument
 *	determines the number of bytes to be scanned.  If numBytes is
 *	negative, the first NUL byte encountered will terminate the scan.
 *	If numBytes is non-negative, then no more than numBytes bytes will
 *	be scanned.  
 *
 *	The argument flags is an input that controls the numeric formats
 *	recognized by the parser.  The flag bits are:
 *
 *	- TCL_PARSE_INTEGER_ONLY:	accept only integer values; reject
 *		strings that denote floating point values (or accept only the
 *		leading portion of them that are integer values).
 *	- TCL_PARSE_SCAN_PREFIXES:	ignore the prefixes 0b and 0o that are
 *		not part of the [scan] command's vocabulary. Use only in
 *		combination with TCL_PARSE_INTEGER_ONLY.
 * 	- TCL_PARSE_OCTAL_ONLY:		parse only in the octal format, whether
 *		or not a prefix is present that would lead to octal parsing. Use
 *		only in combination with TCL_PARSE_INTEGER_ONLY.
 * 	- TCL_PARSE_HEXADECIMAL_ONLY:	parse only in the hexadecimal format,
 *		whether or not a prefix is present that would lead to
 *		hexadecimal parsing. Use only in combination with
 *		TCL_PARSE_INTEGER_ONLY.
 * 	- TCL_PARSE_DECIMAL_ONLY:	parse only in the decimal format, no
 *		matter whether a 0 prefix would normally force a different base.
 *	- TCL_PARSE_NO_WHITESPACE:	reject any leading/trailing whitespace
 *
 *	The arguments interp and expected are inputs that control error message
 * 	generation.  If interp is NULL, no error message will be generated.
 *	If interp is non-NULL, then expected must also be non-NULL.  When
 *	TCL_ERROR is returned, an error message will be left in the result
 *	of interp, and the expected argument will appear in the error message
 *	as the thing TclParseNumber expected, but failed to find in the string.
 *
 *	The arguments objPtr and endPtrPtr as well as the return code are the
 *	outputs.
 *
 *	When the parser cannot find any prefix of the string that matches a
 *	format it is looking for, TCL_ERROR is returned and an error message
 *	may be generated and returned as described above.  The contents of
 *	objPtr will not be changed.  If endPtrPtr is non-NULL, a pointer to
 *	the character in the string that terminated the scan will be written
 *	to *endPtrPtr.
 *
 *	When the parser determines that the entire string matches a format
 *	it is looking for, TCL_OK is returned, and if objPtr is non-NULL,
 *	then the internal rep and Tcl_ObjType of objPtr are set	to the
 *	"canonical" numeric value that matches the scanned string.  If
 *	endPtrPtr is non-NULL, a pointer to the end of the string will be
 *	written to *endPtrPtr (that is, either bytes+numBytes or a pointer
 *	to a terminating NUL byte).
 *
 *	When the parser determines that a partial string matches a format
 *	it is looking for, the value of endPtrPtr determines what happens:
 *
 *	- If endPtrPtr is NULL, then TCL_ERROR is returned, with error message
 *		generation as above.
 *
 *	- If endPtrPtr is non-NULL, then TCL_OK is returned and objPtr
 *		internals are set as above.  Also, a pointer to the first
 *		character following the parsed numeric string is written
 *		to *endPtrPtr.
 *
 *	In some cases where the string being scanned is the string rep of
 *	objPtr, this routine can leave objPtr in an inconsistent state
 *	where its string rep and its internal rep do not agree.  In these
 *	cases the internal rep will be in agreement with only some substring
 *	of the string rep.  This might happen if the caller passes in a
 * 	non-NULL bytes value that points somewhere into the string rep.  It
 *	might happen if the caller passes in a numBytes value that limits the 
 * 	scan to only a prefix of the string rep.  Or it might happen if a
 *	non-NULL value of endPtrPtr permits a TCL_OK return from only a partial
 *	string match.  It is the responsibility of the caller to detect and
 *	correct such inconsistencies when they can and do arise.
 *
 * Results:
 *	Returns a standard Tcl result.
 *
 * Side effects:
 *	The string representaton of objPtr may be generated.
 *	
 *	The internal representation and Tcl_ObjType of objPtr may be changed.
 *	This may involve allocation and/or freeing of memory.
 *
 *----------------------------------------------------------------------
 */

int
TclParseNumber(
    Tcl_Interp *interp,		/* Used for error reporting. May be NULL */
    Tcl_Obj *objPtr,		/* Object to receive the internal rep */
    CONST char *expected,	/* Description of the type of number the caller
				 * expects to be able to parse ("integer", 
				 * "boolean value", etc.). */
    CONST char *bytes,		/* Pointer to the start of the string to scan */
    int numBytes,		/* Maximum number of bytes to scan, see above */
    CONST char **endPtrPtr,	/* Place to store pointer to the character
				 * that terminated the scan */
    int flags)			/* Flags governing the parse */
{
    enum State {
	INITIAL, SIGNUM, ZERO, ZERO_X, 
#ifdef TIP_114_FORMATS
	ZERO_O, ZERO_B, BINARY,
#endif
	HEXADECIMAL, OCTAL, BAD_OCTAL, DECIMAL, 
	LEADING_RADIX_POINT, FRACTION,
	EXPONENT_START, EXPONENT_SIGNUM, EXPONENT,
	sI, sIN, sINF, sINFI, sINFIN, sINFINI, sINFINIT, sINFINITY
#ifdef IEEE_FLOATING_POINT
	, sN, sNA, sNAN, sNANPAREN, sNANHEX, sNANFINISH
#endif
    } state = INITIAL;
    enum State acceptState = INITIAL;

    int signum = 0;		/* Sign of the number being parsed */
    Tcl_WideUInt significandWide = 0;
				/* Significand of the number being parsed (if
				 * no overflow) */
    mp_int significandBig;	/* Significand of the number being parsed (if
				 * it overflows significandWide) */
    int significandOverflow = 0;/* Flag==1 iff significandBig is used */
    Tcl_WideUInt octalSignificandWide = 0;
				/* Significand of an octal number; needed
				 * because we don't know whether a number with
				 * a leading zero is octal or decimal until
				 * we've scanned forward to a '.' or 'e' */
    mp_int octalSignificandBig;	/* Significand of octal number once
				 * octalSignificandWide overflows */
    int octalSignificandOverflow = 0;
				/* Flag==1 if octalSignificandBig is used */
    int numSigDigs = 0;		/* Number of significant digits in the decimal
				 * significand */
    int numTrailZeros = 0;	/* Number of trailing zeroes at the current
				 * point in the parse. */
    int numDigitsAfterDp = 0;	/* Number of digits scanned after the decimal
				 * point */
    int exponentSignum = 0;	/* Signum of the exponent of a floating point
				 * number */
    long exponent = 0;		/* Exponent of a floating point number */
    CONST char *p;		/* Pointer to next character to scan */
    size_t len;			/* Number of characters remaining after p */
    CONST char *acceptPoint;	/* Pointer to position after last character in
				 * an acceptable number */
    size_t acceptLen;		/* Number of characters following that
				 * point. */
    int status = TCL_OK;	/* Status to return to caller */
    char d = 0;			/* Last hexadecimal digit scanned; initialized
				 * to avoid a compiler warning. */
    int shift = 0;		/* Amount to shift when accumulating binary */
#ifdef TIP_114_FORMATS
    int explicitOctal = 0;
#endif

#define ALL_BITS	(~(Tcl_WideUInt)0)
#define MOST_BITS	(ALL_BITS >> 1)

    /* 
     * Initialize bytes to start of the object's string rep if the caller
     * didn't pass anything else.
     */

    if (bytes == NULL) {
	bytes = TclGetString(objPtr);
    }

    p = bytes;
    len = numBytes;
    acceptPoint = p;
    acceptLen = len;
    while (1) {
	char c = len ? *p : '\0';
	switch (state) {

	case INITIAL:
	    /*
	     * Initial state. Acceptable characters are +, -, digits, period,
	     * I, N, and whitespace.
	     */

	    if (isspace(UCHAR(c))) {
		if (flags & TCL_PARSE_NO_WHITESPACE) {
		    goto endgame;
		}
		break;
	    } else if (c == '+') {
		state = SIGNUM;
		break;
	    } else if (c == '-') {
		signum = 1;
		state = SIGNUM;
		break;
	    } 
	    /* FALLTHROUGH */
	    
	case SIGNUM:
	    /*
	     * Scanned a leading + or -. Acceptable characters are digits,
	     * period, I, and N.
	     */

	    if (c == '0') {
		if (flags & TCL_PARSE_DECIMAL_ONLY) {
		    state = DECIMAL;
		} else {
		    state = ZERO;
		}
		break;
	    } else if (flags & TCL_PARSE_HEXADECIMAL_ONLY) {
		goto zerox;
	    } else if (flags & TCL_PARSE_OCTAL_ONLY) {
		goto zeroo;
	    } else if (isdigit(UCHAR(c))) {
		significandWide = c - '0';
		numSigDigs = 1;
		state = DECIMAL;
		break;
	    } else if (flags & TCL_PARSE_INTEGER_ONLY) {
		goto endgame;
	    } else if (c == '.') {
		state = LEADING_RADIX_POINT;
		break;
	    } else if (c == 'I' || c == 'i') {
		state = sI;
		break;
#ifdef IEEE_FLOATING_POINT
	    } else if (c == 'N' || c == 'n') {
		state = sN;
		break;
#endif
	    }
	    goto endgame;

	case ZERO:
	    /*
	     * Scanned a leading zero (perhaps with a + or -). Acceptable
	     * inputs are digits, period, X, and E. If 8 or 9 is encountered,
	     * the number can't be octal. This state and the OCTAL state
	     * differ only in whether they recognize 'X'.
	     */

	    acceptState = state;
	    acceptPoint = p;
	    acceptLen = len;
	    if (c == 'x' || c == 'X') {
		state = ZERO_X;
		break;
	    }
	    if (flags & TCL_PARSE_HEXADECIMAL_ONLY) {
		goto zerox;
	    }
#ifdef TIP_114_FORMATS 
	    if (flags & TCL_PARSE_SCAN_PREFIXES) {
		goto zeroo;
	    }
	    if (c == 'b' || c == 'B') {
		state = ZERO_B;
		break;
	    }
	    if (c == 'o' || c == 'O') {
		explicitOctal = 1;
		state = ZERO_O;
		break;
	    }
#ifdef KILL_OCTAL
	    goto decimal;
#endif
#endif
	    /* FALLTHROUGH */

	case OCTAL:
	    /*
	     * Scanned an optional + or -, followed by a string of octal
	     * digits. Acceptable inputs are more digits, period, or E. If 8
	     * or 9 is encountered, commit to floating point.
	     */

	    acceptState = state;
	    acceptPoint = p;
	    acceptLen = len;
#ifdef TIP_114_FORMATS
	    /* FALLTHROUGH */
	case ZERO_O:
#endif
	zeroo:
	    if (c == '0') {
		++numTrailZeros;
		state = OCTAL;
		break;
	    } else if (c >= '1' && c <= '7') {
		if (objPtr != NULL) {
		    shift = 3 * (numTrailZeros + 1);
		    significandOverflow = AccumulateDecimalDigit(
			    (unsigned)(c-'0'), numTrailZeros,
			    &significandWide, &significandBig,
			    significandOverflow);

		    if (!octalSignificandOverflow) {
			/*
			 * Shifting by more bits than are in the value being
			 * shifted is at least de facto nonportable. Check for
			 * too large shifts first.
			 */

			if ((octalSignificandWide != 0) 
				&& ((shift >= CHAR_BIT*sizeof(Tcl_WideUInt)) 
				|| (octalSignificandWide 
				> (~(Tcl_WideUInt)0 >> shift)))) {
			    octalSignificandOverflow = 1;
			    TclBNInitBignumFromWideUInt(&octalSignificandBig,
				    octalSignificandWide);
			}
		    }
		    if (!octalSignificandOverflow) {
			octalSignificandWide =
				(octalSignificandWide << shift) + (c - '0');
		    } else {
			mp_mul_2d(&octalSignificandBig, shift,
				&octalSignificandBig);
			mp_add_d(&octalSignificandBig, (mp_digit)(c - '0'),
				&octalSignificandBig);
		    }
		}
		if (numSigDigs != 0) {
		    numSigDigs += numTrailZeros+1;
		} else {
		    numSigDigs = 1;
		}
		numTrailZeros = 0;
		state = OCTAL;
		break;
	    }
	    /* FALLTHROUGH */

	case BAD_OCTAL:
#ifdef TIP_114_FORMATS
	    if (explicitOctal) {
		/*
		 * No forgiveness for bad digits in explicitly octal
		 * numbers.
		 */

		goto endgame;
	    }
#endif
	    if (flags & TCL_PARSE_INTEGER_ONLY) {
		/*
		 * No seeking floating point when parsing only integer.
		 */

		goto endgame;
	    }
#ifndef KILL_OCTAL

	    /*
	     * Scanned a number with a leading zero that contains an 8, 9,
	     * radix point or E. This is an invalid octal number, but might
	     * still be floating point.
	     */

	    if (c == '0') {
		++numTrailZeros;
		state = BAD_OCTAL;
		break;
	    } else if (isdigit(UCHAR(c))) {
		if (objPtr != NULL) {
		    significandOverflow = AccumulateDecimalDigit(
			    (unsigned)(c-'0'), numTrailZeros,
			    &significandWide, &significandBig,
			    significandOverflow);
		}
		if (numSigDigs != 0) {
		    numSigDigs += (numTrailZeros + 1);
		} else {
		    numSigDigs = 1;
		}
		numTrailZeros = 0;
		state = BAD_OCTAL;
		break;
	    } else if (c == '.') {
		state = FRACTION;
		break;
	    } else if (c == 'E' || c == 'e') {
		state = EXPONENT_START;
		break;
	    } 
#endif
	    goto endgame;

	    /*
	     * Scanned 0x. If state is HEXADECIMAL, scanned at least one
	     * character following the 0x. The only acceptable inputs are
	     * hexadecimal digits.
	     */

	case HEXADECIMAL:
	    acceptState = state;
	    acceptPoint = p;
	    acceptLen = len;
	    /* FALLTHROUGH */

	case ZERO_X:
	zerox:
	    if (c == '0') {
		++numTrailZeros;
		state = HEXADECIMAL;
		break;
	    } else if (isdigit(UCHAR(c))) {
		d = (c-'0');
	    } else if (c >= 'A' && c <= 'F') {
		d = (c-'A'+10);
	    } else if (c >= 'a' && c <= 'f') {
		d = (c-'a'+10);
	    } else {
		goto endgame;
	    }
	    if (objPtr != NULL) {
		shift = 4 * (numTrailZeros + 1);
		if (!significandOverflow) {
		    /*
		     * Shifting by more bits than are in the value being
		     * shifted is at least de facto nonportable. Check for too
		     * large shifts first.
		     */

		    if (significandWide != 0 &&
			    (shift >= CHAR_BIT*sizeof(Tcl_WideUInt) ||
			    significandWide > (~(Tcl_WideUInt)0 >> shift))) {
			significandOverflow = 1;
			TclBNInitBignumFromWideUInt(&significandBig,
				significandWide);
		    }
		}
		if (!significandOverflow) {
		    significandWide = (significandWide << shift) + d;
		} else {
		    mp_mul_2d(&significandBig, shift, &significandBig);
		    mp_add_d(&significandBig, (mp_digit) d, &significandBig);
		}
	    }
	    numTrailZeros = 0;
	    state = HEXADECIMAL;
	    break;

#ifdef TIP_114_FORMATS
	case BINARY:
	    acceptState = state;
	    acceptPoint = p;
	    acceptLen = len;
	case ZERO_B:
	    if (c == '0') {
		++numTrailZeros;
		state = BINARY;
		break;
	    } else if (c != '1') {
		goto endgame;
	    }
	    if (objPtr != NULL) {
		shift = numTrailZeros + 1;
		if (!significandOverflow) {
		    /*
		     * Shifting by more bits than are in the value being
		     * shifted is at least de facto nonportable. Check for too
		     * large shifts first.
		     */

		    if (significandWide != 0 &&
			    (shift >= CHAR_BIT*sizeof(Tcl_WideUInt) ||
			    significandWide > (~(Tcl_WideUInt)0 >> shift))) {
			significandOverflow = 1;
			TclBNInitBignumFromWideUInt(&significandBig,
				significandWide);
		    }
		}
		if (!significandOverflow) {
		    significandWide = (significandWide << shift) + 1;
		} else {
		    mp_mul_2d(&significandBig, shift, &significandBig);
		    mp_add_d(&significandBig, (mp_digit) 1, &significandBig);
		}
	    }
	    numTrailZeros = 0;
	    state = BINARY;
	    break;
#endif

	case DECIMAL:
	    /*
	     * Scanned an optional + or - followed by a string of decimal
	     * digits.
	     */

#ifdef KILL_OCTAL
	decimal:
#endif
	    acceptState = state;
	    acceptPoint = p;
	    acceptLen = len;
	    if (c == '0') {
		++numTrailZeros;
		state = DECIMAL;
		break;
	    } else if (isdigit(UCHAR(c))) {
		if (objPtr != NULL) {
		    significandOverflow = AccumulateDecimalDigit(
			    (unsigned)(c - '0'), numTrailZeros, 
			    &significandWide, &significandBig, 
			    significandOverflow);
		}
		numSigDigs += numTrailZeros+1;
		numTrailZeros = 0;
		state = DECIMAL;
		break;
	    } else if (flags & TCL_PARSE_INTEGER_ONLY) {
		goto endgame;
	    } else if (c == '.') {
		state = FRACTION;
		break;
	    } else if (c == 'E' || c == 'e') {
		state = EXPONENT_START;
		break;
	    }
	    goto endgame;

	    /* 
	     * Found a decimal point. If no digits have yet been scanned, E is
	     * not allowed; otherwise, it introduces the exponent. If at least
	     * one digit has been found, we have a possible complete number.
	     */

	case FRACTION:
	    acceptState = state;
	    acceptPoint = p;
	    acceptLen = len;
	    if (c == 'E' || c=='e') {
		state = EXPONENT_START;
		break;
	    }
	    /* FALLTHROUGH */

	case LEADING_RADIX_POINT:
	    if (c == '0') {
		++numDigitsAfterDp;
		++numTrailZeros;
		state = FRACTION;
		break;
	    } else if (isdigit(UCHAR(c))) {
		++numDigitsAfterDp;
		if (objPtr != NULL) {
		    significandOverflow = AccumulateDecimalDigit(
			    (unsigned)(c-'0'), numTrailZeros, 
			    &significandWide, &significandBig, 
			    significandOverflow);
		}
		if (numSigDigs != 0) {
		    numSigDigs += numTrailZeros+1;
		} else {
		    numSigDigs = 1;
		}
		numTrailZeros = 0;
		state = FRACTION;
		break;
	    }
	    goto endgame;

	case EXPONENT_START:
	    /* 
	     * Scanned the E at the start of an exponent. Make sure a legal
	     * character follows before using the C library strtol routine,
	     * which allows whitespace.
	     */

	    if (c == '+') {
		state = EXPONENT_SIGNUM;
		break;
	    } else if (c == '-') {
		exponentSignum = 1;
		state = EXPONENT_SIGNUM;
		break;
	    }
	    /* FALLTHROUGH */

	case EXPONENT_SIGNUM:
	    /*
	     * Found the E at the start of the exponent, followed by a sign
	     * character.
	     */

	    if (isdigit(UCHAR(c))) {
		exponent = c - '0';
		state = EXPONENT;
		break;
	    }
	    goto endgame;

	case EXPONENT:
	    /* 
	     * Found an exponent with at least one digit. Accumulate it,
	     * making sure to hard-pin it to LONG_MAX on overflow.
	     */

	    acceptState = state;
	    acceptPoint = p;
	    acceptLen = len;
	    if (isdigit(UCHAR(c))) {
		if (exponent < (LONG_MAX - 9) / 10) {
		    exponent = 10 * exponent + (c - '0');
		} else {
		    exponent = LONG_MAX;
		}
		state = EXPONENT;
		break;
	    }
	    goto endgame;

	    /*
	     * Parse out INFINITY by simply spelling it out. INF is accepted
	     * as an abbreviation; other prefices are not.
	     */

	case sI:
	    if (c == 'n' || c == 'N') {
		state = sIN;
		break;
	    } 
	    goto endgame;
	case sIN:
	    if (c == 'f' || c == 'F') {
		state = sINF;
		break;
	    } 
	    goto endgame;
	case sINF:
	    acceptState = state;
	    acceptPoint = p;
	    acceptLen = len;
	    if (c == 'i' || c == 'I') {
		state = sINFI;
		break;
	    }
	    goto endgame;
	case sINFI:
	    if (c == 'n' || c == 'N') {
		state = sINFIN;
		break;
	    }
	    goto endgame;
	case sINFIN:
	    if (c == 'i' || c == 'I') {
		state = sINFINI;
		break;
	    }
	    goto endgame;
	case sINFINI:
	    if (c == 't' || c == 'T') {
		state = sINFINIT;
		break;
	    }
	    goto endgame;
	case sINFINIT:
	    if (c == 'y' || c == 'Y') {
		state = sINFINITY;
		break;
	    }
	    goto endgame;

	    /*
	     * Parse NaN's.
	     */
#ifdef IEEE_FLOATING_POINT
	case sN:
	    if (c == 'a' || c == 'A') {
		state = sNA;
		break;
	    }
	    goto endgame;
	case sNA:
	    if (c == 'n' || c == 'N') {
		state = sNAN;
		break;
	    }
	    goto endgame;
	case sNAN:
	    acceptState = state;
	    acceptPoint = p;
	    acceptLen = len;
	    if (c == '(') {
		state = sNANPAREN;
		break;
	    }
	    goto endgame;

	    /*
	     * Parse NaN(hexdigits)
	     */
	case sNANHEX:
	    if (c == ')') {
		state = sNANFINISH;
		break;
	    }
	    /* FALLTHROUGH */
	case sNANPAREN:
	    if (isspace(UCHAR(c))) {
		break;
	    }
	    if (numSigDigs < 13) {
		if (c >= '0' && c <= '9') {
		    d = c - '0';
		} else if (c >= 'a' && c <= 'f') {
		    d = 10 + c - 'a';
		} else if (c >= 'A' && c <= 'F') {
		    d = 10 + c - 'A';
		}
		significandWide = (significandWide << 4) + d;
		state = sNANHEX;
		break;
	    }
	    goto endgame;
	case sNANFINISH:
#endif

	case sINFINITY:
	    acceptState = state;
	    acceptPoint = p;
	    acceptLen = len;
	    goto endgame;
	}
	++p; 
	--len;
    }

  endgame:
    if (acceptState == INITIAL) {
	/* No numeric string at all found */
	status = TCL_ERROR;
	if (endPtrPtr != NULL) {
	    *endPtrPtr = p;
	}
    } else {
	/* Back up to the last accepting state in the lexer. */
	p = acceptPoint;
	len = acceptLen;
	if (!(flags & TCL_PARSE_NO_WHITESPACE)) {
	    /* Accept trailing whitespace */
	    while (len != 0 && isspace(UCHAR(*p))) {
		++p;
		--len;
	    }
	}
	if (endPtrPtr == NULL) {
	    if ((len != 0) && ((numBytes > 0) || (*p != '\0'))) {
		status = TCL_ERROR;
	    }
	} else {
	    *endPtrPtr = p;
	}
    }

    /*
     * Generate and store the appropriate internal rep.
     */

    if (status == TCL_OK && objPtr != NULL) {
	TclFreeIntRep(objPtr);
	switch (acceptState) {
	case SIGNUM:
	case BAD_OCTAL:
	case ZERO_X:
#ifdef TIP_114_FORMATS
	case ZERO_O:
	case ZERO_B:
#endif
	case LEADING_RADIX_POINT:
	case EXPONENT_START:
	case EXPONENT_SIGNUM:
	case sI:
	case sIN:
	case sINFI:
	case sINFIN:
	case sINFINI:
	case sINFINIT:
	case sN:
	case sNA:
	case sNANPAREN:
	case sNANHEX:
	    Tcl_Panic("TclParseNumber: bad acceptState %d parsing '%s'",
		    acceptState, bytes);

#ifdef TIP_114_FORMATS
	case BINARY:
	    shift = numTrailZeros;
	    if (!significandOverflow && significandWide != 0 &&
		    (shift >= CHAR_BIT*sizeof(Tcl_WideUInt) ||
		    significandWide > (MOST_BITS + signum) >> shift)) {
		significandOverflow = 1;
		TclBNInitBignumFromWideUInt(&significandBig, significandWide);
	    }
	    if (shift) {
		if (!significandOverflow) {
		    significandWide <<= shift;
		} else {
		    mp_mul_2d(&significandBig, shift, &significandBig);
		}
	    }
	    goto returnInteger;
#endif

	case HEXADECIMAL:
	    /*
	     * Returning a hex integer. Final scaling step.
	     */

	    shift = 4 * numTrailZeros;
	    if (!significandOverflow && significandWide !=0 &&
		    (shift >= CHAR_BIT*sizeof(Tcl_WideUInt) ||
		    significandWide > (MOST_BITS + signum) >> shift)) {
		significandOverflow = 1;
		TclBNInitBignumFromWideUInt(&significandBig, significandWide);
	    }
	    if (shift) {
		if (!significandOverflow) {
		    significandWide <<= shift;
		} else {
		    mp_mul_2d(&significandBig, shift, &significandBig);
		}
	    }
	    goto returnInteger;

	case OCTAL:
	    /*
	     * Returning an octal integer. Final scaling step
	     */

	    shift = 3 * numTrailZeros;
	    if (!octalSignificandOverflow && octalSignificandWide != 0 &&
		    (shift >= CHAR_BIT*sizeof(Tcl_WideUInt) ||
		    octalSignificandWide > (MOST_BITS + signum) >> shift)) {
		octalSignificandOverflow = 1;
		TclBNInitBignumFromWideUInt(&octalSignificandBig,
			octalSignificandWide);
	    }
	    if (shift) {
		if (!octalSignificandOverflow) {
		    octalSignificandWide <<= shift;
		} else {
		    mp_mul_2d(&octalSignificandBig, shift,
			    &octalSignificandBig);
		}
	    }
	    if (!octalSignificandOverflow) {
		if (octalSignificandWide > 
			(Tcl_WideUInt)(((~(unsigned long)0) >> 1) + signum)) {
#ifndef NO_WIDE_TYPE
		    if (octalSignificandWide <= (MOST_BITS + signum)) {
			objPtr->typePtr = &tclWideIntType;
			if (signum) {
			    objPtr->internalRep.wideValue =
				    - (Tcl_WideInt) octalSignificandWide;
			} else {
			    objPtr->internalRep.wideValue =
				    (Tcl_WideInt) octalSignificandWide;
			}
			break;
		    }
#endif
		    TclBNInitBignumFromWideUInt(&octalSignificandBig,
			    octalSignificandWide);
		    octalSignificandOverflow = 1;
		} else {
		    objPtr->typePtr = &tclIntType;
		    if (signum) {
			objPtr->internalRep.longValue =
				- (long) octalSignificandWide;
		    } else {
			objPtr->internalRep.longValue =
				(long) octalSignificandWide;
		    }
		}
	    }
	    if (octalSignificandOverflow) {
		if (signum) {
		    mp_neg(&octalSignificandBig, &octalSignificandBig);
		}
		TclSetBignumIntRep(objPtr, &octalSignificandBig);
	    }		
	    break;

	case ZERO:
	case DECIMAL:
	    significandOverflow = AccumulateDecimalDigit(0, numTrailZeros-1,
		    &significandWide, &significandBig, significandOverflow);
	    if (!significandOverflow && (significandWide > MOST_BITS+signum)) {
		significandOverflow = 1;
		TclBNInitBignumFromWideUInt(&significandBig, significandWide);
	    }
	returnInteger:
	    if (!significandOverflow) {
		if (significandWide > 
			(Tcl_WideUInt)(((~(unsigned long)0) >> 1) + signum)) {
#ifndef NO_WIDE_TYPE
		    if (significandWide <= MOST_BITS+signum) {
			objPtr->typePtr = &tclWideIntType;
			if (signum) {
			    objPtr->internalRep.wideValue =
				    - (Tcl_WideInt) significandWide;
			} else {
			    objPtr->internalRep.wideValue =
				    (Tcl_WideInt) significandWide;
			}
			break;
		    }
#endif
		    TclBNInitBignumFromWideUInt(&significandBig,
			    significandWide);
		    significandOverflow = 1;
		} else {
		    objPtr->typePtr = &tclIntType;
		    if (signum) {
			objPtr->internalRep.longValue =
				- (long) significandWide;
		    } else {
			objPtr->internalRep.longValue =
				(long) significandWide;
		    }
		}
	    }
	    if (significandOverflow) {
		if (signum) {
		    mp_neg(&significandBig, &significandBig);
		}
		TclSetBignumIntRep(objPtr, &significandBig);
	    }
	    break;

	case FRACTION:
	case EXPONENT:

	    /*
	     * Here, we're parsing a floating-point number. 'significandWide'
	     * or 'significandBig' contains the exact significand, according
	     * to whether 'significandOverflow' is set. The desired floating
	     * point value is significand * 10**k, where
	     * k = numTrailZeros+exponent-numDigitsAfterDp.
	     */

	    objPtr->typePtr = &tclDoubleType;
	    if (exponentSignum) {
		exponent = - exponent;
	    }
	    if (!significandOverflow) {
		objPtr->internalRep.doubleValue = MakeLowPrecisionDouble(
			signum, significandWide, numSigDigs,
			(numTrailZeros + exponent - numDigitsAfterDp));
	    } else {
		objPtr->internalRep.doubleValue = MakeHighPrecisionDouble(
			signum, &significandBig, numSigDigs,
			(numTrailZeros + exponent - numDigitsAfterDp));
	    }
	    break;

	case sINF:
	case sINFINITY:
	    if (signum) {
		objPtr->internalRep.doubleValue = -HUGE_VAL;
	    } else {
		objPtr->internalRep.doubleValue = HUGE_VAL;
	    }
	    objPtr->typePtr = &tclDoubleType;
	    break;

	case sNAN:
	case sNANFINISH:
	    objPtr->internalRep.doubleValue = MakeNaN(signum, significandWide);
	    objPtr->typePtr = &tclDoubleType;
	    break;

	case INITIAL:
	    /* This case only to silence compiler warning */
	    Tcl_Panic("TclParseNumber: state INITIAL can't happen here");
	}
    }

    /*
     * Format an error message when an invalid number is encountered.
     */

    if (status != TCL_OK) {
	if (interp != NULL) {
	    Tcl_Obj *msg = Tcl_NewStringObj("expected ", -1);
	    Tcl_AppendToObj(msg, expected, -1);
	    Tcl_AppendToObj(msg, " but got \"", -1);
	    TclAppendLimitedToObj(msg, bytes, numBytes, 50, "");
	    Tcl_AppendToObj(msg, "\"", -1);
	    if (state == BAD_OCTAL) {
		Tcl_AppendToObj(msg, " (looks like invalid octal number)", -1);
	    }
	    Tcl_SetObjResult(interp, msg);
	}
    }

    /*
     * Free memory.
     */

    if (octalSignificandOverflow) {
	mp_clear(&octalSignificandBig);
    }
    if (significandOverflow) {
	mp_clear(&significandBig);
    }
    return status;
}

/*
 *----------------------------------------------------------------------
 *
 * AccumulateDecimalDigit --
 *
 *	Consume a decimal digit in a number being scanned.
 *
 * Results:
 *	Returns 1 if the number has overflowed to a bignum, 0 if it still fits
 *	in a wide integer.
 *
 * Side effects:
 *	Updates either the wide or bignum representation.
 *
 *----------------------------------------------------------------------
 */

static int
AccumulateDecimalDigit(
    unsigned digit,		/* Digit being scanned. */
    int numZeros,		/* Count of zero digits preceding the digit
				 * being scanned. */
    Tcl_WideUInt *wideRepPtr,	/* Representation of the partial number as a
				 * wide integer. */
    mp_int *bignumRepPtr,	/* Representation of the partial number as a
				 * bignum. */
    int bignumFlag)		/* Flag == 1 if the number overflowed previous
				 * to this digit. */
{
    int i, n;

    /*
     * Check if the number still fits in a wide.
     */

    if (!bignumFlag && *wideRepPtr!=0 && ((numZeros >= maxpow10_wide) ||
	    *wideRepPtr > ((~(Tcl_WideUInt)0)-digit)/pow10_wide[numZeros+1])) {
	/*
	 * Oops, it's overflowed, have to allocate a bignum.
	 */

	TclBNInitBignumFromWideUInt (bignumRepPtr, *wideRepPtr);
	bignumFlag = 1;
    }

    /*
     * Multiply the number by 10**numZeros+1 and add in the new digit.
     */

    if (!bignumFlag) {
	/*
	 * Wide multiplication.
	 */

	*wideRepPtr = *wideRepPtr * pow10_wide[numZeros+1] + digit;
    } else if (numZeros < log10_DIGIT_MAX) {
	/*
	 * Up to about 8 zeros - single digit multiplication.
	 */

	mp_mul_d(bignumRepPtr, (mp_digit) pow10_wide[numZeros+1],
		bignumRepPtr);
	mp_add_d(bignumRepPtr, (mp_digit) digit, bignumRepPtr);
    } else {
	/* 
	 * More than single digit multiplication. Multiply by the appropriate
	 * small powers of 5, and then shift. Large strings of zeroes are
	 * eaten 256 at a time; this is less efficient than it could be, but
	 * seems implausible. We presume that DIGIT_BIT is at least 27. The
	 * first multiplication, by up to 10**7, is done with a one-DIGIT
	 * multiply (this presumes that DIGIT_BIT >= 24).
	 */

	n = numZeros + 1;
	mp_mul_d(bignumRepPtr, (mp_digit) pow10_wide[n&0x7], bignumRepPtr);
	for (i=3; i<=7; ++i) {
	    if (n & (1 << i)) {
		mp_mul(bignumRepPtr, pow5+i, bignumRepPtr);
	    }
	}
	while (n >= 256) {
	    mp_mul(bignumRepPtr, pow5+8, bignumRepPtr);
	    n -= 256;
	}
	mp_mul_2d(bignumRepPtr, (int)(numZeros+1)&~0x7, bignumRepPtr);
    }

    return bignumFlag;
}

/*
 *----------------------------------------------------------------------
 *
 * MakeLowPrecisionDouble --
 *
 *	Makes the double precision number, signum*significand*10**exponent.
 *
 * Results:
 *	Returns the constructed number.
 *
 *	Common cases, where there are few enough digits that the number can be
 *	represented with at most roundoff, are handled specially here. If the
 *	number requires more than one rounded operation to compute, the code
 *	promotes the significand to a bignum and calls MakeHighPrecisionDouble
 *	to do it instead.
 *
 *----------------------------------------------------------------------
 */

static double
MakeLowPrecisionDouble(
    int signum,			/* 1 if the number is negative, 0 otherwise */
    Tcl_WideUInt significand,	/* Significand of the number */
    int numSigDigs,		/* Number of digits in the significand */
    int exponent)		/* Power of ten */
{
    double retval;		/* Value of the number */
    mp_int significandBig;	/* Significand expressed as a bignum */

    /*
     * With gcc on x86, the floating point rounding mode is double-extended.
     * This causes the result of double-precision calculations to be rounded
     * twice: once to the precision of double-extended and then again to the
     * precision of double. Double-rounding introduces gratuitous errors of 1
     * ulp, so we need to change rounding mode to 53-bits.
     */

#if defined(__GNUC__) && defined(__i386)
    fpu_control_t roundTo53Bits = 0x027f;
    fpu_control_t oldRoundingMode;
    _FPU_GETCW(oldRoundingMode);
    _FPU_SETCW(roundTo53Bits);
#endif

    /*
     * Test for the easy cases.
     */

    if (numSigDigs <= DBL_DIG) {
	if (exponent >= 0) {
	    if (exponent <= mmaxpow) {
		/* 
		 * The significand is an exact integer, and so is
		 * 10**exponent. The product will be correct to within 1/2 ulp
		 * without special handling.
		 */

		retval = (double)(Tcl_WideInt)significand * pow10[ exponent ];
		goto returnValue;
	    } else {
		int diff = DBL_DIG - numSigDigs;
		if (exponent-diff <= mmaxpow) {
		    /* 
		     * 10**exponent is not an exact integer, but
		     * 10**(exponent-diff) is exact, and so is
		     * significand*10**diff, so we can still compute the value
		     * with only one roundoff.
		     */

		    volatile double factor =
			    (double)(Tcl_WideInt)significand * pow10[diff];
		    retval = factor * pow10[exponent-diff];
		    goto returnValue;
		}
	    }
	} else {
	    if (exponent >= -mmaxpow) {
		/* 
		 * 10**-exponent is an exact integer, and so is the
		 * significand. Compute the result by one division, again with
		 * only one rounding.
		 */

		retval = (double)(Tcl_WideInt)significand / pow10[-exponent];
		goto returnValue;
	    }
	}
    }

    /*
     * All the easy cases have failed. Promote ths significand to bignum and
     * call MakeHighPrecisionDouble to do it the hard way.
     */

    TclBNInitBignumFromWideUInt(&significandBig, significand);
    retval = MakeHighPrecisionDouble(0, &significandBig, numSigDigs,
	    exponent);
    mp_clear(&significandBig);
    
    /*
     * Come here to return the computed value.
     */

  returnValue:
    if (signum) {
	retval = -retval;
    }

    /*
     * On gcc on x86, restore the floating point mode word.
     */

#if defined(__GNUC__) && defined(__i386)
    _FPU_SETCW(oldRoundingMode);
#endif

    return retval;
}

/*
 *----------------------------------------------------------------------
 *
 * MakeHighPrecisionDouble --
 *
 *	Makes the double precision number, signum*significand*10**exponent.
 *
 * Results:
 *	Returns the constructed number.
 *
 *	MakeHighPrecisionDouble is used when arbitrary-precision arithmetic is
 *	needed to ensure correct rounding. It begins by calculating a
 *	low-precision approximation to the desired number, and then refines
 *	the answer in high precision.
 *
 *----------------------------------------------------------------------
 */

static double
MakeHighPrecisionDouble(
    int signum,			/* 1=negative, 0=nonnegative */
    mp_int *significand,	/* Exact significand of the number */
    int numSigDigs,		/* Number of significant digits */
    int exponent)		/* Power of 10 by which to multiply */
{
    double retval;
    int machexp;		/* Machine exponent of a power of 10 */

    /*
     * With gcc on x86, the floating point rounding mode is double-extended.
     * This causes the result of double-precision calculations to be rounded
     * twice: once to the precision of double-extended and then again to the
     * precision of double. Double-rounding introduces gratuitous errors of 1
     * ulp, so we need to change rounding mode to 53-bits.
     */

#if defined(__GNUC__) && defined(__i386)
    fpu_control_t roundTo53Bits = 0x027f;
    fpu_control_t oldRoundingMode;
    _FPU_GETCW(oldRoundingMode);
    _FPU_SETCW(roundTo53Bits);
#endif

    /*
     * Quick checks for over/underflow.
     */

    if (numSigDigs+exponent-1 > maxDigits) {
	retval = HUGE_VAL;
	goto returnValue;
    }
    if (numSigDigs+exponent-1 < minDigits) {
	retval = 0;
	goto returnValue;
    }

    /* 
     * Develop a first approximation to the significand. It is tempting simply
     * to force bignum to double, but that will overflow on input numbers like
     * 1.[string repeat 0 1000]1; while this is a not terribly likely
     * scenario, we still have to deal with it. Use fraction and exponent
     * instead. Once we have the significand, multiply by 10**exponent. Test
     * for overflow. Convert back to a double, and test for underflow.
     */

    retval = BignumToBiasedFrExp(significand, &machexp);
    retval = Pow10TimesFrExp(exponent, retval, &machexp);
    if (machexp > DBL_MAX_EXP*log2FLT_RADIX) {
	retval = HUGE_VAL;
	goto returnValue;
    }
    retval = SafeLdExp(retval, machexp);
    if (retval < tiny) {
	retval = tiny;
    }

    /* 
     * Refine the result twice. (The second refinement should be necessary
     * only if the best approximation is a power of 2 minus 1/2 ulp).
     */

    retval = RefineApproximation(retval, significand, exponent);
    retval = RefineApproximation(retval, significand, exponent);

    /*
     * Come here to return the computed value.
     */

  returnValue:
    if (signum) {
	retval = -retval;
    }

    /*
     * On gcc on x86, restore the floating point mode word.
     */

#if defined(__GNUC__) && defined(__i386)
    _FPU_SETCW(oldRoundingMode);
#endif
    return retval;
}

/*
 *----------------------------------------------------------------------
 *
 * MakeNaN --
 *
 *	Makes a "Not a Number" given a set of bits to put in the tag bits
 *
 *	Note that a signalling NaN is never returned.
 *
 *----------------------------------------------------------------------
 */

#ifdef IEEE_FLOATING_POINT
static double
MakeNaN(
    int signum,			/* Sign bit (1=negative, 0=nonnegative */
    Tcl_WideUInt tags)	 	/* Tag bits to put in the NaN */
{
    union {
	Tcl_WideUInt iv;
	double dv;
    } theNaN;

    theNaN.iv = tags;
    theNaN.iv &= (((Tcl_WideUInt) 1) << 51) - 1;
    if (signum) {
	theNaN.iv |= ((Tcl_WideUInt) (0x8000 | NAN_START)) << 48;
    } else {
	theNaN.iv |= ((Tcl_WideUInt) NAN_START) << 48;
    }

    return theNaN.dv;
}
#endif

/*
 *----------------------------------------------------------------------
 *
 * RefineApproximation --
 *
 *	Given a poor approximation to a floating point number, returns a
 *	better one. (The better approximation is correct to within 1 ulp, and
 *	is entirely correct if the poor approximation is correct to 1 ulp.)
 *
 * Results:
 *	Returns the improved result.
 *
 *----------------------------------------------------------------------
 */

static double
RefineApproximation(
    double approxResult,	/* Approximate result of conversion */
    mp_int *exactSignificand,	/* Integer significand */
    int exponent)		/* Power of 10 to multiply by significand */
{
    int M2, M5;			/* Powers of 2 and of 5 needed to put the
				 * decimal and binary numbers over a common
				 * denominator. */
    double significand;		/* Sigificand of the binary number */
    int binExponent;		/* Exponent of the binary number */
    int msb;			/* Most significant bit position of an
				 * intermediate result */
    int nDigits;		/* Number of mp_digit's in an intermediate
				 * result */
    mp_int twoMv;		/* Approx binary value expressed as an exact
				 * integer scaled by the multiplier 2M */
    mp_int twoMd;		/* Exact decimal value expressed as an exact
				 * integer scaled by the multiplier 2M */
    int scale;			/* Scale factor for M */
    int multiplier;		/* Power of two to scale M */
    double num, den;		/* Numerator and denominator of the correction
				 * term */
    double quot;		/* Correction term */
    double minincr;		/* Lower bound on the absolute value of the
				 * correction term. */
    int i;

    /*
     * The first approximation is always low. If we find that it's HUGE_VAL,
     * we're done.
     */

    if (approxResult == HUGE_VAL) {
	return approxResult;
    }

    /*
     * Find a common denominator for the decimal and binary fractions. The
     * common denominator will be 2**M2 + 5**M5.
     */

    significand = frexp(approxResult, &binExponent);
    i = mantBits - binExponent;
    if (i < 0) {
	M2 = 0;
    } else {
	M2 = i;
    }
    if (exponent > 0) {
	M5 = 0;
    } else {
	M5 = -exponent;
	if ((M5-1) > M2) {
	    M2 = M5-1;
	}
    }

    /* 
     * The floating point number is significand*2**binExponent. Compute the
     * large integer significand*2**(binExponent+M2+1). The 2**-1 bit of the
     * significand (the most significant) corresponds to the
     * 2**(binExponent+M2 + 1) bit of 2*M2*v. Allocate enough digits to hold
     * that quantity, then convert the significand to a large integer, scaled
     * appropriately. Then multiply by the appropriate power of 5.
     */

    msb = binExponent + M2;	/* 1008 */
    nDigits = msb / DIGIT_BIT + 1;
    mp_init_size(&twoMv, nDigits);
    i = (msb % DIGIT_BIT + 1);
    twoMv.used = nDigits;
    significand *= SafeLdExp(1.0, i);
    while (--nDigits >= 0) {
	twoMv.dp[nDigits] = (mp_digit) significand;
	significand -= (mp_digit) significand;
	significand = SafeLdExp(significand, DIGIT_BIT);
    }
    for (i = 0; i <= 8; ++i) {
	if (M5 & (1 << i)) {
	    mp_mul(&twoMv, pow5+i, &twoMv);
	}
    }
    
    /* 
     * Collect the decimal significand as a high precision integer. The least
     * significant bit corresponds to bit M2+exponent+1 so it will need to be
     * shifted left by that many bits after being multiplied by
     * 5**(M5+exponent).
     */

    mp_init_copy(&twoMd, exactSignificand);
    for (i=0; i<=8; ++i) {
	if ((M5+exponent) & (1 << i)) {
	    mp_mul(&twoMd, pow5+i, &twoMd);
	}
    }
    mp_mul_2d(&twoMd, M2+exponent+1, &twoMd);
    mp_sub(&twoMd, &twoMv, &twoMd);

    /*
     * The result, 2Mv-2Md, needs to be divided by 2M to yield a correction
     * term. Because 2M may well overflow a double, we need to scale the
     * denominator by a factor of 2**binExponent-mantBits
     */

    scale = binExponent - mantBits - 1;

    mp_set(&twoMv, 1);
    for (i=0; i<=8; ++i) {
	if (M5 & (1 << i)) {
	    mp_mul(&twoMv, pow5+i, &twoMv);
	}
    }
    multiplier = M2 + scale + 1;
    if (multiplier > 0) {
	mp_mul_2d(&twoMv, multiplier, &twoMv);
    } else if (multiplier < 0) {
	mp_div_2d(&twoMv, -multiplier, &twoMv, NULL);
    }

    /*
     * If the result is less than unity, the error is less than 1/2 unit in
     * the last place, so there's no correction to make.
     */

    if (mp_cmp_mag(&twoMd, &twoMv) == MP_LT) {
        mp_clear(&twoMd);
        mp_clear(&twoMv);
	return approxResult;
    }

    /* 
     * Convert the numerator and denominator of the corrector term accurately
     * to floating point numbers.
     */

    num = TclBignumToDouble(&twoMd);
    den = TclBignumToDouble(&twoMv);

    quot = SafeLdExp(num/den, scale);
    minincr = SafeLdExp(1.0, binExponent-mantBits);

    if (quot<0. && quot>-minincr) {
	quot = -minincr;
    } else if (quot>0. && quot<minincr) {
	quot = minincr;
    }

    mp_clear(&twoMd);
    mp_clear(&twoMv);

    return approxResult + quot;
}

/*
 *----------------------------------------------------------------------
 *
 * TclDoubleDigits --
 *
 *	Converts a double to a string of digits.
 *
 * 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.
 *
 * Side effects:
 *	Stores the digits in the given buffer and sets 'signum' according to
 *	the sign of the number.
 *
 *----------------------------------------------------------------------

 */

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. */
{
    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;

    /*
     * Split the number into absolute value and signum.
     */

    v = AbsoluteValue(v, signum);

    /*
     * Handle zero specially.
     */

    if (v == 0.0) {
	*buffer++ = '0';
	*buffer++ = '\0';
	return 1;
    }

    /* 
     * Find a large integer r, and integer e, such that 
     *         v = r * FLT_RADIX**e
     * and r is as small as possible. Also determine whether the significand
     * is the smallest possible.
     */

    smallestSig = GetIntegerTimesPower(v, &r, &e);

    lowOK = highOK = (mp_iseven(&r));

    /*
     * 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.
     */

    if (e >= 0) {
	int bits = e * log2FLT_RADIX;

	if (!smallestSig) {
	    /*
	     * Normal case, m+ and m- are both FLT_RADIX**e
	     */

	    rfac2 = bits + 1;
	    sfac2 = 1;
	    mplusfac2 = bits;
	    mminusfac2 = bits;
	} else {
	    /*
	     * 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.
	     */

	    rfac2 = bits + log2FLT_RADIX + 1;
	    sfac2 = 1 + log2FLT_RADIX;
	    mplusfac2 = bits + log2FLT_RADIX;
	    mminusfac2 = bits;
	}
    } else {
	/*
	 * v has digits after the binary point
	 */

	if (e <= DBL_MIN_EXP-DBL_MANT_DIG || !smallestSig) {
	    /*
	     * Either f isn't the smallest significand or e is the smallest
	     * exponent. mplus and mminus will both be 1.
	     */

	    rfac2 = 1;
	    sfac2 = 1 - e * log2FLT_RADIX;
	    mplusfac2 = 0;
	    mminusfac2 = 0;
	} else {
	    /*
	     * 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.
	     */

	    rfac2 = 1 + log2FLT_RADIX;
	    sfac2 = 1 + log2FLT_RADIX * (1 - e);
	    mplusfac2 = FLT_RADIX;
	    mminusfac2 = 0;
	}
    }

    /*
     * Estimate the highest power of ten that will be needed to hold the
     * result.
     */

    k = (int) ceil(log(v) / log(10.));
    if (k >= 0) {
	sfac2 += k;
	sfac5 = k;
    } else {
	rfac2 -= k;
	mplusfac2 -= k;
	mminusfac2 -= k;
	rfac5 = -k;
    }

    /*
     * Scale r, s, mplus, mminus by the appropriate powers of 2 and 5.
     */

    mp_init_set(&mplus, 1);
    for (i=0 ; i<=8 ; ++i) {
	if (rfac5 & (1 << i)) {
	    mp_mul(&mplus, pow5+i, &mplus);
	}
    }
    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);
	}
    }
    mp_mul_2d(&s, sfac2, &s);

    /*
     * It is possible for k to be off by one because we used an inexact
     * logarithm.
     */

    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;
	}
    }

    /*
     * 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.
     */

    for (;;) {
	int tc1, tc2;

	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);
	}
	mp_add(&r, &mplus, &temp);
	tc2 = mp_cmp_mag(&temp, &s);
	if (highOK) {
	    tc2 = (tc2 >= 0);
	} else {
	    tc2= (tc2 > 0);
	}
	if (!tc1) {
	    if (!tc2) {
		*buffer++ = '0' + i;
	    } else {
		c = (char) (i + '1');
		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');
		}
	    }
	    break;
	}
    };
    *buffer++ = c;
    *buffer++ = '\0';

    /*
     * Free memory, and return.
     */

    mp_clear_multi(&r, &s, &mplus, &mminus, &temp, NULL);
    return k;
}

/*
 *----------------------------------------------------------------------
 *
 * AbsoluteValue --
 *
 *	Splits a 'double' into its absolute value and sign.
 *
 * Results:
 *	Returns the absolute value.
 *
 * Side effects:
 *	Stores the signum in '*signum'.
 *
 *----------------------------------------------------------------------
 */

static double
AbsoluteValue(
    double v,			/* Number to split */
    int *signum)		/* (Output) Sign of the number 1=-, 0=+ */
{
    /*
     * Take the absolute value of the number, and report the number's sign.
     * Take special steps to preserve signed zeroes in IEEE floating point.
     * (We can't use fpclassify, because that's a C9x feature and we still
     * have to build on C89 compilers.)
     */

#ifndef IEEE_FLOATING_POINT
    if (v >= 0.0) {
	*signum = 0;
    } else {
	*signum = 1;
	v = -v;
    }
#else
    union {
	Tcl_WideUInt iv;
	double dv;
    } bitwhack;
    bitwhack.dv = v;
    if (bitwhack.iv & ((Tcl_WideUInt) 1 << 63)) {
	*signum = 1;
	bitwhack.iv &= ~((Tcl_WideUInt) 1 << 63);
	v = bitwhack.dv;
    } else {
	*signum = 0;
    }
#endif
    return v;
}

/*
 *----------------------------------------------------------------------
 *
 * GetIntegerTimesPower --
 *
 *	Converts a floating point number to an exact integer times a power of
 *	the floating point radix.
 *
 * Results:
 *	Returns 1 if it converted the smallest significand, 0 otherwise.
 *
 * Side effects:
 *	Initializes the integer value (does not just assign it), and stores
 *	the exponent.
 *
 *----------------------------------------------------------------------
 */

static int
GetIntegerTimesPower(
    double v,			/* Value to convert */
    mp_int *rPtr,		/* (Output) Integer value */
    int *ePtr)			/* (Output) Power of FLT_RADIX by which r must
				 * be multiplied to yield v*/
{
    double a, f;
    int e, i, n;

    /*
     * Develop f and e such that v = f * FLT_RADIX**e, with
     * 1.0/FLT_RADIX <= f < 1.
     */

    f = frexp(v, &e);
#if FLT_RADIX > 2
    n = e % log2FLT_RADIX;
    if (n > 0) {
	n -= log2FLT_RADIX;
	e += 1;
	f *= ldexp(1.0, n);
    }
    e = (e - n) / log2FLT_RADIX;
#endif
    if (f == 1.0) {
	f = 1.0 / FLT_RADIX;
	e += 1;
    }

    /*
     * If the original number was denormalized, adjust e and f to be denormal
     * as well.
     */

    if (e < DBL_MIN_EXP) {
	n = mantBits + (e - DBL_MIN_EXP)*log2FLT_RADIX;
	f = ldexp(f, (e - DBL_MIN_EXP)*log2FLT_RADIX);
	e = DBL_MIN_EXP;
	n = (n + DIGIT_BIT - 1) / DIGIT_BIT;
    } else {
	n = mantDIGIT;
    }

    /*
     * Now extract the base-2**DIGIT_BIT digits of f into a multi-precision
     * integer r. Preserve the invariant v = r * 2**rfac2 * FLT_RADIX**e by
     * adjusting e.
     */

    a = f;
    n = mantDIGIT;
    mp_init_size(rPtr, n);
    rPtr->used = n;
    rPtr->sign = MP_ZPOS;
    i = (mantBits % DIGIT_BIT);
    if (i == 0) {
	i = DIGIT_BIT;
    }
    while (n > 0) {
	a *= ldexp(1.0, i);
	i = DIGIT_BIT;
	rPtr->dp[--n] = (mp_digit) a;
	a -= (mp_digit) a;
    }
    *ePtr = e - DBL_MANT_DIG;
    return (f == 1.0 / FLT_RADIX);
}

/*
 *----------------------------------------------------------------------
 *
 * TclInitDoubleConversion --
 *
 *	Initializes constants that are needed for conversions to and from
 *	'double'
 *
 * Results:
 *	None.
 *
 * Side effects:
 *	The log base 2 of the floating point radix, the number of bits in a
 *	double mantissa, and a table of the powers of five and ten are
 *	computed and stored.
 *
 *----------------------------------------------------------------------
 */

void
TclInitDoubleConversion(void)
{
    int i;
    int x;
    Tcl_WideUInt u;
    double d;

    /*
     * Initialize table of powers of 10 expressed as wide integers.
     */

    maxpow10_wide = (int)
	    floor(sizeof(Tcl_WideUInt) * CHAR_BIT * log(2.) / log(10.));
    pow10_wide = (Tcl_WideUInt *)
	    ckalloc((maxpow10_wide + 1) * sizeof(Tcl_WideUInt));
    u = 1;
    for (i = 0; i < maxpow10_wide; ++i) {
	pow10_wide[i] = u;
	u *= 10;
    }
    pow10_wide[i] = u;

    /*
     * Determine how many bits of precision a double has, and how many
     * decimal digits that represents.
     */

    if (frexp((double) FLT_RADIX, &log2FLT_RADIX) != 0.5) {
	Tcl_Panic("This code doesn't work on a decimal machine!");
    }
    --log2FLT_RADIX;
    mantBits = DBL_MANT_DIG * log2FLT_RADIX;
    d = 1.0;

    /* 
     * Initialize a table of powers of ten that can be exactly represented
     * in a double.
     */

    x = (int) (DBL_MANT_DIG * log((double) FLT_RADIX) / log(5.0));
    if (x < MAXPOW) {
	mmaxpow = x;
    } else {
	mmaxpow = MAXPOW;
    }
    for (i=0 ; i<=mmaxpow ; ++i) {
	pow10[i] = d;
	d *= 10.0;
    }

    /*
     * Initialize a table of large powers of five.
     */

    for (i=0; i<9; ++i) {
	mp_init(pow5 + i);
    }
    mp_set(pow5, 5);
    for (i=0; i<8; ++i) {
	mp_sqr(pow5+i, pow5+i+1);
    }

    /* 
     * Determine the number of decimal digits to the left and right of the
     * decimal point in the largest and smallest double, the smallest double
     * that differs from zero, and the number of mp_digits needed to represent 
     * the significand of a double.
     */

    tiny = SafeLdExp(1.0, DBL_MIN_EXP * log2FLT_RADIX - mantBits);
    maxDigits = (int) ((DBL_MAX_EXP * log((double) FLT_RADIX)
	    + 0.5 * log(10.)) / log(10.));
    minDigits = (int) floor((DBL_MIN_EXP - DBL_MANT_DIG)
	    * log((double) FLT_RADIX) / log(10.));
    mantDIGIT = (mantBits + DIGIT_BIT-1) / DIGIT_BIT;
    log10_DIGIT_MAX = (int) floor(DIGIT_BIT * log(2.) / log(10.));
}

/*
 *----------------------------------------------------------------------
 *
 * TclFinalizeDoubleConversion --
 *
 *	Cleans up this file on exit.
 *
 * Results:
 *	None
 *
 * Side effects:
 *	Memory allocated by TclInitDoubleConversion is freed.
 *
 *----------------------------------------------------------------------
 */

void
TclFinalizeDoubleConversion()
{
    int i;

    Tcl_Free((char*)pow10_wide);
    for (i=0; i<9; ++i) {
	mp_clear(pow5 + i);
    }
}

/*
 *----------------------------------------------------------------------
 *
 * TclInitBignumFromDouble --
 *
 *	Extracts the integer part of a double and converts it to an arbitrary
 *	precision integer.
 *
 * Results:
 *	None.
 *
 * Side effects:
 *	Initializes the bignum supplied, and stores the converted number * in
 *	it.
 *
 *----------------------------------------------------------------------
 */

int
TclInitBignumFromDouble(
    Tcl_Interp *interp,		/* For error message */
    double d,			/* Number to convert */
    mp_int *b)			/* Place to store the result */
{
    double fract;
    int expt;

    /*
     * Infinite values can't convert to bignum.
     */

    if (TclIsInfinite(d)) {
	if (interp != NULL) {
	    char *s = "integer value too large to represent";
	    Tcl_SetObjResult(interp, Tcl_NewStringObj(s, -1));
	    Tcl_SetErrorCode(interp, "ARITH", "IOVERFLOW", s, NULL);
	}
	return TCL_ERROR;
    }

    fract = frexp(d,&expt);
    if (expt <= 0) {
	mp_init(b);
	mp_zero(b);
    } else {
	Tcl_WideInt w = (Tcl_WideInt) ldexp(fract, mantBits);
	int shift = expt - mantBits;

	TclBNInitBignumFromWideInt(b, w);
	if (shift < 0) {
	    mp_div_2d(b, -shift, b, NULL);
	} else if (shift > 0) {
	    mp_mul_2d(b, shift, b);
	}
    }
    return TCL_OK;
}

/*
 *----------------------------------------------------------------------
 *
 * TclBignumToDouble --
 *
 *	Convert an arbitrary-precision integer to a native floating point
 *	number.
 *
 * Results:
 *	Returns the converted number. Sets errno to ERANGE if the number is
 *	too large to convert.
 *
 *----------------------------------------------------------------------
 */

double
TclBignumToDouble(
    mp_int *a)			/* Integer to convert. */
{
    mp_int b;
    int bits;
    int shift;
    int i;
    double r;

    /*
     * Determine how many bits we need, and extract that many from the input.
     * Round to nearest unit in the last place.
     */

    bits = mp_count_bits(a);
    if (bits > DBL_MAX_EXP*log2FLT_RADIX) {
	errno = ERANGE;
	if (a->sign == MP_ZPOS) {
	    return HUGE_VAL;
	} else {
	    return -HUGE_VAL;
	}
    }
    shift = mantBits + 1 - bits;
    mp_init(&b);
    if (shift > 0) {
	mp_mul_2d(a, shift, &b);
    } else if (shift < 0) {
	mp_div_2d(a, -shift, &b, NULL);
    } else {
	mp_copy(a, &b);
    }
    mp_add_d(&b, 1, &b);
    mp_div_2d(&b, 1, &b, NULL);

    /*
     * Accumulate the result, one mp_digit at a time.
     */

    r = 0.0;
    for (i=b.used-1 ; i>=0 ; --i) {
	r = ldexp(r, DIGIT_BIT) + b.dp[i];
    }
    mp_clear(&b);

    /*
     * Scale the result to the correct number of bits.
     */

    r = ldexp(r, bits - mantBits);

    /*
     * Return the result with the appropriate sign.
     */

    if (a->sign == MP_ZPOS) {
	return r;
    } else {
	return -r;
    }
}

double
TclCeil(
    mp_int *a)			/* Integer to convert. */
{
    double r = 0.0;
    mp_int b;

    mp_init(&b);
    if (mp_cmp_d(a, 0) == MP_LT) {
	mp_neg(a, &b);
	r = -TclFloor(&b);
    } else {
	int bits = mp_count_bits(a);

	if (bits > DBL_MAX_EXP*log2FLT_RADIX) {
	    r = HUGE_VAL;
	} else {
	    int i, exact = 1, shift = mantBits - bits;

	    if (shift > 0) {
		mp_mul_2d(a, shift, &b);
	    } else if (shift < 0) {
		mp_int d;
		mp_init(&d);
		mp_div_2d(a, -shift, &b, &d);
		exact = mp_iszero(&d);
		mp_clear(&d);
	    } else {
		mp_copy(a, &b);
	    }
	    if (!exact) {
		mp_add_d(&b, 1, &b);
	    }
	    for (i=b.used-1 ; i>=0 ; --i) {
		r = ldexp(r, DIGIT_BIT) + b.dp[i];
	    }
	    r = ldexp(r, bits - mantBits);
	}
    }
    mp_clear(&b);
    return r;
}

double
TclFloor(
    mp_int *a)			/* Integer to convert. */
{
    double r = 0.0;
    mp_int b;

    mp_init(&b);
    if (mp_cmp_d(a, 0) == MP_LT) {
	mp_neg(a, &b);
	r = -TclCeil(&b);
    } else {
	int bits = mp_count_bits(a);

	if (bits > DBL_MAX_EXP*log2FLT_RADIX) {
	    r = DBL_MAX;
	} else {
	    int i, shift = mantBits - bits;

	    if (shift > 0) {
		mp_mul_2d(a, shift, &b);
	    } else if (shift < 0) {
		mp_div_2d(a, -shift, &b, NULL);
	    } else {
		mp_copy(a, &b);
	    }
	    for (i=b.used-1 ; i>=0 ; --i) {
		r = ldexp(r, DIGIT_BIT) + b.dp[i];
	    }
	    r = ldexp(r, bits - mantBits);
	}
    }
    mp_clear(&b);
    return r;
}

/*
 *----------------------------------------------------------------------
 *
 * BignumToBiasedFrExp --
 *
 *	Convert an arbitrary-precision integer to a native floating point
 *	number in the range [0.5,1) times a power of two. NOTE: Intentionally
 *	converts to a number that's a few ulp too small, so that
 *	RefineApproximation will not overflow near the high end of the
 *	machine's arithmetic range.
 *
 * Results:
 *	Returns the converted number.
 *
 * Side effects:
 *	Stores the exponent of two in 'machexp'.
 *
 *----------------------------------------------------------------------
 */

static double
BignumToBiasedFrExp(
    mp_int *a,			/* Integer to convert */
    int *machexp)		/* Power of two */
{
    mp_int b;
    int bits;
    int shift;
    int i;
    double r;

    /*
     * Determine how many bits we need, and extract that many from the input.
     * Round to nearest unit in the last place.
     */

    bits = mp_count_bits(a);
    shift = mantBits - 2 - bits;
    mp_init(&b);
    if (shift > 0) {
	mp_mul_2d(a, shift, &b);
    } else if (shift < 0) {
	mp_div_2d(a, -shift, &b, NULL);
    } else {
	mp_copy(a, &b);
    }

    /*
     * Accumulate the result, one mp_digit at a time.
     */

    r = 0.0;
    for (i=b.used-1; i>=0; --i) {
	r = ldexp(r, DIGIT_BIT) + b.dp[i];
    }
    mp_clear(&b);

    /*
     * Return the result with the appropriate sign.
     */

    *machexp = bits - mantBits + 2;
    return ((a->sign == MP_ZPOS) ? r : -r);
}

/*
 *----------------------------------------------------------------------
 *
 * Pow10TimesFrExp --
 *
 *	Multiply a power of ten by a number expressed as fraction and
 *	exponent.
 *
 * Results:
 *	Returns the significand of the result.
 *
 * Side effects:
 *	Overwrites the 'machexp' parameter with the exponent of the result.
 *
 * Assumes that 'exponent' is such that 10**exponent would be a double, even
 * though 'fraction*10**(machexp+exponent)' might overflow.
 *
 *----------------------------------------------------------------------
 */

static double
Pow10TimesFrExp(
    int exponent,	 	/* Power of 10 to multiply by */
    double fraction,		/* Significand of multiplicand */
    int *machexp)		/* On input, exponent of multiplicand. On
				 * output, exponent of result. */
{
    int i, j;
    int expt = *machexp;
    double retval = fraction;

    if (exponent > 0) {
	/*
	 * Multiply by 10**exponent
	 */
	
	retval = frexp(retval * pow10[exponent&0xf], &j);
	expt += j;
	for (i=4; i<9; ++i) {
	    if (exponent & (1<<i)) {
		retval = frexp(retval * pow_10_2_n[i], &j);
		expt += j;
	    }
	}
    } else if (exponent < 0) {
	/*
	 * Divide by 10**-exponent
	 */

	retval = frexp(retval / pow10[(-exponent) & 0xf], &j);
	expt += j;
	for (i=4; i<9; ++i) {
	    if ((-exponent) & (1<<i)) {
		retval = frexp(retval / pow_10_2_n[i], &j);
		expt += j;
	    }
	}
    }

    *machexp = expt;
    return retval;
}

/*
 *----------------------------------------------------------------------
 *
 * SafeLdExp --
 *
 *	Do an 'ldexp' operation, but handle denormals gracefully.
 *
 * Results:
 *	Returns the appropriately scaled value.
 *
 *	On some platforms, 'ldexp' fails when presented with a number too
 *	small to represent as a normalized double. This routine does 'ldexp'
 *	in two steps for those numbers, to return correctly denormalized
 *	values.
 *
 *----------------------------------------------------------------------
 */

static double
SafeLdExp(
    double fract,
    int expt)
{
    int minexpt = DBL_MIN_EXP * log2FLT_RADIX;
    volatile double a, b, retval;

    if (expt < minexpt) {
	a = ldexp(fract, expt - mantBits - minexpt);
	b = ldexp(1.0, mantBits + minexpt);
	retval = a * b;
    } else {
	retval = ldexp(fract, expt);
    }
    return retval;
}

/*
 *----------------------------------------------------------------------
 *
 * TclFormatNaN --
 *
 *	Makes the string representation of a "Not a Number"
 *
 * Results:
 *	None.
 *
 * Side effects:
 *	Stores the string representation in the supplied buffer, which must be
 *	at least TCL_DOUBLE_SPACE characters.
 *
 *----------------------------------------------------------------------
 */

void
TclFormatNaN(
    double value,		/* The Not-a-Number to format. */
    char *buffer)		/* String representation. */
{
#ifndef IEEE_FLOATING_POINT
    strcpy(buffer, "NaN");
    return;
#else
    union {
	double dv;
	Tcl_WideUInt iv;
    } bitwhack;

    bitwhack.dv = value;
    if (bitwhack.iv & ((Tcl_WideUInt) 1 << 63)) {
	bitwhack.iv &= ~ ((Tcl_WideUInt) 1 << 63);
	*buffer++ = '-';
    }
    *buffer++ = 'N';
    *buffer++ = 'a';
    *buffer++ = 'N';
    bitwhack.iv &= (((Tcl_WideUInt) 1) << 51) - 1;
    if (bitwhack.iv != 0) {
	sprintf(buffer, "(%" TCL_LL_MODIFIER "x)", bitwhack.iv);
    } else {
	*buffer = '\0';
    }
#endif /* IEEE_FLOATING_POINT */
}

/*
 * Local Variables:
 * mode: c
 * c-basic-offset: 4
 * fill-column: 78
 * End:
 */