summaryrefslogtreecommitdiffstats
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)
commitaf46fd5b12c1265e04beeab0871ce9207e86c7a9 (patch)
treebbd9845eb037231f2b5cb390052f77e653fd6ceb
parent69ead1ea0999daf76ab2e38733cdd63fc979ad8f (diff)
downloadtcl-af46fd5b12c1265e04beeab0871ce9207e86c7a9.zip
tcl-af46fd5b12c1265e04beeab0871ce9207e86c7a9.tar.gz
tcl-af46fd5b12c1265e04beeab0871ce9207e86c7a9.tar.bz2
Make sure that an input lying precisely 1/2 ULP between two floating point values is rounded to even.
-rw-r--r--generic/tclStrToD.c39
-rw-r--r--tests/expr.test4
2 files changed, 35 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.
diff --git a/tests/expr.test b/tests/expr.test
index 6ad7208..4c03262 100644
--- a/tests/expr.test
+++ b/tests/expr.test
@@ -7174,6 +7174,10 @@ test expr-50.1 {test sqrt() of bignums with non-Inf answer} {
expr {sqrt("1[string repeat 0 616]") == 1e308}
} 1
+test expr-51.1 {test round-to-even on input} {
+ expr 6.9294956446009195e15
+} 6929495644600920.0
+
# cleanup