*
 * See the file "license.terms" for information on usage and redistribution of
 * this file, and for a DISCLAIMER OF ALL WARRANTIES.
 */

#include "tclInt.h"
#include "tclTomMath.h"
#include <float.h>
#include <math.h>

/*
 * Older MSVC has no copysign function, but it's available at least since
 * MSVC++ 12.0 (that is Visual Studio 2013).
 */

#ifdef _WIN32
#if (defined(_MSC_VER) && (_MSC_VER < 1800))
inline static double
copysign(double a, double b) {
    return _copysign(a, b);
#define copysign _copysign
}
#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.
 */

			    Tcl_WideUInt significand, int nSigDigs,
			    long exponent);
#ifdef IEEE_FLOATING_POINT
static double		MakeNaN(int signum, Tcl_WideUInt tag);
#endif
static double		RefineApproximation(double approx,
			    mp_int *exactSignificand, int exponent);
static void		MulPow5(mp_int *, unsigned, mp_int *);
static mp_err	MulPow5(mp_int *, unsigned, mp_int *) MP_WUR;
static int 		NormalizeRightward(Tcl_WideUInt *);
static int		RequiredPrecision(Tcl_WideUInt);
static void		DoubleToExpAndSig(double, Tcl_WideUInt *, int *,
			    int *);
static void		TakeAbsoluteValue(Double *, int *);
static char *		FormatInfAndNaN(Double *, int *, char **);
static char *		FormatZero(int *, char **);

    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 */
    int explicitOctal = 0;
    mp_err err = MP_OKAY;

#define ALL_BITS	((Tcl_WideUInt)-1)
#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 ((octalSignificandWide != 0)
				&& (((size_t)shift >=
					CHAR_BIT*sizeof(Tcl_WideUInt))
				|| (octalSignificandWide >
					((Tcl_WideUInt)-1 >> shift)))) {
			    octalSignificandOverflow = 1;
			    mp_init_u64(&octalSignificandBig,
			    err = mp_init_u64(&octalSignificandBig,
				    octalSignificandWide);
			}
		    }
		    if (!octalSignificandOverflow) {
			octalSignificandWide =
				(octalSignificandWide << shift) + (c - '0');
		    } else {
			if (err == MP_OKAY) {
			mp_mul_2d(&octalSignificandBig, shift,
				&octalSignificandBig);
			mp_add_d(&octalSignificandBig, (mp_digit)(c - '0'),
				&octalSignificandBig);
			    err = mp_mul_2d(&octalSignificandBig, shift,
				    &octalSignificandBig);
			}
			if (err == MP_OKAY) {
			    err = mp_add_d(&octalSignificandBig, (mp_digit)(c - '0'),
				    &octalSignificandBig);
			}
		    }
		    if (err != MP_OKAY) {
			return TCL_ERROR;
		    }
		}
		if (numSigDigs != 0) {
		    numSigDigs += numTrailZeros+1;
		} else {
		    numSigDigs = 1;
		}

		     * large shifts first.
		     */

		    if (significandWide != 0 &&
			    ((size_t)shift >= CHAR_BIT*sizeof(Tcl_WideUInt) ||
			    significandWide > ((Tcl_WideUInt)-1 >> shift))) {
			significandOverflow = 1;
			mp_init_u64(&significandBig,
			err = mp_init_u64(&significandBig,
				significandWide);
		    }
		}
		if (!significandOverflow) {
		    significandWide = (significandWide << shift) + d;
		} else {
		    mp_mul_2d(&significandBig, shift, &significandBig);
		    mp_add_d(&significandBig, (mp_digit) d, &significandBig);
		}
		} else if (err == MP_OKAY) {
		    err = mp_mul_2d(&significandBig, shift, &significandBig);
		    if (err == MP_OKAY) {
			err = mp_add_d(&significandBig, (mp_digit) d, &significandBig);
		    }
		}
	    }
	    if (err != MP_OKAY) {
		return TCL_ERROR;
	    }
	    numTrailZeros = 0;
	    state = HEXADECIMAL;
	    break;

	case BINARY:
	    acceptState = state;

		     * large shifts first.
		     */

		    if (significandWide != 0 &&
			    ((size_t)shift >= CHAR_BIT*sizeof(Tcl_WideUInt) ||
			    significandWide > ((Tcl_WideUInt)-1 >> shift))) {
			significandOverflow = 1;
			mp_init_u64(&significandBig,
			err = mp_init_u64(&significandBig,
				significandWide);
		    }
		}
		if (!significandOverflow) {
		    significandWide = (significandWide << shift) + 1;
		} else {
		    mp_mul_2d(&significandBig, shift, &significandBig);
		    mp_add_d(&significandBig, (mp_digit) 1, &significandBig);
		}
		} else if (err == MP_OKAY) {
		    err = mp_mul_2d(&significandBig, shift, &significandBig);
		    if (err == MP_OKAY) {
			err = mp_add_d(&significandBig, (mp_digit) 1, &significandBig);
		    }
		}
	    }
	    if (err != MP_OKAY) {
		return TCL_ERROR;
	    }
	    numTrailZeros = 0;
	    state = BINARY;
	    break;

	case ZERO_D:
	    if (c == '0') {

		    acceptState, bytes);
	case BINARY:
	    shift = numTrailZeros;
	    if (!significandOverflow && significandWide != 0 &&
		    ((size_t)shift >= CHAR_BIT*sizeof(Tcl_WideUInt) ||
		    significandWide > (MOST_BITS + signum) >> shift)) {
		significandOverflow = 1;
		mp_init_u64(&significandBig, significandWide);
		err = mp_init_u64(&significandBig, significandWide);
	    }
	    if (shift) {
		if (!significandOverflow) {
		    significandWide <<= shift;
		} else {
		    mp_mul_2d(&significandBig, shift, &significandBig);
		} else if (err == MP_OKAY) {
		    err = mp_mul_2d(&significandBig, shift, &significandBig);
		}
	    }
	    if (err != MP_OKAY) {
		return TCL_ERROR;
	    }
	    goto returnInteger;

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

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

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

	    shift = 3 * numTrailZeros;
	    if (!octalSignificandOverflow && octalSignificandWide != 0 &&
		    ((size_t)shift >= CHAR_BIT*sizeof(Tcl_WideUInt) ||
		    octalSignificandWide > (MOST_BITS + signum) >> shift)) {
		octalSignificandOverflow = 1;
		mp_init_u64(&octalSignificandBig,
		err = mp_init_u64(&octalSignificandBig,
			octalSignificandWide);
	    }
	    if (shift) {
		if (!octalSignificandOverflow) {
		    octalSignificandWide <<= shift;
		} else {
		    mp_mul_2d(&octalSignificandBig, shift,
		} else if (err == MP_OKAY) {
		    err = mp_mul_2d(&octalSignificandBig, shift,
			    &octalSignificandBig);
		}
	    }
	    if (!octalSignificandOverflow) {
		if (octalSignificandWide > (MOST_BITS + signum)) {
		    mp_init_u64(&octalSignificandBig,
		if ((err == MP_OKAY) && (octalSignificandWide > (MOST_BITS + signum))) {
		    err = mp_init_u64(&octalSignificandBig,
			    octalSignificandWide);
		    octalSignificandOverflow = 1;
		} else {
		    objPtr->typePtr = &tclIntType;
		    if (signum) {
			objPtr->internalRep.wideValue =
				- (Tcl_WideInt) octalSignificandWide;
		    } else {
			objPtr->internalRep.wideValue =
				(Tcl_WideInt) octalSignificandWide;
		    }
		}
	    }
	    if (octalSignificandOverflow) {
	    if ((err == MP_OKAY) && octalSignificandOverflow) {
		if (signum) {
		    (void)mp_neg(&octalSignificandBig, &octalSignificandBig);
		    err = mp_neg(&octalSignificandBig, &octalSignificandBig);
		}
		TclSetBignumIntRep(objPtr, &octalSignificandBig);
	    }
	    if (err != MP_OKAY) {
		return TCL_ERROR;
	    }
	    break;

	case ZERO:
	case DECIMAL:
	    significandOverflow = AccumulateDecimalDigit(0, numTrailZeros-1,
		    &significandWide, &significandBig, significandOverflow);
	    if (!significandOverflow && (significandWide > MOST_BITS+signum)){
	    if ((err == MP_OKAY) && !significandOverflow && (significandWide > MOST_BITS+signum)) {
		significandOverflow = 1;
		mp_init_u64(&significandBig, significandWide);
		err = mp_init_u64(&significandBig, significandWide);
	    }
	returnInteger:
	    if (!significandOverflow) {
		if (significandWide > MOST_BITS+signum) {
		    mp_init_u64(&significandBig,
		if ((err == MP_OKAY) && (significandWide > MOST_BITS+signum)) {
		    err = mp_init_u64(&significandBig,
			    significandWide);
		    significandOverflow = 1;
		} else {
		    objPtr->typePtr = &tclIntType;
		    if (signum) {
			objPtr->internalRep.wideValue =
				- (Tcl_WideInt) significandWide;
		    } else {
			objPtr->internalRep.wideValue =
				(Tcl_WideInt) significandWide;
		    }
		}
	    }
	    if (significandOverflow) {
	    if ((err == MP_OKAY) && significandOverflow) {
		if (signum) {
		    (void)mp_neg(&significandBig, &significandBig);
		    err = mp_neg(&significandBig, &significandBig);
		}
		TclSetBignumIntRep(objPtr, &significandBig);
	    }
	    if (err != MP_OKAY) {
		return TCL_ERROR;
	    }
	    break;

	case FRACTION:
	case EXPONENT:

	    /*

	    }
	    objPtr->typePtr = &tclDoubleType;
	    break;

#ifdef IEEE_FLOATING_POINT
	case sNAN:
	case sNANFINISH:
	    objPtr->internalRep.doubleValue = MakeNaN(signum,significandWide);
	    objPtr->internalRep.doubleValue = MakeNaN(signum, significandWide);
	    objPtr->typePtr = &tclDoubleType;
	    break;
#endif
	case INITIAL:
	    /* This case only to silence compiler warning. */
	    Tcl_Panic("TclParseNumber: state INITIAL can't happen here");
	}

	} else if (numZeros >= maxpow10_wide
		|| w > ((Tcl_WideUInt)-1-digit)/pow10_wide[numZeros+1]) {
	    /*
	     * Wide multiplication will overflow.  Expand the number to a
	     * bignum and fall through into the bignum case.
	     */

	    mp_init_u64(bignumRepPtr, w);
	    if (mp_init_u64(bignumRepPtr, w) != MP_OKAY) {
		return 0;
	    }
	} else {
	    /*
	     * Wide multiplication.
	     */

	    *wideRepPtr = w * pow10_wide[numZeros+1] + digit;
	    return 0;
	}
    }

    /*
     * Bignum multiplication.
     */

    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);
	if ((mp_mul_d(bignumRepPtr, (mp_digit) pow10_wide[numZeros+1],
		bignumRepPtr) != MP_OKAY)
		|| (mp_add_d(bignumRepPtr, (mp_digit) digit, bignumRepPtr) != MP_OKAY))
	return 0;
    } else {
	mp_err err;
	/*
	 * 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 MP_DIGIT_BIT is at least 27. The
	 * first multiplication, by up to 10**7, is done with a one-DIGIT
	 * multiply (this presumes that MP_DIGIT_BIT >= 24).
	 */

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

    return 1;
}

/*
 *----------------------------------------------------------------------

    }

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

    mp_init_u64(&significandBig, significand);
    if (mp_init_u64(&significandBig, significand) != MP_OKAY) {
	return 0.0;
    }
    retval = MakeHighPrecisionDouble(0, &significandBig, numSigDigs,
	    exponent);
    mp_clear(&significandBig);

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

MakeHighPrecisionDouble(
    int signum,			/* 1=negative, 0=nonnegative */
    mp_int *significand,	/* Exact significand of the number */
    int numSigDigs,		/* Number of significant digits */
    long exponent)		/* Power of 10 by which to multiply */
{
    double retval;
    int machexp;		/* Machine exponent of a power of 10. */
    int machexp = 0;		/* 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.

    int rteExponent;		/* Exponent of the round-to-even result */
    int shift;			/* Shift count for converting numerator
				 * and denominator of corrector to floating
				 * point */
    Tcl_WideInt rteSigWide;	/* Wide integer version of the significand
				 * for testing evenness */
    int i;
    mp_err err = MP_OKAY;

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

    if (approxResult == HUGE_VAL) {

     * Compute twoMv as 2*M*v, where v is the approximate value.
     * This is done by bit-whacking to calculate 2**(M2+1)*significand,
     * and then multiplying by 5**M5.
     */

    msb = binExponent + M2;	/* 1008 */
    nDigits = msb / MP_DIGIT_BIT + 1;
    mp_init_size(&twoMv, nDigits);
    if (mp_init_size(&twoMv, nDigits) != MP_OKAY) {
	return approxResult;
    }
    i = (msb % MP_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, MP_DIGIT_BIT);
    }
    for (i = 0; i <= 8; ++i) {
	if (M5 & (1 << i)) {
	    mp_mul(&twoMv, pow5+i, &twoMv);
	if (M5 & (1 << i) && (mp_mul(&twoMv, pow5+i, &twoMv) != MP_OKAY)) {
	    mp_clear(&twoMv);
	    return approxResult;
	}
    }

    /*
     * Compute twoMd as 2*M*d, where d is the exact value.
     * This is done by multiplying by 5**(M5+exponent) and then multiplying
     * by 2**(M5+exponent+1), which is, of couse, a left shift.
     */

    mp_init_copy(&twoMd, exactSignificand);
    for (i=0; i<=8; ++i) {
    if (mp_init_copy(&twoMd, exactSignificand) != MP_OKAY) {
	mp_clear(&twoMv);
	return approxResult;
    }
    for (i = 0; (i <= 8); ++i) {
	if ((M5 + exponent) & (1 << i)) {
	    mp_mul(&twoMd, pow5+i, &twoMd);
	    err = mp_mul(&twoMd, pow5+i, &twoMd);
	}
    }
    if (err == MP_OKAY) {
    mp_mul_2d(&twoMd, M2+exponent+1, &twoMd);
	err = mp_mul_2d(&twoMd, M2+exponent+1, &twoMd);
    }

    /*
     * Now let twoMd = twoMd - twoMv, the difference between the exact and
     * approximate values.
     */

    if (err == MP_OKAY) {
    mp_sub(&twoMd, &twoMv, &twoMd);
	err = 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. Place that factor
     * times 1/2 ULP into twoMd.
     */

    scale = binExponent - mantBits - 1;
    mp_set_u64(&twoMv, 1);
    for (i=0; i<=8; ++i) {
    for (i = 0; (i <= 8) && (err == MP_OKAY); ++i) {
	if (M5 & (1 << i)) {
	    mp_mul(&twoMv, pow5+i, &twoMv);
	    err = mp_mul(&twoMv, pow5+i, &twoMv);
	}
    }
    multiplier = M2 + scale + 1;
    if (err != MP_OKAY) {
	mp_clear(&twoMd);
	mp_clear(&twoMv);
	return approxResult;
    if (multiplier > 0) {
	mp_mul_2d(&twoMv, multiplier, &twoMv);
    } else if (multiplier > 0) {
	err = mp_mul_2d(&twoMv, multiplier, &twoMv);
    } else if (multiplier < 0) {
	mp_div_2d(&twoMv, -multiplier, &twoMv, NULL);
	err = mp_div_2d(&twoMv, -multiplier, &twoMv, NULL);
    }
    if (err != MP_OKAY) {
	mp_clear(&twoMd);
	mp_clear(&twoMv);
	return approxResult;
    }

    /*
     * Will the eventual correction term be less than, equal to, or
     * greater than 1/2 ULP?
     */

    /*
     * Reduce the numerator and denominator of the corrector term so that
     * they will fit in the floating point precision.
     */
    shift = mp_count_bits(&twoMv) - FP_PRECISION - 1;
    if (shift > 0) {
	mp_div_2d(&twoMv, shift, &twoMv, NULL);
	mp_div_2d(&twoMd, shift, &twoMd, NULL);
	err = mp_div_2d(&twoMv, shift, &twoMv, NULL);
	if (err == MP_OKAY) {
	    err = mp_div_2d(&twoMd, shift, &twoMd, NULL);
	}
    }
    if (err != MP_OKAY) {
	mp_clear(&twoMd);
	mp_clear(&twoMv);
	return approxResult;
    }

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

 *
 * Side effects:
 *	Stores base*5**n in result.
 *
 *----------------------------------------------------------------------
 */

static inline void
static inline mp_err
MulPow5(
    mp_int *base, 		/* Number to multiply. */
    unsigned n,			/* Power of 5 to multiply by. */
    mp_int *result)		/* Place to store the result. */
{
    mp_int *p = base;
    int n13 = n / 13;
    int r = n % 13;
    mp_err err = MP_OKAY;

    if (r != 0) {
	mp_mul_d(p, dpow5[r], result);
	err = mp_mul_d(p, dpow5[r], result);
	p = result;
    }
    r = 0;
    while (n13 != 0) {
    while ((err == MP_OKAY) && (n13 != 0)) {
	if (n13 & 1) {
	    mp_mul(p, pow5_13+r, result);
	    err = mp_mul(p, pow5_13+r, result);
	    p = result;
	}
	n13 >>= 1;
	++r;
    }
    if (p != result) {
	mp_copy(p, result);
    if ((err == MP_OKAY) && (p != result)) {
	err = mp_copy(p, result);
    }
    return err;
}

/*
 *----------------------------------------------------------------------
 *
 * NormalizeRightward --
 *

    /*
     * Compare B and S-m - which is the same as comparing B+m and S - which we
     * do by computing b+m and doing a bitwhack compare against
     * 2**(MP_DIGIT_BIT*sd)
     */

    mp_add(b, m, temp);
    if (temp->used <= sd) {	/* Too few digits to be > s */
    if ((mp_add(b, m, temp) != MP_OKAY) || (temp->used <= sd)) {	/* Too few digits to be > s */
	return 0;
    }
    if (temp->used > sd+1 || temp->dp[sd] > 1) {
				/* >= 2s */
	return 1;
    }
    for (i = sd-1; i >= 0; --i) {

				 * converted. */
    mp_int mplus, mminus;	/* Bounds for roundoff. */
    mp_digit digit;		/* Current output digit. */
    char *s = retval;		/* Cursor in the output buffer. */
    int i;			/* Index in the output buffer. */
    mp_int temp;
    int r1;
    mp_err err = MP_OKAY;

    /*
     * b = bw * 2**b2 * 5**b5
     * mminus = 5**m5
     */

    mp_init_u64(&b, bw);
    mp_init_set(&mminus, 1);
    MulPow5(&b, b5, &b);
    mp_mul_2d(&b, b2, &b);
    if ((retval == NULL) || (mp_init_u64(&b, bw) != MP_OKAY)) {
	return NULL;
    }
    if (mp_init_set(&mminus, 1) != MP_OKAY) {
	mp_clear(&b);
	return NULL;
    }
    err = MulPow5(&b, b5, &b);
    if (err == MP_OKAY) {
	err = mp_mul_2d(&b, b2, &b);
    }

    /*
     * Adjust if the logarithm was guessed wrong.
     */

    if (b.used <= sd) {
	mp_mul_d(&b, 10, &b);
    if ((err == MP_OKAY) && (b.used <= sd)) {
	err = mp_mul_d(&b, 10, &b);
	++m2plus; ++m2minus; ++m5;
	ilim = ilim1;
	--k;
    }

    /*
     * mminus = 5**m5 * 2**m2minus
     * mplus = 5**m5 * 2**m2plus
     */

    if (err == MP_OKAY) {
    mp_mul_2d(&mminus, m2minus, &mminus);
    MulPow5(&mminus, m5, &mminus);
    if (m2plus > m2minus) {
	mp_init_copy(&mplus, &mminus);
	mp_mul_2d(&mplus, m2plus-m2minus, &mplus);
    }
    mp_init(&temp);
	err = mp_mul_2d(&mminus, m2minus, &mminus);
    }
    if (err == MP_OKAY) {
	err = MulPow5(&mminus, m5, &mminus);
    }
    if ((err == MP_OKAY) && (m2plus > m2minus)) {
	err = mp_init_copy(&mplus, &mminus);
	if (err == MP_OKAY) {
	    err = mp_mul_2d(&mplus, m2plus-m2minus, &mplus);
	}
    }
    if (err == MP_OKAY) {
	err = mp_init(&temp);
    }

    /*
     * Loop through the digits. Do division and mod by s == 2**(sd*MP_DIGIT_BIT)
     * by mp_digit extraction.
     */

    i = 0;

	    break;
	}

	/*
	 * Advance to the next digit.
	 */

	if (err == MP_OKAY) {
	mp_mul_d(&b, 10, &b);
	mp_mul_d(&mminus, 10, &mminus);
	if (m2plus > m2minus) {
	    mp_mul_2d(&mminus, m2plus-m2minus, &mplus);
	    err = mp_mul_d(&b, 10, &b);
	}
	if (err == MP_OKAY) {
	    err = mp_mul_d(&mminus, 10, &mminus);
	}
	if ((err == MP_OKAY) && (m2plus > m2minus)) {
	    err = mp_mul_2d(&mminus, m2plus-m2minus, &mplus);
	}
	++i;
    }

    /*
     * Endgame - store the location of the decimal point and the end of the
     * string.
     */

    if (m2plus > m2minus) {
	mp_clear(&mplus);
    }
    mp_clear_multi(&b, &mminus, &temp, NULL);
    *s = '\0';
    *decpt = k;
    if (endPtr) {
	*endPtr = s;
    }
    return retval;
    return (err == MP_OKAY) ? retval : NULL;
}

/*
 *----------------------------------------------------------------------
 *
 * StrictBignumConversionPowD --
 *

    char *retval = (char *)ckalloc(len + 1);
				/* Output buffer. */
    mp_int b;			/* Numerator of the fraction being
				 * converted. */
    mp_digit digit;		/* Current output digit. */
    char *s = retval;		/* Cursor in the output buffer. */
    int i;			/* Index in the output buffer. */
    mp_err err;
    (void)dPtr;

    /*
     * b = bw * 2**b2 * 5**b5
     */

    mp_init_u64(&b, bw);
    MulPow5(&b, b5, &b);
    mp_mul_2d(&b, b2, &b);
    if (mp_init_u64(&b, bw) != MP_OKAY) {
	return NULL;
    }
    err = MulPow5(&b, b5, &b);
    if (err == MP_OKAY) {
	err = mp_mul_2d(&b, b2, &b);
    }

    /*
     * Adjust if the logarithm was guessed wrong.
     */

    if (b.used <= sd) {
	mp_mul_d(&b, 10, &b);
    if ((err == MP_OKAY) && (b.used <= sd)) {
	err = mp_mul_d(&b, 10, &b);
	ilim = ilim1;
	--k;
    }

    /*
     * Loop through the digits. Do division and mod by s == 2**(sd*MP_DIGIT_BIT)
     * by mp_digit extraction.
     */

    i = 1;
    for (;;) {
    while (err == MP_OKAY) {
	if (b.used <= sd) {
	    digit = 0;
	} else {
	    digit = b.dp[sd];
	    if (b.used > sd+1 || digit >= 10) {
		Tcl_Panic("wrong digit!");
	    }

	    break;
	}

	/*
	 * Advance to the next digit.
	 */

	mp_mul_d(&b, 10, &b);
	err = mp_mul_d(&b, 10, &b);
	++i;
    }

    /*
     * Endgame - store the location of the decimal point and the end of the
     * string.
     */

    int r;
    mp_int temp;

    /*
     * Compare b and S-m: this is the same as comparing B+m and S.
     */

    mp_init(&temp);
    mp_add(b, m, &temp);
    if ((mp_init(&temp) != MP_OKAY) || (mp_add(b, m, &temp) != MP_OKAY)) {
	return 0;
    }
    r = mp_cmp_mag(&temp, S);
    mp_clear(&temp);
    switch(r) {
    case MP_LT:
	return 0;
    case MP_EQ:
	return isodd;

    mp_int mplus;		/* 1/2 ulp above the result. */
    mp_int S;			/* Denominator of the result. */
    mp_int dig;			/* Current digit of the result. */
    int digit;			/* Current digit of the result. */
    int minit = 1;		/* Fudge factor for when we misguess k. */
    int i;
    int r1;
    mp_err err;

    /*
     * b = bw * 2**b2 * 5**b5
     * S = 2**s2 * 5*s5
     */

    mp_init_u64(&b, bw);
    mp_mul_2d(&b, b2, &b);
    mp_init_set(&S, 1);
    MulPow5(&S, s5, &S); mp_mul_2d(&S, s2, &S);
    if ((retval == NULL) || (mp_init_u64(&b, bw) != MP_OKAY)) {
	return NULL;
    }
    err = mp_mul_2d(&b, b2, &b);
    if (err == MP_OKAY) {
	err = mp_init_set(&S, 1);
    }
    if (err == MP_OKAY) {
	err = MulPow5(&S, s5, &S);
    }
    if (err == MP_OKAY) {
	err = mp_mul_2d(&S, s2, &S);
    }

    /*
     * Handle the case where we guess the position of the decimal point wrong.
     */

    if (mp_cmp_mag(&b, &S) == MP_LT) {
	mp_mul_d(&b, 10, &b);
    if ((err == MP_OKAY) && (mp_cmp_mag(&b, &S) == MP_LT)) {
	err = mp_mul_d(&b, 10, &b);
	minit = 10;
	ilim =ilim1;
	--k;
    }

    /*
     * mminus = 2**m2minus * 5**m5
     */

    if (err == MP_OKAY) {
    mp_init_set(&mminus, minit);
    mp_mul_2d(&mminus, m2minus, &mminus);
    if (m2plus > m2minus) {
	mp_init_copy(&mplus, &mminus);
	mp_mul_2d(&mplus, m2plus-m2minus, &mplus);
	err = mp_init_set(&mminus, minit);
    }
    if (err == MP_OKAY) {
	err = mp_mul_2d(&mminus, m2minus, &mminus);
    }
    if ((err == MP_OKAY) && (m2plus > m2minus)) {
	err = mp_init_copy(&mplus, &mminus);
	if (err == MP_OKAY) {
	    err = mp_mul_2d(&mplus, m2plus-m2minus, &mplus);
	}
    }

    /*
     * Loop through the digits.
     */

    if (err == MP_OKAY) {
    mp_init(&dig);
	err = mp_init(&dig);
    }
    i = 1;
    for (;;) {
	mp_div(&b, &S, &dig, &b);
    while (err == MP_OKAY) {
	err = mp_div(&b, &S, &dig, &b);
	if (dig.used > 1 || dig.dp[0] >= 10) {
	    Tcl_Panic("wrong digit!");
	}
	digit = dig.dp[0];

	/*
	 * Does the current digit leave us with a remainder small enough to
	 * round to it?
	 */

	r1 = mp_cmp_mag(&b, (m2plus > m2minus)? &mplus : &mminus);
	if (r1 == MP_LT || (r1 == MP_EQ && (dPtr->w.word1 & 1) == 0)) {
	    mp_mul_2d(&b, 1, &b);
	    err = mp_mul_2d(&b, 1, &b);
	    if (ShouldBankerRoundUp(&b, &S, digit&1)) {
		++digit;
		if (digit == 10) {
		    *s++ = '9';
		    s = BumpUp(s, retval, &k);
		    break;
		}

	}

	/*
	 * Have we converted all the requested digits?
	 */

	*s++ = '0' + digit;
	if (i == ilim) {
	    mp_mul_2d(&b, 1, &b);
	if ((err == MP_OKAY) && (i == ilim)) {
	    err = mp_mul_2d(&b, 1, &b);
	    if (ShouldBankerRoundUp(&b, &S, digit&1)) {
		s = BumpUp(s, retval, &k);
	    }
	    break;
	}

	/*
	 * Advance to the next digit.
	 */

	if (s5 > 0) {
	if ((err == MP_OKAY) && (s5 > 0)) {
	    /*
	     * Can possibly shorten the denominator.
	     */

	    mp_mul_2d(&b, 1, &b);
	    mp_mul_2d(&mminus, 1, &mminus);
	    if (m2plus > m2minus) {
		mp_mul_2d(&mplus, 1, &mplus);
	    err = mp_mul_2d(&b, 1, &b);
	    if (err == MP_OKAY) {
		err = mp_mul_2d(&mminus, 1, &mminus);
	    }
	    if ((err == MP_OKAY) && (m2plus > m2minus)) {
		err = mp_mul_2d(&mplus, 1, &mplus);
	    }
	    if (err == MP_OKAY) {
	    mp_div_d(&S, 5, &S, NULL);
		err = mp_div_d(&S, 5, &S, NULL);
	    }
	    --s5;

	    /*
	     * IDEA: It might possibly be a win to fall back to int64_t
	     *       arithmetic here if S < 2**64/10. But it's a win only for
	     *       a fairly narrow range of magnitudes so perhaps not worth
	     *       bothering.  We already know that we shorten the

	     * 10**38  12 trips
	     * 10**39  13 trips
	     * 10**40  14 trips
	     * 10**41  15 trips
	     * 10**42  16 trips
	     * thereafter no gain.
	     */
	} else {
	    mp_mul_d(&b, 10, &b);
	    mp_mul_d(&mminus, 10, &mminus);
	    if (m2plus > m2minus) {
		mp_mul_2d(&mplus, 10, &mplus);
	} else if (err == MP_OKAY) {
	    err = mp_mul_d(&b, 10, &b);
	    if (err == MP_OKAY) {
		err = mp_mul_d(&mminus, 10, &mminus);
	    }
	    if ((err == MP_OKAY) && (m2plus > m2minus)) {
		err = mp_mul_2d(&mplus, 10, &mplus);
	    }
	}

	++i;
    }

    /*

    char *s = retval;		/* Cursor in the return value. */
    mp_int b;			/* Numerator of the result. */
    mp_int S;			/* Denominator of the result. */
    mp_int dig;			/* Current digit of the result. */
    int digit;			/* Current digit of the result. */
    int g;			/* Size of the current digit ground. */
    int i, j;
    mp_err err;
    (void)dPtr;

    /*
     * b = bw * 2**b2 * 5**b5
     * S = 2**s2 * 5*s5
     */

    if (mp_init(&dig) != MP_OKAY) {
    mp_init_multi(&dig, NULL);
    mp_init_u64(&b, bw);
    mp_mul_2d(&b, b2, &b);
    mp_init_set(&S, 1);
    MulPow5(&S, s5, &S); mp_mul_2d(&S, s2, &S);
	return NULL;
    }
    if (mp_init_u64(&b, bw) != MP_OKAY) {
	mp_clear(&dig);
	return NULL;
    }
    err = mp_mul_2d(&b, b2, &b);
    if (err == MP_OKAY) {
 	err = mp_init_set(&S, 1);
    }
    if (err == MP_OKAY) {
	err = MulPow5(&S, s5, &S);
	if (err == MP_OKAY) {
	    err = mp_mul_2d(&S, s2, &S);
	}
    }

    /*
     * Handle the case where we guess the position of the decimal point wrong.
     */

    if (mp_cmp_mag(&b, &S) == MP_LT) {
    if ((mp_cmp_mag(&b, &S) == MP_LT) && (mp_mul_d(&b, 10, &b) == MP_OKAY)) {
	mp_mul_d(&b, 10, &b);
	ilim =ilim1;
	--k;
    }

    /*
     * Convert the leading digit.
     */

    i = 0;
    mp_div(&b, &S, &dig, &b);
    err = mp_div(&b, &S, &dig, &b);
    if (dig.used > 1 || dig.dp[0] >= 10) {
	Tcl_Panic("wrong digit!");
    }
    digit = dig.dp[0];

    /*
     * Is a single digit all that was requested?
     */

    *s++ = '0' + digit;
    if (++i >= ilim) {
	mp_mul_2d(&b, 1, &b);
	if (ShouldBankerRoundUp(&b, &S, digit&1)) {
	if ((mp_mul_2d(&b, 1, &b) == MP_OKAY) && ShouldBankerRoundUp(&b, &S, digit&1)) {
	    s = BumpUp(s, retval, &k);
	}
    } else {
	for (;;) {
	while (err == MP_OKAY) {
	    /*
	     * Shift by a group of digits.
	     */

	    g = ilim - i;
	    if (g > DIGIT_GROUP) {
		g = DIGIT_GROUP;
	    }
	    if (s5 >= g) {
		mp_div_d(&S, dpow5[g], &S, NULL);
		err = mp_div_d(&S, dpow5[g], &S, NULL);
		s5 -= g;
	    } else if (s5 > 0) {
		mp_div_d(&S, dpow5[s5], &S, NULL);
		mp_mul_d(&b, dpow5[g - s5], &b);
		err = mp_div_d(&S, dpow5[s5], &S, NULL);
		if (err == MP_OKAY) {
		    err = mp_mul_d(&b, dpow5[g - s5], &b);
		}
		s5 = 0;
	    } else {
		mp_mul_d(&b, dpow5[g], &b);
		err = mp_mul_d(&b, dpow5[g], &b);
	    }
	    if (err == MP_OKAY) {
	    mp_mul_2d(&b, g, &b);
		err = mp_mul_2d(&b, g, &b);
	    }

	    /*
	     * As with the shortening bignum conversion, it's possible at this
	     * point that we will have reduced the denominator to less than
	     * 2**64/10, at which point it would be possible to fall back to
	     * to int64_t arithmetic. But the potential payoff is tremendously
	     * less - unless we're working in F format - because we know that
	     * three groups of digits will always suffice for %#.17e, the
	     * longest format that doesn't introduce empty precision.
	     *
	     * Extract the next group of digits.
	     */

	    mp_div(&b, &S, &dig, &b);
	    if (dig.used > 1) {

	    if ((err != MP_OKAY) || (mp_div(&b, &S, &dig, &b) != MP_OKAY) || (dig.used > 1)) {
		Tcl_Panic("wrong digit!");
	    }
	    digit = dig.dp[0];
	    for (j = g-1; j >= 0; --j) {
		int t = itens[j];

		*s++ = digit / t + '0';
		digit %= t;
	    }
	    i += g;

	    /*
	     * Have we converted all the requested digits?
	     */

	    if (i == ilim) {
		mp_mul_2d(&b, 1, &b);
		if (ShouldBankerRoundUp(&b, &S, digit&1)) {
		if ((mp_mul_2d(&b, 1, &b) == MP_OKAY) && ShouldBankerRoundUp(&b, &S, digit&1)) {
		    s = BumpUp(s, retval, &k);
		}
		break;
	    }
	}
    }
    while (*--s == '0') {

    double d;
#ifdef IEEE_FLOATING_POINT
    union {
	double dv;
	Tcl_WideUInt iv;
    } bitwhack;
#endif
    mp_err err = MP_OKAY;
#if defined(__sgi) && defined(_COMPILER_VERSION)
    union fpc_csr mipsCR;

    mipsCR.fc_word = get_fpc_csr();
    mipsCR.fc_struct.flush = 0;
    set_fpc_csr(mipsCR.fc_word);
#endif

    /*
     * 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));
    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;

    }

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

    for (i=0; i<9; ++i) {
	mp_init(pow5 + i);
	err = err || mp_init(pow5 + i);
    }
    mp_set_u64(pow5, 5);
    for (i=0; i<8; ++i) {
	mp_sqr(pow5+i, pow5+i+1);
	err = err || mp_sqr(pow5+i, pow5+i+1);
    }
    mp_init_u64(pow5_13, 1220703125);
    err = err || mp_init_u64(pow5_13, 1220703125);
    for (i = 1; i < 5; ++i) {
	mp_init(pow5_13 + i);
	mp_sqr(pow5_13 + i - 1, pow5_13 + i);
	err = err || mp_init(pow5_13 + i);
	err = err || mp_sqr(pow5_13 + i - 1, pow5_13 + i);
    }
    if (err != MP_OKAY) {
	Tcl_Panic("out of memory");
    }

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

 *----------------------------------------------------------------------
 */

int
Tcl_InitBignumFromDouble(
    Tcl_Interp *interp,		/* For error message. */
    double d,			/* Number to convert. */
    mp_int *b)			/* Place to store the result. */
    void *big)			/* Place to store the result. */
{
    double fract;
    int expt;
    mp_err err;
    mp_int *b = (mp_int *)big;

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

    if (TclIsInfinite(d)) {
	if (interp != NULL) {
	    const 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);
	err = mp_init(b);
	mp_zero(b);
    } else {
	Tcl_WideInt w = (Tcl_WideInt) ldexp(fract, mantBits);
	int shift = expt - mantBits;

	mp_init_i64(b, w);
	if (shift < 0) {
	    mp_div_2d(b, -shift, b, NULL);
	err = mp_init_i64(b, w);
	if (err != MP_OKAY) {
		/* just skip */
	} else if (shift < 0) {
	    err = mp_div_2d(b, -shift, b, NULL);
	} else if (shift > 0) {
	    mp_mul_2d(b, shift, b);
	    err = mp_mul_2d(b, shift, b);
	}
    }
    if (err != MP_OKAY) {
	return TCL_ERROR;
    }
    return TCL_OK;
}

/*
 *----------------------------------------------------------------------
 *

 *	too large to convert.
 *
 *----------------------------------------------------------------------
 */

double
TclBignumToDouble(
    const mp_int *a)			/* Integer to convert. */
    const void *big)			/* Integer to convert. */
{
    mp_int b;
    int bits, shift, i, lsb;
    double r;
    mp_err err;
    const mp_int *a = (const mp_int *)big;


    /*
     * We need a 'mantBits'-bit significand.  Determine what shift will
     * give us that.
     */

     * in length.  If shift < 0, we will need to shift the significand right
     * by the requisite number of bits, and round it. If the '1-shift'
     * least significant bits are 0, but the 'shift'th bit is nonzero,
     * then the significand lies exactly between two values and must be
     * 'rounded to even'.
     */

    mp_init(&b);
    if (shift == 0) {
	mp_copy(a, &b);
    err = mp_init(&b);
    if (err != MP_OKAY) {
	/* just skip */
    } else if (shift == 0) {
	err = mp_copy(a, &b);
    } else if (shift > 0) {
	mp_mul_2d(a, shift, &b);
	err = mp_mul_2d(a, shift, &b);
    } else if (shift < 0) {
	lsb = mp_cnt_lsb(a);
	if (lsb == -1-shift) {

	    /*
	     * Round to even
	     */

	    mp_div_2d(a, -shift, &b, NULL);
	    if (mp_isodd(&b)) {
	    err = mp_div_2d(a, -shift, &b, NULL);
	    if ((err == MP_OKAY) && mp_isodd(&b)) {
		if (mp_isneg(&b)) {
		    mp_sub_d(&b, 1, &b);
		    err = mp_sub_d(&b, 1, &b);
		} else {
		    mp_add_d(&b, 1, &b);
		    err = mp_add_d(&b, 1, &b);
		}
	    }
	} else {

	    /*
	     * Ordinary rounding
	     */

	    mp_div_2d(a, -1-shift, &b, NULL);
	    if (mp_isneg(&b)) {
		mp_sub_d(&b, 1, &b);
	    err = mp_div_2d(a, -1-shift, &b, NULL);
	    if (err != MP_OKAY) {
		/* just skip */
	    } else if (mp_isneg(&b)) {
		err = mp_sub_d(&b, 1, &b);
	    } else {
		mp_add_d(&b, 1, &b);
		err = mp_add_d(&b, 1, &b);
	    }
	    mp_div_2d(&b, 1, &b, NULL);
	    err = mp_div_2d(&b, 1, &b, NULL);
	}
    }

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

    if (err != MP_OKAY) {
	return 0.0;
    }
    r = 0.0;
    for (i=b.used-1 ; i>=0 ; --i) {
    for (i = b.used-1; i>=0; --i) {
	r = ldexp(r, MP_DIGIT_BIT) + b.dp[i];
    }
    mp_clear(&b);

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

 *	Returns the floating point number.
 *
 *----------------------------------------------------------------------
 */

double
TclCeil(
    const mp_int *a)			/* Integer to convert. */
    const void *big)			/* Integer to convert. */
{
    double r = 0.0;
    mp_int b;
    mp_err err;
    const mp_int *a = (const mp_int *)big;

    mp_init(&b);
    if (mp_isneg(a)) {
	mp_neg(a, &b);
    err = mp_init(&b);
    if ((err == MP_OKAY) && mp_isneg(a)) {
	err = 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 (err != MP_OKAY) {
		/* just skip */
	    if (shift > 0) {
		mp_mul_2d(a, shift, &b);
	    } else if (shift > 0) {
		err = mp_mul_2d(a, shift, &b);
	    } else if (shift < 0) {
		mp_int d;
		mp_init(&d);
		mp_div_2d(a, -shift, &b, &d);
		exact = d.used == 0;
		err = mp_init(&d);
		if (err == MP_OKAY) {
		    err = mp_div_2d(a, -shift, &b, &d);
		}
		exact = mp_iszero(&d);
		mp_clear(&d);
	    } else {
		mp_copy(a, &b);
		err = mp_copy(a, &b);
	    }
	    if (!exact) {
		mp_add_d(&b, 1, &b);
	    if ((err == MP_OKAY) && !exact) {
		err = mp_add_d(&b, 1, &b);
	    }
	    if (err != MP_OKAY) {
		return 0.0;
	    }
	    for (i=b.used-1 ; i>=0 ; --i) {
		r = ldexp(r, MP_DIGIT_BIT) + b.dp[i];
	    }
	    r = ldexp(r, bits - mantBits);
	}
    }

 *	Returns the floating point value.
 *
 *----------------------------------------------------------------------
 */

double
TclFloor(
    const mp_int *a)			/* Integer to convert. */
    const void *big)			/* Integer to convert. */
{
    double r = 0.0;
    mp_int b;
    mp_err err;
    const mp_int *a = (const mp_int *)big;

    mp_init(&b);
    if (mp_isneg(a)) {
	mp_neg(a, &b);
    err = mp_init(&b);
    if ((err == MP_OKAY) && mp_isneg(a)) {
	err = 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);
		err = mp_mul_2d(a, shift, &b);
	    } else if (shift < 0) {
		mp_div_2d(a, -shift, &b, NULL);
		err = mp_div_2d(a, -shift, &b, NULL);
	    } else {
		mp_copy(a, &b);
		err = mp_copy(a, &b);
	    }
	    if (err != MP_OKAY) {
		return 0.0;
	    }
	    for (i=b.used-1 ; i>=0 ; --i) {
		r = ldexp(r, MP_DIGIT_BIT) + b.dp[i];
	    }
	    r = ldexp(r, bits - mantBits);
	}
    }

    int *machexp)		/* Power of two. */
{
    mp_int b;
    int bits;
    int shift;
    int i;
    double r;
    mp_err err = MP_OKAY;

    /*
     * 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 (mp_init(&b)) {
	return 0.0;
    }
    if (shift > 0) {
	mp_mul_2d(a, shift, &b);
	err = mp_mul_2d(a, shift, &b);
    } else if (shift < 0) {
	mp_div_2d(a, -shift, &b, NULL);
	err = mp_div_2d(a, -shift, &b, NULL);
    } else {
	mp_copy(a, &b);
	err = mp_copy(a, &b);
    }

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

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

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