diff options
Diffstat (limited to 'libtommath/bn_mp_sqrt.c')
| -rw-r--r-- | libtommath/bn_mp_sqrt.c | 83 |
1 files changed, 72 insertions, 11 deletions
diff --git a/libtommath/bn_mp_sqrt.c b/libtommath/bn_mp_sqrt.c index 8fd057c..016b8ba 100644 --- a/libtommath/bn_mp_sqrt.c +++ b/libtommath/bn_mp_sqrt.c @@ -1,4 +1,5 @@ #include <tommath.h> + #ifdef BN_MP_SQRT_C /* LibTomMath, multiple-precision integer library -- Tom St Denis * @@ -12,14 +13,23 @@ * The library is free for all purposes without any express * guarantee it works. * - * Tom St Denis, tomstdenis@gmail.com, http://libtom.org + * Tom St Denis, tomstdenis@gmail.com, http://math.libtomcrypt.com */ +#ifndef NO_FLOATING_POINT +#include <math.h> +#endif + /* this function is less generic than mp_n_root, simpler and faster */ int mp_sqrt(mp_int *arg, mp_int *ret) { int res; mp_int t1,t2; + int i, j, k; +#ifndef NO_FLOATING_POINT + volatile double d; + mp_digit dig; +#endif /* must be positive */ if (arg->sign == MP_NEG) { @@ -31,17 +41,72 @@ int mp_sqrt(mp_int *arg, mp_int *ret) mp_zero(ret); return MP_OKAY; } - - if ((res = mp_init_copy(&t1, arg)) != MP_OKAY) { - return res; + + i = (arg->used / 2) - 1; + j = 2 * i; + if ((res = mp_init_size(&t1, i+2)) != MP_OKAY) { + return res; } - + if ((res = mp_init(&t2)) != MP_OKAY) { goto E2; } - /* First approx. (not very bad for large arg) */ - mp_rshd (&t1,t1.used/2); + for (k = 0; k < i; ++k) { + t1.dp[k] = (mp_digit) 0; + } + +#ifndef NO_FLOATING_POINT + + /* Estimate the square root using the hardware floating point unit. */ + + d = 0.0; + 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 >= 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; + } + +#else + + /* Estimate the square root as having 1 in the most significant place. */ + + t1.used = i + 2; + t1.dp[i+1] = (mp_digit) 1; + t1.dp[i] = (mp_digit) 0; + +#endif /* t1 > 0 */ if ((res = mp_div(arg,&t1,&t2,NULL)) != MP_OKAY) { @@ -75,7 +140,3 @@ E2: mp_clear(&t1); } #endif - -/* $Source: /cvs/libtom/libtommath/bn_mp_sqrt.c,v $ */ -/* $Revision: 1.4 $ */ -/* $Date: 2006/12/28 01:25:13 $ */ |