summaryrefslogtreecommitdiffstats
path: root/generic/tclStrToD.c
diff options
context:
space:
mode:
authordgp <dgp@users.sourceforge.net>2015-06-24 13:55:46 (GMT)
committerdgp <dgp@users.sourceforge.net>2015-06-24 13:55:46 (GMT)
commit448da15b7fffbb7e10b589a22c7ee059a629bf60 (patch)
tree5ccaf989644eb1335f2de7a6ce5b0c14f471ec70 /generic/tclStrToD.c
parent447722a6ad983008b5be213d3e39685262398019 (diff)
downloadtcl-448da15b7fffbb7e10b589a22c7ee059a629bf60.zip
tcl-448da15b7fffbb7e10b589a22c7ee059a629bf60.tar.gz
tcl-448da15b7fffbb7e10b589a22c7ee059a629bf60.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.c38
1 files changed, 30 insertions, 8 deletions
diff --git a/generic/tclStrToD.c b/generic/tclStrToD.c
index 76adf75..ec5e764 100644
--- a/generic/tclStrToD.c
+++ b/generic/tclStrToD.c
@@ -1515,7 +1515,7 @@ MakeLowPrecisionDouble(
* Test for the easy cases.
*/
- if (numSigDigs <= DBL_DIG) {
+ if (numSigDigs <= QUICK_MAX) {
if (exponent >= 0) {
if (exponent <= mmaxpow) {
/*
@@ -1527,7 +1527,7 @@ MakeLowPrecisionDouble(
retval = (double)(Tcl_WideInt)significand * pow10vals[exponent];
goto returnValue;
} else {
- int diff = DBL_DIG - numSigDigs;
+ int diff = QUICK_MAX - numSigDigs;
if (exponent-diff <= mmaxpow) {
/*
* 10**exponent is not an exact integer, but
@@ -1779,6 +1779,12 @@ RefineApproximation(
double quot; /* Correction term */
double minincr; /* Lower bound on the absolute value of the
* correction term. */
+ int roundToEven = 0; /* 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;
/*
@@ -1874,17 +1880,33 @@ 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:
+ 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.