diff options
Diffstat (limited to 'libtommath/bn_mp_sqrt.c')
| -rw-r--r-- | libtommath/bn_mp_sqrt.c | 35 |
1 files changed, 27 insertions, 8 deletions
diff --git a/libtommath/bn_mp_sqrt.c b/libtommath/bn_mp_sqrt.c index 72a9ff5..016b8ba 100644 --- a/libtommath/bn_mp_sqrt.c +++ b/libtommath/bn_mp_sqrt.c @@ -13,7 +13,7 @@ * The library is free for all purposes without any express * guarantee it works. * - * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org + * Tom St Denis, tomstdenis@gmail.com, http://math.libtomcrypt.com */ #ifndef NO_FLOATING_POINT @@ -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,16 +64,39 @@ 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; - t1.dp[i+1] = dig; d -= ldexp((double) dig, DIGIT_BIT); + if (d >= 1.0) { + t1.dp[i+1] = dig; + t1.dp[i] = ((mp_digit) d) - 1; + } else { + t1.dp[i+1] = dig-1; + t1.dp[i] = MP_DIGIT_MAX; + } } else { t1.used = i+1; + t1.dp[i] = ((mp_digit) d) - 1; } - t1.dp[i] = (mp_digit) d; #else @@ -117,7 +140,3 @@ E2: mp_clear(&t1); } #endif - -/* $Source: /root/tcl/repos-to-convert/tcl/libtommath/bn_mp_sqrt.c,v $ */ -/* $Revision: 1.2 $ */ -/* $Date: 2005/12/27 17:39:02 $ */ |