path: root/generic/tclStrToD.c

/*
 *----------------------------------------------------------------------
 *
 * tclStrToD.c --
 *
 *	This file contains a TclStrToD procedure that handles conversion of
 *	string to double, with correct rounding even where extended precision
 *	is needed to achieve that.  It also contains a TclDoubleDigits
 *	procedure that handles conversion of double to string (at least the
 *	significand), and several utility functions for interconverting
 *	'double' and the integer types.
 *
 * Copyright (c) 2005 by Kevin B. Kenny.  All rights reserved.
 *
 * See the file "license.terms" for information on usage and redistribution
 * of this file, and for a DISCLAIMER OF ALL WARRANTIES.
 *
 * RCS: @(#) $Id: tclStrToD.c,v 1.8 2005/08/24 15:15:45 kennykb Exp $
 *
 *----------------------------------------------------------------------
 */

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

/*
 * The stuff below is a bit of a hack so that this file can be used in
 * environments that include no UNIX, i.e. no errno: just arrange to use the
 * errno from tclExecute.c here.
 */

#ifdef TCL_GENERIC_ONLY
#   define NO_ERRNO_H
#endif

#ifdef NO_ERRNO_H
extern int errno;			/* Use errno from tclExecute.c. */
#   define ERANGE 34
#endif

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

/*
 * gcc on x86 needs access to rounding controls.  It is tempting to include
 * fpu_control.h, but that file exists only on Linux; it is missing on Cygwin
 * and MinGW. Most gcc-isms and ix86-isms are factored out here.
 */

#if defined(__GNUC__) && defined(__i386)
typedef unsigned int fpu_control_t __attribute__ ((__mode__ (__HI__)));
#   define _FPU_GETCW(cw)	__asm__ ("fnstcw %0" : "=m" (*&cw))
#   define _FPU_SETCW(cw)	__asm__ ("fldcw %0" : : "m" (*&cw))
#   define FPU_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

/*
 * There now follows a lot of static variables that are shared across all
 * threads but which are not guarded by mutexes. This is OK, because they are
 * only ever assigned _ONCE_ during Tcl's library initialization sequence.
 */

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
};
#define MAXPOW	22		/* Num of exactly representable powers of 10 */
static double pow10[MAXPOW+1];	/* The powers of ten that can be represented
				 * exactly as IEEE754 doubles. */
static int mmaxpow;		/* Largest power of ten that can be
				 * represented exactly in a 'double'. */
static int log2FLT_RADIX;	/* Logarithm of the floating point radix. */
static int mantBits;		/* Number of bits in a double's significand. */
static mp_int pow5[9];		/* Table of powers of 5**(2**n), up to
				 * 5**256. */
static double tiny;		/* The smallest representable double. */
static int maxDigits;		/* The maximum number of digits to the left of
				 * the decimal point of a double. */
static int minDigits;		/* The maximum number of digits to the right
				 * of the decimal point in a double. */
static int mantDIGIT;		/* Number of mp_digit's needed to hold the
				 * significand of a double. */

/* Static functions defined in this file */

static double	RefineResult(double approx, CONST char *start, int nDigits,
		    long exponent);
static double	ParseNaN(int signum, CONST char **end);
static double	SafeLdExp(double fraction, int exponent);

/*
 *----------------------------------------------------------------------
 *
 * TclStrToD --
 *
 *	Scans a double from a string.
 *
 * Results:
 *	Returns the scanned number. In the case of underflow, returns an
 *	appropriately signed zero; in the case of overflow, returns an
 *	appropriately signed HUGE_VAL.
 *
 * Side effects:
 *	Stores a pointer to the end of the scanned number in '*endPtr', if
 *	endPtr is not NULL. If '*endPtr' is equal to 's' on return from this
 *	function, it indicates that the input string could not be recognized
 *	as a number.  In the case of underflow or overflow, 'errno' is set to
 *	ERANGE.
 *
 *------------------------------------------------------------------------
 */

double
TclStrToD(CONST char *s,	/* String to scan. */
	  CONST char **endPtr)	/* Pointer to the end of the scanned number. */
{
    const char *p = s;
    const char *startOfSignificand = NULL;
				/* Start of the significand in the string. */
    int signum = 0;		/* Sign of the significand. */
    double exactSignificand = 0.0;
				/* Significand, represented exactly as a
				 * floating-point number. */
    int seenDigit = 0;		/* Flag == 1 if a digit has been seen. */
    int nSigDigs = 0;		/* Number of significant digits presented. */
    int nDigitsAfterDp = 0;	/* Number of digits after the decimal point. */
    int nTrailZero = 0;		/* Number of trailing zeros in the
				 * significand. */
    long exponent = 0;		/* Exponent. */
    int seenDp = 0;		/* Flag == 1 if decimal point has been seen. */
    char c;			/* One character extracted from the input. */
    volatile double v;		/* Scanned value; must be 'volatile double' on
				 * gc-ix86 to force correct rounding to IEEE
				 * double and not Intel double-extended. */
    int machexp;		/* Exponent of the machine rep of the scanned
				 * value. */
    int expt2;			/* Exponent for computing first approximation
				 * to the true value. */
    int i, j;

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

#ifdef ADJUST_FPU_CONTROL_WORD
    fpu_control_t roundTo53Bits = FPU_IEEE_ROUNDING;
    fpu_control_t oldRoundingMode;
    _FPU_GETCW(oldRoundingMode);
    _FPU_SETCW(roundTo53Bits);
#   define RestoreRoundingMode()	_FPU_SETCW(oldRoundingMode)
#else
#   define RestoreRoundingMode()	(void) 0 /* Do nothing */
#endif

    /*
     * Discard leading whitespace from input.
     */

    while (isspace(UCHAR(*p))) {
	++p;
    }

    /*
     * Determine the sign of the significand.
     */

    switch (*p) {
	case '-':
	    signum = 1;
	    /* FALLTHROUGH */
	case '+':
	    ++p;
    }

    /*
     * Discard leading zeroes from input.
     */

    while (*p == '0') {
	seenDigit = 1;
	++p;
    }

    /*
     * Scan digits from the significand. Simultaneously, keep track of the
     * number of digits after the decimal point. Maintain a pointer to the
     * start of the significand. Keep "exactSignificand" equal to the
     * conversion of the DBL_DIG most significant digits.
     */

    for (;;) {
	c = *p;
	if (c == '.' && !seenDp) {
	    seenDp = 1;
	    ++p;
	} else if (isdigit(UCHAR(c))) {
	    if (c == '0') {
		if (startOfSignificand != NULL) {
		    ++nTrailZero;
		}
	    } else {
		if (startOfSignificand == NULL) {
		    startOfSignificand = p;
		} else if (nTrailZero) {
		    if (nTrailZero + nSigDigs < DBL_DIG) {
			exactSignificand *= pow10[nTrailZero];
		    } else if (nSigDigs < DBL_DIG) {
			exactSignificand *= pow10[DBL_DIG - nSigDigs];
		    }
		    nSigDigs += nTrailZero;
		}
		if (nSigDigs < DBL_DIG) {
		    exactSignificand = 10. * exactSignificand + (c - '0');
		}
		++nSigDigs;
		nTrailZero = 0;
	    }
	    if (seenDp) {
		++nDigitsAfterDp;
	    }
	    seenDigit = 1;
	    ++p;
	} else {
	    break;
	}
    }

    /*
     * At this point, we've scanned the significand, and p points to the
     * character beyond it. "startOfSignificand" is the first non-zero
     * character in the significand. "nSigDigs" is the number of significant
     * digits of the significand, not including any trailing zeroes.
     * "exactSignificand" is a floating point number that represents, without
     * loss of precision, the first min(DBL_DIG,n) digits of the significand.
     * "nDigitsAfterDp" is the number of digits after the decimal point, again
     * excluding trailing zeroes.
     *
     * Now scan 'E' format
     */

    exponent = 0;
    if (seenDigit && (*p == 'e' || *p == 'E')) {
	const char* stringSave = p;
	++p;
	c = *p;
	if (isdigit(UCHAR(c)) || c == '+' || c == '-') {
	    errno = 0;
	    exponent = strtol(p, (char**)&p, 10);
	    if (errno == ERANGE) {
		if (exponent > 0) {
		    v = HUGE_VAL;
		} else {
		    v = 0.0;
		}
		*endPtr = p;
		goto returnValue;
	    }
	}
	if (p == stringSave+1) {
	    p = stringSave;
	    exponent = 0;
	}
    }
    exponent += nTrailZero - nDigitsAfterDp;

    /*
     * If we come here with no significant digits, we might still be looking
     * at Inf or NaN.  Go parse them.
     */

    if (!seenDigit) {
	/*
	 * Test for Inf or Infinity (in any case).
	 */

	if (c == 'I' || c == 'i') {
	    if ((p[1] == 'N' || p[1] == 'n')
		    && (p[2] == 'F' || p[2] == 'f')) {
		p += 3;
		if ((p[0] == 'I' || p[0] == 'i')
			&& (p[1] == 'N' || p[1] == 'n')
			&& (p[2] == 'I' || p[2] == 'i')
			&& (p[3] == 'T' || p[3] == 't')
			&& (p[4] == 'Y' || p[1] == 'y')) {
		    p += 5;
		}
		errno = ERANGE;
		v = HUGE_VAL;
		if (endPtr != NULL) {
		    *endPtr = p;
		}
		goto returnValue;
	    }

#ifdef IEEE_FLOATING_POINT
	    /*
	     * Only IEEE floating point supports NaN
	     */
	} else if ((c == 'N' || c == 'n')
		&& (sizeof(Tcl_WideUInt) == sizeof(double))) {
	    if ((p[1] == 'A' || p[1] == 'a')
		    && (p[2] == 'N' || p[2] == 'n')) {
		p += 3;

		if (endPtr != NULL) {
		    *endPtr = p;
		}

		/*
		 * Restore FPU mode word.
		 */

		RestoreRoundingMode();
		return ParseNaN(signum, endPtr);
	    }
#endif
	}

	goto error;
    }

    /*
     * We've successfully scanned; update the end-of-element pointer.
     */

    if (endPtr != NULL) {
	*endPtr = p;
    }

    /*
     * Test for zero.
     */

    if (nSigDigs == 0) {
	v = 0.0;
	goto returnValue;
    }

    /*
     * The easy cases are where we have an exact significand and the exponent
     * is small enough that we can compute the value with only one roundoff.
     * In addition to the cases where we can multiply or divide an
     * exact-integer significand by an exact-integer power of 10, there is
     * also David Gay's case where we can scale the significand by a power of
     * 10 (still keeping it exact) and then multiply by an exact power of 10.
     * The last case enables combinations like 83e25 that would otherwise
     * require high precision arithmetic.
     */

    if (nSigDigs <= DBL_DIG) {
	if (exponent >= 0) {
	    if (exponent <= mmaxpow) {
		v = exactSignificand * pow10[exponent];
		goto returnValue;
	    } else {
		int diff = DBL_DIG - nSigDigs;
		if (exponent - diff <= mmaxpow) {
		    volatile double factor = exactSignificand * pow10[diff];
		    v = factor * pow10[exponent - diff];
		    goto returnValue;
		}
	    }
	} else if (exponent >= -mmaxpow) {
	    v = exactSignificand / pow10[-exponent];
	    goto returnValue;
	}
    }

    /*
     * We don't have one of the easy cases, so we can't compute the scanned
     * number exactly, and have to do it in multiple precision.  Begin by
     * testing for obvious overflows and underflows.
     */

    if (nSigDigs + exponent - 1 > maxDigits) {
	v = HUGE_VAL;
	errno = ERANGE;
	goto returnValue;
    }
    if (nSigDigs + exponent - 1 < minDigits) {
	errno = ERANGE;
	v = 0.;
	goto returnValue;
    }

    /*
     * Nothing exceeds the boundaries of the tables, at least.  Compute an
     * approximate value for the number, with no possibility of overflow
     * because we manage the exponent separately.
     */

    if (nSigDigs > DBL_DIG) {
	expt2 = exponent + nSigDigs - DBL_DIG;
    } else {
	expt2 = exponent;
    }
    v = frexp(exactSignificand, &machexp);
    if (expt2 > 0) {
	v = frexp(v * pow10[expt2 & 0xf], &j);
	machexp += j;
	for (i=4 ; i<9 ; ++i) {
	    if (expt2 & (1 << i)) {
		v = frexp(v * pow_10_2_n[i], &j);
		machexp += j;
	    }
	}
    } else {
	v = frexp(v / pow10[(-expt2) & 0xf], &j);
	machexp += j;
	for (i=4 ; i<9 ; ++i) {
	    if ((-expt2) & (1 << i)) {
		v = frexp(v / pow_10_2_n[i], &j);
		machexp += j;
	    }
	}
    }

    /*
     * A first approximation is that the result will be v * 2 ** machexp. v is
     * greater than or equal to 0.5 and less than 1. If machexp >
     * DBL_MAX_EXP*log2(FLT_RADIX), there is an overflow. Constrain the result
     * to the smallest representible number to avoid premature underflow.
     */

    if (machexp > DBL_MAX_EXP * log2FLT_RADIX) {
	v = HUGE_VAL;
	errno = ERANGE;
	goto returnValue;
    }

    v = SafeLdExp(v, machexp);
    if (v < tiny) {
	v = tiny;
    }

    /*
     * We have a first approximation in v. Now we need to refine it.
     */

    v = RefineResult(v, startOfSignificand, nSigDigs, exponent);

    /*
     * In a very few cases, a second iteration is needed. e.g., 457e-102
     */

    v = RefineResult(v, startOfSignificand, nSigDigs, exponent);

    /*
     * Handle underflow.
     */

  returnValue:
    if (nSigDigs != 0 && v == 0.0) {
	errno = ERANGE;
    }

    /*
     * Return a number with correct sign.
     */

    if (signum) {
	v = -v;
    }

    /*
     * Restore FPU mode word and return.
     */

    RestoreRoundingMode();
    return v;

    /*
     * Come here on an invalid input.
     */

  error:
    if (endPtr != NULL) {
	*endPtr = s;
    }

    /*
     * Restore FPU mode word and return.
     */

    RestoreRoundingMode();
    return 0.0;
}

/*
 *----------------------------------------------------------------------
 *
 * RefineResult --
 *
 *	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
RefineResult(double approxResult, /* Approximate result of conversion. */
	     CONST char* sigStart,
				/* Pointer to start of significand in input
				 * string. */
	     int nSigDigs,	/* Number of significant digits. */
	     long exponent)	/* Power of ten to multiply by significand. */
{
    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;
    const char* p;

    /*
     * 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. 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(&twoMd);
    mp_zero(&twoMd);
    i = nSigDigs;
    for (p=sigStart ;; ++p) {
	char c = *p;
	if (isdigit(UCHAR(c))) {
	    mp_mul_d(&twoMd, (unsigned) 10, &twoMd);
	    mp_add_d(&twoMd, (unsigned) (c - '0'), &twoMd);
	    --i;
	    if (i == 0) {
		break;
	    }
	}
    }
    for (i=0 ; i<=8 ; ++i) {
	if ((M5+exponent) & (1 << i)) {
	    mp_mul(&twoMd, pow5+i, &twoMd);
	}
    }
    mp_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;
}

/*
 *----------------------------------------------------------------------
 *
 * ParseNaN --
 *
 *	Parses a "not a number" from an input string, and returns the double
 *	precision NaN corresponding to it.
 *
 * Side effects:
 *	Advances endPtr to follow any (hex) in the input string.
 *
 *	If the NaN is followed by a left paren, a string of spaes and
 *	hexadecimal digits, and a right paren, endPtr is advanced to follow
 *	it.
 *
 *	The string of hexadecimal digits is OR'ed into the resulting NaN, and
 *	the signum is set as well.  Note that a signalling NaN is never
 *	returned.
 *
 *----------------------------------------------------------------------
 */

static double
ParseNaN(int signum,		/* Flag == 1 if minus sign has been seen in
				 * front of NaN. */
	 CONST char** endPtr)	/* Pointer-to-pointer to char following "NaN"
				 * in the input string. */
{
    const char* p = *endPtr;
    char c;
    union {
	Tcl_WideUInt iv;
	double dv;
    } theNaN;

    /*
     * Scan off a hex number in parentheses.  Embedded blanks are ok.
     */

    theNaN.iv = 0;
    if (*p == '(') {
	++p;
	for (;;) {
	    c = *p++;
	    if (isspace(UCHAR(c))) {
		continue;
	    } else if (c == ')') {
		*endPtr = p;
		break;
	    } else if (isdigit(UCHAR(c))) {
		c -= '0';
	    } else if (c >= 'A' && c <= 'F') {
		c -= 'A' + 10;
	    } else if (c >= 'a' && c <= 'f') {
		c -= 'a' + 10;
	    } else {
		theNaN.iv = (((Tcl_WideUInt) NAN_START) << 48)
			| (((Tcl_WideUInt) signum) << 63);
		return theNaN.dv;
	    }
	    theNaN.iv = (theNaN.iv << 4) | c;
	}
    }

    /*
     * Mask the hex number down to the least significant 51 bits.
     */

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

    *endPtr = p;
    return theNaN.dv;
}

/*
 *----------------------------------------------------------------------
 *
 * TclDoubleDigits --
 *
 *	Converts a double to a string of digits.
 *
 * 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 * strPtr,	/* Buffer in which to store the result, must
				 * have at least 18 chars. */
		double v,	/* Number to convert. Must be finite, and not
				 * NaN. */
		int *signum)	/* Output: 1 if the number is negative.
				 * Should handle -0 correctly on the IEEE
				 * architecture. */
{
    double f;			/* Significand of v. */
    int e;			/* Power of FLT_RADIX that satisfies
				 * v = f * FLT_RADIX**e */
    int lowOK, highOK;
    mp_int r;			/* Scaled significand. */
    mp_int s;			/* Divisor such that v = r / s */
    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;
    double a;
    char c;
    int i, k, n;

    /*
     * Take the absolute value of the number, and report the number's sign.
     * Take special steps to preserve signed zeroes in IEEE floating point.
     * (We can't use fpclassify, because that's a C9x feature and we still
     * have to build on C89 compilers.)
     */

#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

    /*
     * Handle zero specially.
     */

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

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

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

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

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

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

    a = f;
    n = mantDIGIT;
    mp_init_size(&r, n);
    r.used = n;
    r.sign = MP_ZPOS;
    i = (mantBits % DIGIT_BIT);
    if (i == 0) {
	i = DIGIT_BIT;
    }
    while (n > 0) {
	a *= ldexp(1.0, i);
	i = DIGIT_BIT;
	r.dp[--n] = (mp_digit) a;
	a -= (mp_digit) a;
    }
    e -= DBL_MANT_DIG;

    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 (f != 1.0/FLT_RADIX) {
	    /*
	     * 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 || f != 1.0/FLT_RADIX) {
	    /*
	     * 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) {
		*strPtr++ = '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;
	}
    };
    *strPtr++ = c;
    *strPtr++ = '\0';

    /*
     * Free memory, and return.
     */

    mp_clear_multi(&r, &s, &mplus, &mminus, &temp, NULL);
    return k;
}

/*
 *----------------------------------------------------------------------
 *
 * 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;
    double d;

    if (frexp((double) FLT_RADIX, &log2FLT_RADIX) != 0.5) {
	Tcl_Panic("This code doesn't work on a decimal machine!");
    }
    --log2FLT_RADIX;
    mantBits = DBL_MANT_DIG * log2FLT_RADIX;
    d = 1.0;
    x = (int) (DBL_MANT_DIG * log((double) FLT_RADIX) / log(5.0));
    if (x < MAXPOW) {
	mmaxpow = x;
    } else {
	mmaxpow = MAXPOW;
    }
    for (i=0 ; i<=mmaxpow ; ++i) {
	pow10[i] = d;
	d *= 10.0;
    }
    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);
    }
    tiny = SafeLdExp(1.0, DBL_MIN_EXP * log2FLT_RADIX - mantBits);
    maxDigits = (int)
	    ((DBL_MAX_EXP * log((double) FLT_RADIX) + log(10.)/2) / log(10.));
    minDigits = (int)
	    floor((DBL_MIN_EXP-DBL_MANT_DIG)*log((double)FLT_RADIX)/log(10.));
    mantDIGIT = (mantBits + DIGIT_BIT - 1) / DIGIT_BIT;
}

/*
 *----------------------------------------------------------------------
 *
 * TclFinalizeDoubleConversion --
 *
 *	Cleans up this file on exit.
 *
 * Results:
 *	None
 *
 * Side effects:
 *	Memory allocated by TclInitDoubleConversion is freed.
 *
 *----------------------------------------------------------------------
 */

void
TclFinalizeDoubleConversion()
{
    int i;
    for (i=0 ; i<9 ; ++i) {
	mp_clear(pow5 + i);
    }
}

/*
 *----------------------------------------------------------------------
 *
 * 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;
	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;
    }
}

/*
 *----------------------------------------------------------------------
 *
 * 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 */
}