diff options
Diffstat (limited to 'libtommath/bn_mp_sqrt.c')
-rw-r--r-- | libtommath/bn_mp_sqrt.c | 175 |
1 files changed, 89 insertions, 86 deletions
diff --git a/libtommath/bn_mp_sqrt.c b/libtommath/bn_mp_sqrt.c index 016b8ba..d342a32 100644 --- a/libtommath/bn_mp_sqrt.c +++ b/libtommath/bn_mp_sqrt.c @@ -1,5 +1,4 @@ -#include <tommath.h> - +#include "tommath_private.h" #ifdef BN_MP_SQRT_C /* LibTomMath, multiple-precision integer library -- Tom St Denis * @@ -12,8 +11,6 @@ * * The library is free for all purposes without any express * guarantee it works. - * - * Tom St Denis, tomstdenis@gmail.com, http://math.libtomcrypt.com */ #ifndef NO_FLOATING_POINT @@ -21,122 +18,128 @@ #endif /* this function is less generic than mp_n_root, simpler and faster */ -int mp_sqrt(mp_int *arg, mp_int *ret) +int mp_sqrt(const mp_int *arg, mp_int *ret) { - int res; - mp_int t1,t2; - int i, j, k; + int res; + mp_int t1, t2; + int i, j, k; #ifndef NO_FLOATING_POINT - volatile double d; - mp_digit dig; + volatile double d; + mp_digit dig; #endif - /* must be positive */ - if (arg->sign == MP_NEG) { - return MP_VAL; - } - - /* easy out */ - if (mp_iszero(arg) == MP_YES) { - mp_zero(ret); - return MP_OKAY; - } - - i = (arg->used / 2) - 1; - j = 2 * i; - if ((res = mp_init_size(&t1, i+2)) != MP_OKAY) { + /* must be positive */ + if (arg->sign == MP_NEG) { + return MP_VAL; + } + + /* easy out */ + if (mp_iszero(arg) == MP_YES) { + mp_zero(ret); + return MP_OKAY; + } + + 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; - } + } + + if ((res = mp_init(&t2)) != MP_OKAY) { + goto E2; + } - for (k = 0; k < i; ++k) { + 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. */ + /* 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]); - } + 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); + /* + * 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. + */ - /* dig is the most significant mp_digit of the square root */ + d = sqrt(d); - dig = (mp_digit) ldexp(d, -DIGIT_BIT); + /* dig is the most significant mp_digit of the square root */ - /* - * 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. - */ + dig = (mp_digit) ldexp(d, -DIGIT_BIT); - if (dig) { + /* + * 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; + 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; + t1.dp[i+1] = dig-1; + t1.dp[i] = MP_DIGIT_MAX; } - } else { + } 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. */ + /* 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; + 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) { - goto E1; - } - if ((res = mp_add(&t1,&t2,&t1)) != MP_OKAY) { - goto E1; - } - if ((res = mp_div_2(&t1,&t1)) != MP_OKAY) { - goto E1; - } - /* And now t1 > sqrt(arg) */ - do { - if ((res = mp_div(arg,&t1,&t2,NULL)) != MP_OKAY) { + /* t1 > 0 */ + if ((res = mp_div(arg, &t1, &t2, NULL)) != MP_OKAY) { goto E1; - } - if ((res = mp_add(&t1,&t2,&t1)) != MP_OKAY) { + } + if ((res = mp_add(&t1, &t2, &t1)) != MP_OKAY) { goto E1; - } - if ((res = mp_div_2(&t1,&t1)) != MP_OKAY) { + } + if ((res = mp_div_2(&t1, &t1)) != MP_OKAY) { goto E1; - } - /* t1 >= sqrt(arg) >= t2 at this point */ - } while (mp_cmp_mag(&t1,&t2) == MP_GT); + } + /* And now t1 > sqrt(arg) */ + do { + if ((res = mp_div(arg, &t1, &t2, NULL)) != MP_OKAY) { + goto E1; + } + if ((res = mp_add(&t1, &t2, &t1)) != MP_OKAY) { + goto E1; + } + if ((res = mp_div_2(&t1, &t1)) != MP_OKAY) { + goto E1; + } + /* t1 >= sqrt(arg) >= t2 at this point */ + } while (mp_cmp_mag(&t1, &t2) == MP_GT); - mp_exch(&t1,ret); + mp_exch(&t1, ret); -E1: mp_clear(&t2); -E2: mp_clear(&t1); - return res; +E1: + mp_clear(&t2); +E2: + mp_clear(&t1); + return res; } #endif + +/* ref: $Format:%D$ */ +/* git commit: $Format:%H$ */ +/* commit time: $Format:%ai$ */ |