summaryrefslogtreecommitdiffstats
path: root/generic/tclStrToD.c
diff options
context:
space:
mode:
authorKevin B Kenny <kennykb@acm.org>2015-06-20 20:44:17 (GMT)
committerKevin B Kenny <kennykb@acm.org>2015-06-20 20:44:17 (GMT)
commit9c52f7a98a7419823cdb5387061db6b887131cec (patch)
treebbd9845eb037231f2b5cb390052f77e653fd6ceb /generic/tclStrToD.c
parent0985259acdc0129dd9e410e90262ba79d1fa76fe (diff)
downloadtcl-9c52f7a98a7419823cdb5387061db6b887131cec.zip
tcl-9c52f7a98a7419823cdb5387061db6b887131cec.tar.gz
tcl-9c52f7a98a7419823cdb5387061db6b887131cec.tar.bz2
Make sure that an input lying precisely 1/2 ULP between two floating point values is rounded to even.
Diffstat (limited to 'generic/tclStrToD.c')
-rw-r--r--generic/tclStrToD.c39
1 files changed, 31 insertions, 8 deletions
diff --git a/generic/tclStrToD.c b/generic/tclStrToD.c
index 883e2ea..c1b3e99 100644
--- a/generic/tclStrToD.c
+++ b/generic/tclStrToD.c
@@ -1548,7 +1548,7 @@ MakeLowPrecisionDouble(
* Test for the easy cases.
*/
- if (numSigDigs <= DBL_DIG) {
+ if (numSigDigs <= QUICK_MAX) {
if (exponent >= 0) {
if (exponent <= mmaxpow) {
/*
@@ -1561,7 +1561,7 @@ MakeLowPrecisionDouble(
((Tcl_WideInt)significand * pow10vals[exponent]);
goto returnValue;
} else {
- int diff = DBL_DIG - numSigDigs;
+ int diff = QUICK_MAX - numSigDigs;
if (exponent-diff <= mmaxpow) {
/*
@@ -1798,6 +1798,12 @@ RefineApproximation(
double quot; /* Correction term. */
double minincr; /* Lower bound on the absolute value of the
* correction term. */
+ int roundToEven; /* 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 */
+ Tcl_WideInt rteSigWide; /* Wide integer version of the significand
+ * for testing evenness */
int i;
/*
@@ -1893,17 +1899,34 @@ RefineApproximation(
mp_div_2d(&twoMv, -multiplier, &twoMv, NULL);
}
- /*
- * 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 (mp_cmp_mag(&twoMd, &twoMv) == MP_LT) {
+ 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.
+ */
mp_clear(&twoMd);
mp_clear(&twoMv);
return approxResult;
+ case MP_EQ:
+ /*
+ * If the result is exactly unity, we need to round to even.
+ */
+ roundToEven = 1;
+ break;
+ case MP_GT:
+ roundToEven = 0;
+ break;
}
+ if (roundToEven) {
+ rteSignificand = frexp(approxResult, &rteExponent);
+ rteSigWide = (Tcl_WideInt) ldexp(rteSignificand, FP_PRECISION);
+ if ((rteSigWide & 1) == 0) {
+ return approxResult;
+ }
+ }
+
/*
* Convert the numerator and denominator of the corrector term accurately
* to floating point numbers.