diff options
| -rw-r--r-- | generic/tclStrToD.c | 142 | ||||
| -rw-r--r-- | tests/expr.test | 66 |
2 files changed, 172 insertions, 36 deletions
diff --git a/generic/tclStrToD.c b/generic/tclStrToD.c index eb71276..236fe59 100644 --- a/generic/tclStrToD.c +++ b/generic/tclStrToD.c @@ -289,10 +289,10 @@ static const Tcl_WideUInt wuipow5[27] = { static int AccumulateDecimalDigit(unsigned, int, Tcl_WideUInt *, mp_int *, int); static double MakeHighPrecisionDouble(int signum, - mp_int *significand, int nSigDigs, int exponent); + mp_int *significand, int nSigDigs, long exponent); static double MakeLowPrecisionDouble(int signum, Tcl_WideUInt significand, int nSigDigs, - int exponent); + long exponent); #ifdef IEEE_FLOATING_POINT static double MakeNaN(int signum, Tcl_WideUInt tag); #endif @@ -1338,16 +1338,45 @@ TclParseNumber( 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, - numTrailZeros + exponent - numDigitsAfterDp); + signum, significandWide, numSigDigs, exponent); } else { objPtr->internalRep.doubleValue = MakeHighPrecisionDouble( - signum, &significandBig, numSigDigs, - numTrailZeros + exponent - numDigitsAfterDp); + signum, &significandBig, numSigDigs, exponent); } break; @@ -1531,9 +1560,9 @@ AccumulateDecimalDigit( 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. */ + 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. */ @@ -1552,6 +1581,9 @@ MakeLowPrecisionDouble( * Test for the easy cases. */ + if (significand == 0) { + return copysign(0.0, -signum); + } if (numSigDigs <= QUICK_MAX) { if (exponent >= 0) { if (exponent <= mmaxpow) { @@ -1644,10 +1676,10 @@ MakeLowPrecisionDouble( 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. */ + 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. */ @@ -1663,15 +1695,18 @@ MakeHighPrecisionDouble( TCL_IEEE_DOUBLE_ROUNDING; /* - * Quick checks for over/underflow. + * Quick checks for zero, and over/underflow. Be careful to avoid + * integer overflow when calculating with 'exponent'. */ - if (numSigDigs+exponent-1 > maxDigits) { + if (mp_iszero(significand)) { + return copysign(0.0, -signum); + } + if (exponent >= 0 && exponent-1 > maxDigits-numSigDigs) { retval = HUGE_VAL; goto returnValue; - } - if (numSigDigs+exponent-1 < minDigits) { - retval = 0; + } else if (exponent < 0 && numSigDigs+exponent < minDigits+1) { + retval = 0.0; goto returnValue; } @@ -1806,6 +1841,9 @@ RefineApproximation( * "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; @@ -1818,13 +1856,22 @@ RefineApproximation( if (approxResult == HUGE_VAL) { return approxResult; } + significand = frexp(approxResult, &binExponent); /* - * Find a common denominator for the decimal and binary fractions. The - * common denominator will be 2**M2 + 5**M5. + * 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. */ - significand = frexp(approxResult, &binExponent); i = mantBits - binExponent; if (i < 0) { M2 = 0; @@ -1841,12 +1888,9 @@ RefineApproximation( } /* - * 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. + * 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 */ @@ -1867,10 +1911,9 @@ RefineApproximation( } /* - * 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). + * 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); @@ -1880,16 +1923,22 @@ RefineApproximation( } } 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. + * 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)) { @@ -1903,25 +1952,36 @@ RefineApproximation( 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 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 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 result is exactly unity, we need to round to even. + * 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); @@ -1932,6 +1992,16 @@ RefineApproximation( } } + /* + * 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. diff --git a/tests/expr.test b/tests/expr.test index bc01c03..4724ee1 100644 --- a/tests/expr.test +++ b/tests/expr.test @@ -6822,6 +6822,72 @@ test expr-41.2 {exponent underflow} { expr 1.0e-2147483630 } 0.0 +test expr-41.3 {exponent overflow} { + expr 1e2147483647 +} Inf +test expr-41.4 {exponent overflow} { + expr 1e2147483648 +} Inf +test expr-41.5 {exponent overflow} { + expr 100e2147483645 +} Inf +test expr-41.6 {exponent overflow} { + expr 100e2147483646 +} Inf +test expr-41.7 {exponent overflow} { + expr 1.0e2147483647 +} Inf +test expr-41.8 {exponent overflow} { + expr 1.0e2147483648 +} Inf +test expr-41.9 {exponent overflow} { + expr 1.2e2147483647 +} Inf +test expr-41.10 {exponent overflow} { + expr 1.2e2147483648 +} Inf + +test expr-41.11 {exponent overflow} { + expr 1e-2147483648 +} 0.0 +test expr-41.12 {exponent overflow} { + expr 1e-2147483649 +} 0.0 +test expr-41.13 {exponent overflow} { + expr 100e-2147483650 +} 0.0 +test expr-41.14 {exponent overflow} { + expr 100e-2147483651 +} 0.0 +test expr-41.15 {exponent overflow} { + expr 1.0e-2147483648 +} 0.0 +test expr-41.16 {exponent overflow} { + expr 1.0e-2147483649 +} 0.0 +test expr-41.17 {exponent overflow} { + expr 1.23e-2147483646 +} 0.0 +test expr-41.18 {exponent overflow} { + expr 1.23e-2147483647 +} 0.0 + +test expr-41.19 {numSigDigs == 0} { + expr 0e309 +} 0.0 +test expr-41.20 {numSigDigs == 0} { + expr 0e310 +} 0.0 +test expr-41.21 {negative zero, large exponent} { + expr -0e309 +} -0.0 +test expr-41.22 {negative zero, large exponent} { + expr -0e310 +} -0.0 +test expr-41.23 {floating point overflow on significand (Bug 1de6b0629e)} { + expr 123[string repeat 0 309]1e-310 +} 123.0 + test expr-42.1 {denormals} ieeeFloatingPoint { expr 7e-324 } 5e-324 |