Fix roundoff gaffe in bignum-to-double conversion [Bug 3349507]

12 files changed, 115 insertions, 12 deletions

diff --git a/ChangeLog b/ChangeLog
index 5bd1462..f12e4ac 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,18 @@
+2011-07-02  Kevin B. Kenny  <kennykb@acm.org>
+
+	* generic/tclStrToD.c:
+	* generic/tclTomMath.decls:
+	* generic/tclTomMathDecls.h:
+	* macosx/Tcl.xcode/project.pbxproj:
+	* macosx/Tcl.xcodeproj/project.pbxproj:
+	* tests/util.test:
+	* unix/Makefile.in:
+	* win/Makefile.in:
+	* win/Makefile.vc:
+	Fix a bug where bignum->double conversion is "round up" and
+	not "round to nearest" (causing expr double(1[string repeat 0 23])
+	not to be 1e+23). [Bug 3349507]
+
 2011-06-28  Reinhard Max  <max@suse.de>
 
 	* unix/tclUnixSock.c (CreateClientSocket): Fix and simplify
diff --git a/generic/tclStrToD.c b/generic/tclStrToD.c
index 16f11d2..15bff3e 100755
--- a/generic/tclStrToD.c
+++ b/generic/tclStrToD.c
@@ -4563,12 +4563,13 @@ TclBignumToDouble(
     const mp_int *a)			/* Integer to convert. */
 {
     mp_int b;
-    int bits, shift, i;
+    int bits, shift, i, lsb;
     double r;
 
+
     /*
-     * Determine how many bits we need, and extract that many from the input.
-     * Round to nearest unit in the last place.
+     * We need a 'mantBits'-bit significand.  Determine what shift will 
+     * give us that.
      */
 
     bits = mp_count_bits(a);
@@ -4580,17 +4581,54 @@ TclBignumToDouble(
 	    return -HUGE_VAL;
 	}
     }
-    shift = mantBits + 1 - bits;
+    shift = mantBits - bits;
+
+    /* 
+     * If shift > 0, shift the significand left by the requisite number of
+     * bits.  If shift == 0, the significand is already exactly 'mantBits'
+     * in length.  If shift < 0, we will need to shift the significand right
+     * by the requisite number of bits, and round it. If the '1-shift'
+     * least significant bits are 0, but the 'shift'th bit is nonzero,
+     * then the significand lies exactly between two values and must be
+     * 'rounded to even'.
+     */
+
     mp_init(&b);
-    if (shift > 0) {
+    if (shift == 0) {
+	mp_copy(a, &b);
+    } else if (shift > 0) {
 	mp_mul_2d(a, shift, &b);
     } else if (shift < 0) {
-	mp_div_2d(a, -shift, &b, NULL);
-    } else {
-	mp_copy(a, &b);
+	lsb = mp_cnt_lsb(a);
+	if (lsb == -1-shift) {
+
+	    /*
+	     * Round to even
+	     */
+
+	    mp_div_2d(a, -shift, &b, NULL);
+	    if (mp_isodd(&b)) {
+		if (b.sign == MP_ZPOS) {
+		    mp_add_d(&b, 1, &b);
+		} else {
+		    mp_sub_d(&b, 1, &b);
+		}
+	    }
+	} else {
+
+	    /*
+	     * Ordinary rounding
+	     */
+
+	    mp_div_2d(a, -1-shift, &b, NULL);
+	    if (b.sign == MP_ZPOS) {
+		mp_add_d(&b, 1, &b);
+	    } else {
+		mp_sub_d(&b, 1, &b);
+	    }
+	    mp_div_2d(&b, 1, &b, NULL);
+	}
     }
-    mp_add_d(&b, 1, &b);
-    mp_div_2d(&b, 1, &b, NULL);
 
     /*
      * Accumulate the result, one mp_digit at a time.
diff --git a/generic/tclStubInit.c b/generic/tclStubInit.c
index 054ece5..5b47e95 100644
--- a/generic/tclStubInit.c
+++ b/generic/tclStubInit.c
@@ -464,6 +464,7 @@ const TclTomMathStubs tclTomMathStubs = {
     TclBN_s_mp_sub, /* 60 */
     TclBN_mp_init_set_int, /* 61 */
     TclBN_mp_set_int, /* 62 */
+    TclBN_mp_cnt_lsb, /* 63 */
 };
 
 static const TclStubHooks tclStubHooks = {
diff --git a/generic/tclTomMath.decls b/generic/tclTomMath.decls
index 5bf338e..29a6a03 100644
--- a/generic/tclTomMath.decls
+++ b/generic/tclTomMath.decls
@@ -218,3 +218,6 @@ declare 61 {
 declare 62 {
     int TclBN_mp_set_int(mp_int* a, unsigned long i)
 }
+declare 63 {
+    int TclBN_mp_cnt_lsb(const mp_int* a)
+}
diff --git a/generic/tclTomMathDecls.h b/generic/tclTomMathDecls.h
index 7df0d90..feaefb3 100644
--- a/generic/tclTomMathDecls.h
+++ b/generic/tclTomMathDecls.h
@@ -63,6 +63,7 @@
 #define mp_cmp TclBN_mp_cmp
 #define mp_cmp_d TclBN_mp_cmp_d
 #define mp_cmp_mag TclBN_mp_cmp_mag
+#define mp_cnt_lsb TclBN_mp_cnt_lsb
 #define mp_copy TclBN_mp_copy
 #define mp_count_bits TclBN_mp_count_bits
 #define mp_div TclBN_mp_div
@@ -272,6 +273,8 @@ EXTERN int		TclBN_s_mp_sub(mp_int *a, mp_int *b, mp_int *c);
 EXTERN int		TclBN_mp_init_set_int(mp_int*a, unsigned long i);
 /* 62 */
 EXTERN int		TclBN_mp_set_int(mp_int*a, unsigned long i);
+/* 63 */
+EXTERN int		TclBN_mp_cnt_lsb(const mp_int*a);
 
 typedef struct TclTomMathStubs {
     int magic;
@@ -340,6 +343,7 @@ typedef struct TclTomMathStubs {
     int (*tclBN_s_mp_sub) (mp_int *a, mp_int *b, mp_int *c); /* 60 */
     int (*tclBN_mp_init_set_int) (mp_int*a, unsigned long i); /* 61 */
     int (*tclBN_mp_set_int) (mp_int*a, unsigned long i); /* 62 */
+    int (*tclBN_mp_cnt_lsb) (const mp_int*a); /* 63 */
 } TclTomMathStubs;
 
 #ifdef __cplusplus
@@ -482,6 +486,8 @@ extern const TclTomMathStubs *tclTomMathStubsPtr;
 	(tclTomMathStubsPtr->tclBN_mp_init_set_int) /* 61 */
 #define TclBN_mp_set_int \
 	(tclTomMathStubsPtr->tclBN_mp_set_int) /* 62 */
+#define TclBN_mp_cnt_lsb \
+	(tclTomMathStubsPtr->tclBN_mp_cnt_lsb) /* 63 */
 
 #endif /* defined(USE_TCL_STUBS) */
 
diff --git a/libtommath/bn_mp_cnt_lsb.c b/libtommath/bn_mp_cnt_lsb.c
index 6447a1f..f205e8c 100644
--- a/libtommath/bn_mp_cnt_lsb.c
+++ b/libtommath/bn_mp_cnt_lsb.c
@@ -20,7 +20,7 @@ static const int lnz[16] = {
 };
 
 /* Counts the number of lsbs which are zero before the first zero bit */
-int mp_cnt_lsb(mp_int *a)
+int mp_cnt_lsb(const mp_int *a)
 {
    int x;
    mp_digit q, qq;
diff --git a/macosx/Tcl.xcode/project.pbxproj b/macosx/Tcl.xcode/project.pbxproj
index 54d9e02..6801d54 100644
--- a/macosx/Tcl.xcode/project.pbxproj
+++ b/macosx/Tcl.xcode/project.pbxproj
@@ -104,6 +104,7 @@
 		F96D48ED08F272C3004A47F5 /* bn_mp_clear_multi.c in Sources */ = {isa = PBXBuildFile; fileRef = F96D426F08F272B3004A47F5 /* bn_mp_clear_multi.c */; };
 		F96D48EE08F272C3004A47F5 /* bn_mp_cmp.c in Sources */ = {isa = PBXBuildFile; fileRef = F96D427008F272B3004A47F5 /* bn_mp_cmp.c */; };
 		F96D48F008F272C3004A47F5 /* bn_mp_cmp_mag.c in Sources */ = {isa = PBXBuildFile; fileRef = F96D427208F272B3004A47F5 /* bn_mp_cmp_mag.c */; };
+		F96D48F208F272C3004A47F5 /* bn_mp_cnt_lsb.c in Sources */ = {isa = PBXBuildFile; fileRef = F96D427408F272B3004A47F5 /* bn_mp_cnt_lsb.c */; };
 		F96D48F208F272C3004A47F5 /* bn_mp_copy.c in Sources */ = {isa = PBXBuildFile; fileRef = F96D427408F272B3004A47F5 /* bn_mp_copy.c */; };
 		F96D48F308F272C3004A47F5 /* bn_mp_count_bits.c in Sources */ = {isa = PBXBuildFile; fileRef = F96D427508F272B3004A47F5 /* bn_mp_count_bits.c */; };
 		F96D48F408F272C3004A47F5 /* bn_mp_div.c in Sources */ = {isa = PBXBuildFile; fileRef = F96D427608F272B3004A47F5 /* bn_mp_div.c */; };
@@ -2068,6 +2069,7 @@
 				F96D48EE08F272C3004A47F5 /* bn_mp_cmp.c in Sources */,
 				F9E61D28090A481F002B3151 /* bn_mp_cmp_d.c in Sources */,
 				F96D48F008F272C3004A47F5 /* bn_mp_cmp_mag.c in Sources */,
+				F96D48F208F272C3004A47F5 /* bn_mp_cnt_lsb.c in Sources */,
 				F96D48F208F272C3004A47F5 /* bn_mp_copy.c in Sources */,
 				F96D48F308F272C3004A47F5 /* bn_mp_count_bits.c in Sources */,
 				F96D48F408F272C3004A47F5 /* bn_mp_div.c in Sources */,
diff --git a/macosx/Tcl.xcodeproj/project.pbxproj b/macosx/Tcl.xcodeproj/project.pbxproj
index 3cc34d7..b37f2e3 100644
--- a/macosx/Tcl.xcodeproj/project.pbxproj
+++ b/macosx/Tcl.xcodeproj/project.pbxproj
@@ -104,6 +104,7 @@
 		F96D48ED08F272C3004A47F5 /* bn_mp_clear_multi.c in Sources */ = {isa = PBXBuildFile; fileRef = F96D426F08F272B3004A47F5 /* bn_mp_clear_multi.c */; };
 		F96D48EE08F272C3004A47F5 /* bn_mp_cmp.c in Sources */ = {isa = PBXBuildFile; fileRef = F96D427008F272B3004A47F5 /* bn_mp_cmp.c */; };
 		F96D48F008F272C3004A47F5 /* bn_mp_cmp_mag.c in Sources */ = {isa = PBXBuildFile; fileRef = F96D427208F272B3004A47F5 /* bn_mp_cmp_mag.c */; };
+		F96D48F208F272C3004A47F5 /* bn_mp_cnt_lsb.c in Sources */ = {isa = PBXBuildFile; fileRef = F96D427408F272B3004A47F5 /* bn_mp_cnt_lsb.c */; };
 		F96D48F208F272C3004A47F5 /* bn_mp_copy.c in Sources */ = {isa = PBXBuildFile; fileRef = F96D427408F272B3004A47F5 /* bn_mp_copy.c */; };
 		F96D48F308F272C3004A47F5 /* bn_mp_count_bits.c in Sources */ = {isa = PBXBuildFile; fileRef = F96D427508F272B3004A47F5 /* bn_mp_count_bits.c */; };
 		F96D48F408F272C3004A47F5 /* bn_mp_div.c in Sources */ = {isa = PBXBuildFile; fileRef = F96D427608F272B3004A47F5 /* bn_mp_div.c */; };
@@ -2068,6 +2069,7 @@
 				F96D48EE08F272C3004A47F5 /* bn_mp_cmp.c in Sources */,
 				F9E61D28090A481F002B3151 /* bn_mp_cmp_d.c in Sources */,
 				F96D48F008F272C3004A47F5 /* bn_mp_cmp_mag.c in Sources */,
+				F96D48F208F272C3004A47F5 /* bn_mp_cnt_lsb.c in Sources */,
 				F96D48F208F272C3004A47F5 /* bn_mp_copy.c in Sources */,
 				F96D48F308F272C3004A47F5 /* bn_mp_count_bits.c in Sources */,
 				F96D48F408F272C3004A47F5 /* bn_mp_div.c in Sources */,
diff --git a/tests/util.test b/tests/util.test
index 1f214b5..d06925b 100644
--- a/tests/util.test
+++ b/tests/util.test
@@ -3858,6 +3858,35 @@ test util-16.1.17.307 {8.4 compatible formatting of doubles} \
     {expr 1e307} \
     9.9999999999999999e+306
 
+test util-17.1 {bankers' rounding [Bug 3349507]} {ieeeFloatingPoint} {
+    set r {}
+    foreach {input} {
+	0x1ffffffffffffc000
+	0x1ffffffffffffc800
+	0x1ffffffffffffd000
+	0x1ffffffffffffd800
+	0x1ffffffffffffe000
+	0x1ffffffffffffe800
+	0x1fffffffffffff000
+	0x1fffffffffffff800
+    } {
+	binary scan [binary format q [expr double($input)]] wu x
+	lappend r [format %#llx $x]
+	binary scan [binary format q [expr double(-$input)]] wu x
+	lappend r [format %#llx $x]
+    }
+    set r
+} [list {*}{
+    0x43fffffffffffffc 0xc3fffffffffffffc 
+    0x43fffffffffffffc 0xc3fffffffffffffc
+    0x43fffffffffffffd 0xc3fffffffffffffd
+    0x43fffffffffffffe 0xc3fffffffffffffe
+    0x43fffffffffffffe 0xc3fffffffffffffe
+    0x43fffffffffffffe 0xc3fffffffffffffe
+    0x43ffffffffffffff 0xc3ffffffffffffff
+    0x4400000000000000 0xc400000000000000
+}]
+
 set ::tcl_precision $saved_precision
 
 # cleanup
diff --git a/unix/Makefile.in b/unix/Makefile.in
index ce8e198..b3507ba 100644
--- a/unix/Makefile.in
+++ b/unix/Makefile.in
@@ -315,7 +315,8 @@ OO_OBJS = tclOO.o tclOOBasic.o tclOOCall.o tclOODefineCmds.o tclOOInfo.o \
 TOMMATH_OBJS = bncore.o bn_reverse.o bn_fast_s_mp_mul_digs.o \
 	bn_fast_s_mp_sqr.o bn_mp_add.o bn_mp_and.o \
         bn_mp_add_d.o bn_mp_clamp.o bn_mp_clear.o bn_mp_clear_multi.o \
-        bn_mp_cmp.o bn_mp_cmp_d.o bn_mp_cmp_mag.o bn_mp_copy.o \
+        bn_mp_cmp.o bn_mp_cmp_d.o bn_mp_cmp_mag.o \
+	bn_mp_cnt_lsb.o bn_mp_copy.o \
 	bn_mp_count_bits.o bn_mp_div.o bn_mp_div_d.o bn_mp_div_2.o \
 	bn_mp_div_2d.o bn_mp_div_3.o \
         bn_mp_exch.o bn_mp_expt_d.o bn_mp_grow.o bn_mp_init.o \
@@ -484,6 +485,7 @@ TOMMATH_SRCS = \
 	$(TOMMATH_DIR)/bn_mp_cmp_d.c \
 	$(TOMMATH_DIR)/bn_mp_cmp_mag.c \
 	$(TOMMATH_DIR)/bn_mp_copy.c \
+	$(TOMMATH_DIR)/bn_mp_cnt_lsb.c \
 	$(TOMMATH_DIR)/bn_mp_count_bits.c \
 	$(TOMMATH_DIR)/bn_mp_div.c \
 	$(TOMMATH_DIR)/bn_mp_div_d.c \
@@ -1337,6 +1339,9 @@ bn_mp_cmp_d.o: $(TOMMATH_DIR)/bn_mp_cmp_d.c $(MATHHDRS)
 bn_mp_cmp_mag.o: $(TOMMATH_DIR)/bn_mp_cmp_mag.c $(MATHHDRS)
 	$(CC) -c $(CC_SWITCHES) $(TOMMATH_DIR)/bn_mp_cmp_mag.c
 
+bn_mp_cnt_lsb.o: $(TOMMATH_DIR)/bn_mp_cnt_lsb.c $(MATHHDRS)
+	$(CC) -c $(CC_SWITCHES) $(TOMMATH_DIR)/bn_mp_cnt_lsb.c
+
 bn_mp_copy.o: $(TOMMATH_DIR)/bn_mp_copy.c $(MATHHDRS)
 	$(CC) -c $(CC_SWITCHES) $(TOMMATH_DIR)/bn_mp_copy.c
 
diff --git a/win/Makefile.in b/win/Makefile.in
index 7271c03..8a9359b 100644
--- a/win/Makefile.in
+++ b/win/Makefile.in
@@ -304,6 +304,7 @@ TOMMATH_OBJS = \
 	bn_mp_cmp.${OBJEXT} \
 	bn_mp_cmp_d.${OBJEXT} \
 	bn_mp_cmp_mag.${OBJEXT} \
+	bn_mp_cnt_lsb.${OBJEXT} \
 	bn_mp_copy.${OBJEXT} \
 	bn_mp_count_bits.${OBJEXT} \
 	bn_mp_div.${OBJEXT} \
diff --git a/win/makefile.vc b/win/makefile.vc
index 98d04a3..bc16584 100644
--- a/win/makefile.vc
+++ b/win/makefile.vc
@@ -362,6 +362,7 @@ TOMMATHOBJS = \
 	$(TMP_DIR)\bn_mp_cmp.obj \
 	$(TMP_DIR)\bn_mp_cmp_d.obj \
 	$(TMP_DIR)\bn_mp_cmp_mag.obj \
+	$(TMP_DIR)\bn_mp_cnt_lsb.obj \
 	$(TMP_DIR)\bn_mp_copy.obj \
 	$(TMP_DIR)\bn_mp_count_bits.obj \
 	$(TMP_DIR)\bn_mp_div.obj \