summaryrefslogtreecommitdiffstats
path: root/generic/tclStrToD.c
diff options
context:
space:
mode:
Diffstat (limited to 'generic/tclStrToD.c')
-rw-r--r--generic/tclStrToD.c201
1 files changed, 139 insertions, 62 deletions
diff --git a/generic/tclStrToD.c b/generic/tclStrToD.c
index d624f87..e0c9f15 100644
--- a/generic/tclStrToD.c
+++ b/generic/tclStrToD.c
@@ -14,9 +14,14 @@
*/
#include "tclInt.h"
-#include "tommath.h"
+#include "tclTomMath.h"
+#include <float.h>
#include <math.h>
+#ifdef _WIN32
+#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
@@ -289,10 +294,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
@@ -706,7 +711,7 @@ TclParseNumber(
|| (octalSignificandWide >
((Tcl_WideUInt)-1 >> shift)))) {
octalSignificandOverflow = 1;
- mp_init_ull(&octalSignificandBig,
+ mp_init_u64(&octalSignificandBig,
octalSignificandWide);
}
}
@@ -771,7 +776,7 @@ TclParseNumber(
((size_t)shift >= CHAR_BIT*sizeof(Tcl_WideUInt) ||
significandWide > ((Tcl_WideUInt)-1 >> shift))) {
significandOverflow = 1;
- mp_init_ull(&significandBig,
+ mp_init_u64(&significandBig,
significandWide);
}
}
@@ -790,6 +795,7 @@ TclParseNumber(
acceptState = state;
acceptPoint = p;
acceptLen = len;
+ /* FALLTHRU */
case ZERO_B:
zerob:
if (c == '0') {
@@ -812,7 +818,7 @@ TclParseNumber(
((size_t)shift >= CHAR_BIT*sizeof(Tcl_WideUInt) ||
significandWide > ((Tcl_WideUInt)-1 >> shift))) {
significandOverflow = 1;
- mp_init_ull(&significandBig,
+ mp_init_u64(&significandBig,
significandWide);
}
}
@@ -1154,7 +1160,7 @@ TclParseNumber(
((size_t)shift >= CHAR_BIT*sizeof(Tcl_WideUInt) ||
significandWide > (MOST_BITS + signum) >> shift)) {
significandOverflow = 1;
- mp_init_ull(&significandBig, significandWide);
+ mp_init_u64(&significandBig, significandWide);
}
if (shift) {
if (!significandOverflow) {
@@ -1175,7 +1181,7 @@ TclParseNumber(
((size_t)shift >= CHAR_BIT*sizeof(Tcl_WideUInt) ||
significandWide > (MOST_BITS + signum) >> shift)) {
significandOverflow = 1;
- mp_init_ull(&significandBig, significandWide);
+ mp_init_u64(&significandBig, significandWide);
}
if (shift) {
if (!significandOverflow) {
@@ -1196,7 +1202,7 @@ TclParseNumber(
((size_t)shift >= CHAR_BIT*sizeof(Tcl_WideUInt) ||
octalSignificandWide > (MOST_BITS + signum) >> shift)) {
octalSignificandOverflow = 1;
- mp_init_ull(&octalSignificandBig,
+ mp_init_u64(&octalSignificandBig,
octalSignificandWide);
}
if (shift) {
@@ -1209,7 +1215,7 @@ TclParseNumber(
}
if (!octalSignificandOverflow) {
if (octalSignificandWide > (MOST_BITS + signum)) {
- mp_init_ull(&octalSignificandBig,
+ mp_init_u64(&octalSignificandBig,
octalSignificandWide);
octalSignificandOverflow = 1;
} else {
@@ -1225,7 +1231,7 @@ TclParseNumber(
}
if (octalSignificandOverflow) {
if (signum) {
- mp_neg(&octalSignificandBig, &octalSignificandBig);
+ (void)mp_neg(&octalSignificandBig, &octalSignificandBig);
}
TclSetBignumIntRep(objPtr, &octalSignificandBig);
}
@@ -1235,14 +1241,14 @@ TclParseNumber(
case DECIMAL:
significandOverflow = AccumulateDecimalDigit(0, numTrailZeros-1,
&significandWide, &significandBig, significandOverflow);
- if (!significandOverflow && (significandWide > MOST_BITS+signum)){
+ if (!significandOverflow && (significandWide > MOST_BITS+signum)) {
significandOverflow = 1;
- mp_init_ull(&significandBig, significandWide);
+ mp_init_u64(&significandBig, significandWide);
}
returnInteger:
if (!significandOverflow) {
if (significandWide > MOST_BITS+signum) {
- mp_init_ull(&significandBig,
+ mp_init_u64(&significandBig,
significandWide);
significandOverflow = 1;
} else {
@@ -1258,7 +1264,7 @@ TclParseNumber(
}
if (significandOverflow) {
if (signum) {
- mp_neg(&significandBig, &significandBig);
+ (void)mp_neg(&significandBig, &significandBig);
}
TclSetBignumIntRep(objPtr, &significandBig);
}
@@ -1277,16 +1283,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;
@@ -1303,7 +1338,7 @@ TclParseNumber(
#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
@@ -1394,7 +1429,7 @@ AccumulateDecimalDigit(
* bignum and fall through into the bignum case.
*/
- mp_init_ull(bignumRepPtr, w);
+ mp_init_u64(bignumRepPtr, w);
} else {
/*
* Wide multiplication.
@@ -1467,9 +1502,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. */
@@ -1488,6 +1523,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) {
@@ -1537,7 +1575,7 @@ MakeLowPrecisionDouble(
* call MakeHighPrecisionDouble to do it the hard way.
*/
- mp_init_ull(&significandBig, significand);
+ mp_init_u64(&significandBig, significand);
retval = MakeHighPrecisionDouble(0, &significandBig, numSigDigs,
exponent);
mp_clear(&significandBig);
@@ -1580,10 +1618,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. */
@@ -1599,15 +1637,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;
}
@@ -1742,6 +1783,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;
@@ -1754,13 +1798,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;
@@ -1777,12 +1830,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 */
@@ -1803,10 +1853,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);
@@ -1816,17 +1865,23 @@ 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(&twoMv, 1);
+ mp_set_u64(&twoMv, 1);
for (i=0; i<=8; ++i) {
if (M5 & (1 << i)) {
mp_mul(&twoMv, pow5+i, &twoMv);
@@ -1839,25 +1894,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);
@@ -1869,6 +1935,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.
*/
@@ -3204,7 +3280,7 @@ ShorteningBignumConversionPowD(
* mminus = 5**m5
*/
- mp_init_ull(&b, bw);
+ mp_init_u64(&b, bw);
mp_init_set(&mminus, 1);
MulPow5(&b, b5, &b);
mp_mul_2d(&b, b2, &b);
@@ -3388,7 +3464,7 @@ StrictBignumConversionPowD(
* b = bw * 2**b2 * 5**b5
*/
- mp_init_ull(&b, bw);
+ mp_init_u64(&b, bw);
MulPow5(&b, b5, &b);
mp_mul_2d(&b, b2, &b);
@@ -3588,7 +3664,7 @@ ShorteningBignumConversion(
* S = 2**s2 * 5*s5
*/
- mp_init_ull(&b, bw);
+ 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);
@@ -3797,7 +3873,7 @@ StrictBignumConversion(
*/
mp_init_multi(&dig, NULL);
- mp_init_ull(&b, bw);
+ 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);
@@ -4274,7 +4350,8 @@ TclInitDoubleConversion(void)
maxpow10_wide = (int)
floor(sizeof(Tcl_WideUInt) * CHAR_BIT * log(2.) / log(10.));
- pow10_wide = Tcl_Alloc((maxpow10_wide + 1) * sizeof(Tcl_WideUInt));
+ pow10_wide = (Tcl_WideUInt *)
+ Tcl_Alloc((maxpow10_wide + 1) * sizeof(Tcl_WideUInt));
u = 1;
for (i = 0; i < maxpow10_wide; ++i) {
pow10_wide[i] = u;
@@ -4317,11 +4394,11 @@ TclInitDoubleConversion(void)
for (i=0; i<9; ++i) {
mp_init(pow5 + i);
}
- mp_set(pow5, 5);
+ mp_set_u64(pow5, 5);
for (i=0; i<8; ++i) {
mp_sqr(pow5+i, pow5+i+1);
}
- mp_init_ul(pow5_13, 1220703125);
+ 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);
@@ -4439,7 +4516,7 @@ Tcl_InitBignumFromDouble(
Tcl_WideInt w = (Tcl_WideInt) ldexp(fract, mantBits);
int shift = expt - mantBits;
- mp_init_ll(b, w);
+ mp_init_i64(b, w);
if (shift < 0) {
mp_div_2d(b, -shift, b, NULL);
} else if (shift > 0) {
@@ -4602,7 +4679,7 @@ TclCeil(
mp_int d;
mp_init(&d);
mp_div_2d(a, -shift, &b, &d);
- exact = d.used == 0;
+ exact = mp_iszero(&d);
mp_clear(&d);
} else {
mp_copy(a, &b);
generated by cgit v0.12 at 2025-12-15 04:53:41 (GMT)