*
 * 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 <math.h>

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

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

	    }
	    break;

	case ZERO:
	case DECIMAL:
	    significandOverflow = AccumulateDecimalDigit(0, numTrailZeros-1,
		    &significandWide, &significandBig, significandOverflow);
	    if (!significandOverflow && (significandWide > MOST_BITS+signum)){
		significandOverflow = 1;
		mp_init_u64(&significandBig, significandWide);
	    }
	returnInteger:
	    if (!significandOverflow) {
		if (significandWide > MOST_BITS+signum) {
		    mp_init_u64(&significandBig,

	     * 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

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

	    }
	    objPtr->typePtr = &tclDoubleType;
	    break;

#ifdef IEEE_FLOATING_POINT
	case sNAN:
	case sNANFINISH:
	    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");
	}

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

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

    /*

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

    /*
     * 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_Alloc((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;

	    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 = d.used == 0;
		mp_clear(&d);
	    } else {
		mp_copy(a, &b);
	    }
	    if (!exact) {
		mp_add_d(&b, 1, &b);
	    }