diff options
| -rw-r--r-- | ChangeLog | 8 | ||||
| -rw-r--r-- | libtommath/bn_mp_sqrt.c | 25 | ||||
| -rw-r--r-- | tests/expr.test | 23 |
3 files changed, 51 insertions, 5 deletions
@@ -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-02 Joe Mistachkin <joe@mistachkin.com> * tools/man2help2.tcl: Integrated patches from Harald Oehlmann. diff --git a/libtommath/bn_mp_sqrt.c b/libtommath/bn_mp_sqrt.c index 3aa0f38..aab357a 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.5.4.1 $ */ +/* $Date: 2008/10/05 21:25:23 $ */ diff --git a/tests/expr.test b/tests/expr.test index d0aaadb..30aebc6 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.72.2.1 2008/10/05 21:25:23 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 |