/*
 * Static functions defined in this file.
 */

static int		AccumulateDecimalDigit(unsigned, int,
			    Tcl_WideUInt *, mp_int *, int);
static double		MakeHighPrecisionDouble(int signum,
			    mp_int *significand, int nSigDigs, long exponent);
static double		MakeLowPrecisionDouble(int signum,
			    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 int 		NormalizeRightward(Tcl_WideUInt *);

	     * to whether 'significandOverflow' is set. The desired floating
	     * point value is significand * 10**k, where
	     * k = numTrailZeros+exponent-numDigitsAfterDp.
	     */

	    objPtr->typePtr = &tclDoubleType;
	    if (exponentSignum) {
		/* 
		 * At this point exponent>=0, so the following calculation
		 * cannot underflow. 
		 */
		exponent = -exponent;
	    }

	    /*
	     * Adjust the exponent for the number of trailing zeros that
	     * have not been accumulated, and the number of digits after
	     * the decimal point. Pin any overflow to LONG_MAX/LONG_MIN
	     * respectively.
	     */

	    if (exponent >= 0) {
		if (exponent - numDigitsAfterDp > LONG_MAX - numTrailZeros) {
		    exponent = LONG_MAX;
		} else {
		    exponent = exponent - numDigitsAfterDp + numTrailZeros;
		}
	    } else {
		if (exponent + numTrailZeros < LONG_MIN + numDigitsAfterDp) {
		    exponent = LONG_MIN;
		} else {
		    exponent = exponent + numTrailZeros - numDigitsAfterDp;
		}
	    }

	    /* 
	     * The desired number is now significandWide * 10**exponent
	     * or significandBig * 10**exponent, depending on whether
	     * the significand has overflowed a wide int.
	     */
	    if (!significandOverflow) {
		objPtr->internalRep.doubleValue = MakeLowPrecisionDouble(
			signum, significandWide, numSigDigs, exponent);

	    } else {
		objPtr->internalRep.doubleValue = MakeHighPrecisionDouble(
			signum, &significandBig, numSigDigs, exponent);

	    }
	    break;

	case sINF:
	case sINFINITY:
	    if (signum) {
		objPtr->internalRep.doubleValue = -HUGE_VAL;

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

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

    TCL_IEEE_DOUBLE_ROUNDING;

    /*
     * Test for the easy cases.
     */

    if (significand == 0) {
	return copysign(0.0, -signum);
    }
    if (numSigDigs <= QUICK_MAX) {
	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.

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

    TCL_IEEE_DOUBLE_ROUNDING;

    /*
     * Quick checks for zero, and over/underflow. Be careful to avoid
     * integer overflow when calculating with 'exponent'.
     */

    if (mp_iszero(significand)) {
	return copysign(0.0, -signum);
    }
    if (exponent >= 0 && exponent-1 > maxDigits-numSigDigs) {
	retval = HUGE_VAL;
	goto returnValue;

    } else if (exponent < 0 && numSigDigs+exponent < minDigits+1) {
	retval = 0.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

    double quot;		/* Correction term. */
    double minincr;		/* Lower bound on the absolute value of the
				 * correction term. */
    int roundToEven = 0;	/* Flag == TRUE if we need to invoke
				 * "round to even" functionality */
    double rteSignificand;	/* Significand of the round-to-even result */
    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;

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

    if (approxResult == HUGE_VAL) {
	return approxResult;
    }
    significand = frexp(approxResult, &binExponent);

    /*

     * We are trying to compute a corrector term that, when added to the
     * approximate result, will yield close to the exact result.
     * The exact result is exactSignificand * 10**exponent.
     * The approximate result is significand * 2**binExponent
     * If exponent<0, we need to multiply the exact value by 10**-exponent
     * to make it an integer, plus another factor of 2 to decide on rounding.
     *  Similarly if binExponent<FP_PRECISION, we need
     * to multiply by 2**FP_PRECISION to make the approximate value an integer.
     *
     * Let M = 2**M2 * 5**M5 be the least common multiple of these two
     * multipliers.
     */


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

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

    /*
     * 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 ((M5 + exponent) & (1 << i)) {
	    mp_mul(&twoMd, pow5+i, &twoMd);
	}
    }
    mp_mul_2d(&twoMd, M2+exponent+1, &twoMd);

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

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

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

    switch (mp_cmp_mag(&twoMd, &twoMv)) {
    case MP_LT:
	/*

	 * If the error is less than 1/2 ULP, there's no correction to make.
	 */
	mp_clear(&twoMd);
	mp_clear(&twoMv);
	return approxResult;
    case MP_EQ:
	/*
	 * If the error is exactly 1/2 ULP, we need to round to even.
	 */
	roundToEven = 1;
	break;
    case MP_GT:
	/*
	 * We need to correct the result if the error exceeds 1/2 ULP.
	 */
	break;
    }

    /*
     * If we're in the 'round to even' case, and the significand is already
     * even, we're done. Return the approximate result.
     */
    if (roundToEven) {
	rteSignificand = frexp(approxResult, &rteExponent);
	rteSigWide = (Tcl_WideInt) ldexp(rteSignificand, FP_PRECISION);
	if ((rteSigWide & 1) == 0) {
	    mp_clear(&twoMd);
	    mp_clear(&twoMv);
	    return approxResult;
	}
    }

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

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

    num = TclBignumToDouble(&twoMd);