From 3ea3c1f2f99b0ca58203368be9c7baeff3f6c5e4 Mon Sep 17 00:00:00 2001
From: dgp <dgp@users.sourceforge.net>
Date: Wed, 24 Jun 2015 13:55:46 +0000
Subject: Make sure that an input lying precisely 1/2 ULP between two floating
 point values is rounded to even.

---
 generic/tclStrToD.c | 38 ++++++++++++++++++++++++++++++--------
 tests/expr.test     |  8 ++++++++
 2 files changed, 38 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.
diff --git a/tests/expr.test b/tests/expr.test
index 42d0c79..bd5ed8f 100644
--- a/tests/expr.test
+++ b/tests/expr.test
@@ -7165,6 +7165,14 @@ test expr-48.1 {Bug 1770224} {
     expr {-0x8000000000000001 >> 0x8000000000000000}
 } -1
 
+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
 if {[info exists a]} {
     unset a
-- 
cgit v0.12