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

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.