Make sure that an input lying precisely 1/2 ULP between two floating point values is rounded to even.

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.