* libtommath/bn_mp_sqrt.c (bn_mp_sqrt): Handle the case where

* tests/expr.test (expr-47.13):         a number's square root
is between n<<DIGIT_BIT and n<<DIGIT_BIT+1. [Bug 2143288]
Thanks to Malcolm Boffey (malcolm.boffey@virgin.net) for
the patch.

3 files changed, 51 insertions, 5 deletions

diff --git a/ChangeLog b/ChangeLog
index 6191d35..324dfba 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,11 @@
+2008-10-05  Kevin B, Kenny  <kennykb@acm.org>
+
+	* libtommath/bn_mp_sqrt.c (bn_mp_sqrt): Handle the case where
+	* tests/expr.test (expr-47.13):         a number's square root
+	is between n<<DIGIT_BIT and n<<DIGIT_BIT+1. [Bug 2143288]
+	Thanks to Malcolm Boffey (malcolm.boffey@virgin.net) for
+	the patch.
+	
 2008-10-04  Jan Nijtmans  <nijtmans@users.sf.net>
 
 	* generic/tclInt.decls:  CONSTified the AuxDataType argument
diff --git a/libtommath/bn_mp_sqrt.c b/libtommath/bn_mp_sqrt.c
index 3aa0f38..0bba337 100644
--- a/libtommath/bn_mp_sqrt.c
+++ b/libtommath/bn_mp_sqrt.c
@@ -27,7 +27,7 @@ int mp_sqrt(mp_int *arg, mp_int *ret)
   mp_int t1,t2;
   int i, j, k;
 #ifndef NO_FLOATING_POINT
-  double d;
+  volatile double d;
   mp_digit dig;
 #endif
 
@@ -64,12 +64,29 @@ int mp_sqrt(mp_int *arg, mp_int *ret)
   for (k = arg->used-1; k >= j; --k) {
       d = ldexp(d, DIGIT_BIT) + (double) (arg->dp[k]);
   }
+
+  /* 
+   * At this point, d is the nearest floating point number to the most
+   * significant 1 or 2 mp_digits of arg. Extract its square root.
+   */
+     
   d = sqrt(d);
+
+  /* dig is the most significant mp_digit of the square root */
+
   dig = (mp_digit) ldexp(d, -DIGIT_BIT);
+
+  /* 
+   * If the most significant digit is nonzero, find the next digit down
+   * by subtracting DIGIT_BIT times thie most significant digit. 
+   * Subtract one from the result so that our initial estimate is always
+   * low.
+   */
+
   if (dig) {
       t1.used = i+2;
       d -= ldexp((double) dig, DIGIT_BIT);
-      if (d != 0.0) {
+      if (d >= 1.0) {
 	  t1.dp[i+1] = dig;
 	  t1.dp[i] = ((mp_digit) d) - 1;
       } else {
@@ -126,5 +143,5 @@ E2: mp_clear(&t1);
 
 /* $Source: /root/tcl/repos-to-convert/tcl/libtommath/bn_mp_sqrt.c,v $ */
 /* Based on Tom's 1.3 */
-/* $Revision: 1.5 $ */
-/* $Date: 2006/12/01 05:48:23 $ */
+/* $Revision: 1.6 $ */
+/* $Date: 2008/10/05 21:27:07 $ */
diff --git a/tests/expr.test b/tests/expr.test
index d0aaadb..e51e4c1 100644
--- a/tests/expr.test
+++ b/tests/expr.test
@@ -10,7 +10,7 @@
 # See the file "license.terms" for information on usage and redistribution
 # of this file, and for a DISCLAIMER OF ALL WARRANTIES.
 #
-# RCS: @(#) $Id: expr.test,v 1.72 2008/03/10 16:18:55 dgp Exp $
+# RCS: @(#) $Id: expr.test,v 1.73 2008/10/05 21:27:07 kennykb Exp $
 
 if {[lsearch [namespace children] ::tcltest] == -1} {
     package require tcltest 2.1
@@ -6770,6 +6770,7 @@ test expr-47.11 {isqrt of zero} {
 
 test expr-47.12 {isqrt of various sizes of integer} {
     set faults 0
+    set trouble {}
     for {set i 0} {$faults < 10 && $i <= 1024} {incr i} {
 	set root [expr {1 << $i}]
 	set rm1 [expr {$root - 1}]
@@ -6796,6 +6797,26 @@ test expr-47.12 {isqrt of various sizes of integer} {
     set trouble
 } {}
 
+test expr-47.13 {isqrt and floating point rounding (Bug 2143288)} {
+    set trouble {}
+    set faults 0
+    for {set i 0} {$i < 29 && $faults < 10} {incr i} {
+	for {set j 0} {$j <= $i} {incr j} {
+	    set k [expr {isqrt((1<<56)+(1<<$i)+(1<<$j))}]
+	    if {$k != (1<<28)} {
+		append trouble "i = $i, j = $j, k = $k\n"
+		incr faults
+	    }
+	}
+	set k [expr {isqrt((1<<56)+(1<<29)+(1<<$i))}]
+	if {$k != (1<<28)+1} {
+	    append trouble "i = $i, k = $k\n"
+	    incr faults
+	}
+    }
+    set trouble
+} {}
+
 test expr-48.1 {Bug 1770224} {
     expr {-0x8000000000000001 >> 0x8000000000000000}
 } -1