diff options
Diffstat (limited to 'libtommath')
-rw-r--r-- | libtommath/bn_mp_cmp.c | 2 | ||||
-rw-r--r-- | libtommath/bn_mp_cmp_d.c | 2 | ||||
-rw-r--r-- | libtommath/bn_mp_cmp_mag.c | 2 | ||||
-rw-r--r-- | libtommath/bn_mp_cnt_lsb.c | 2 | ||||
-rw-r--r-- | libtommath/bn_mp_copy.c | 2 | ||||
-rw-r--r-- | libtommath/bn_mp_count_bits.c | 2 | ||||
-rw-r--r-- | libtommath/bn_mp_div_2d.c | 2 | ||||
-rw-r--r-- | libtommath/bn_mp_div_d.c | 7 | ||||
-rw-r--r-- | libtommath/bn_mp_init_copy.c | 2 | ||||
-rw-r--r-- | libtommath/bn_mp_mod_2d.c | 2 | ||||
-rw-r--r-- | libtommath/bn_mp_mul_2d.c | 2 | ||||
-rw-r--r-- | libtommath/bn_mp_neg.c | 2 | ||||
-rw-r--r-- | libtommath/bn_mp_radix_size.c | 2 | ||||
-rw-r--r-- | libtommath/bn_mp_sqrt.c | 76 | ||||
-rw-r--r-- | libtommath/tommath.h | 28 | ||||
-rw-r--r-- | libtommath/tommath_private.h | 2 | ||||
-rwxr-xr-x | libtommath/updatemakes.sh | 33 |
17 files changed, 101 insertions, 69 deletions
diff --git a/libtommath/bn_mp_cmp.c b/libtommath/bn_mp_cmp.c index e757ddf..09cafa0 100644 --- a/libtommath/bn_mp_cmp.c +++ b/libtommath/bn_mp_cmp.c @@ -17,7 +17,7 @@ /* compare two ints (signed)*/ int -mp_cmp (mp_int * a, mp_int * b) +mp_cmp (const mp_int * a, const mp_int * b) { /* compare based on sign */ if (a->sign != b->sign) { diff --git a/libtommath/bn_mp_cmp_d.c b/libtommath/bn_mp_cmp_d.c index 3f5ebae..3768f06 100644 --- a/libtommath/bn_mp_cmp_d.c +++ b/libtommath/bn_mp_cmp_d.c @@ -16,7 +16,7 @@ */ /* compare a digit */ -int mp_cmp_d(mp_int * a, mp_digit b) +int mp_cmp_d(const mp_int * a, mp_digit b) { /* compare based on sign */ if (a->sign == MP_NEG) { diff --git a/libtommath/bn_mp_cmp_mag.c b/libtommath/bn_mp_cmp_mag.c index 7ceda97..3b7f1b4 100644 --- a/libtommath/bn_mp_cmp_mag.c +++ b/libtommath/bn_mp_cmp_mag.c @@ -16,7 +16,7 @@ */ /* compare maginitude of two ints (unsigned) */ -int mp_cmp_mag (mp_int * a, mp_int * b) +int mp_cmp_mag (const mp_int * a, const mp_int * b) { int n; mp_digit *tmpa, *tmpb; diff --git a/libtommath/bn_mp_cnt_lsb.c b/libtommath/bn_mp_cnt_lsb.c index bf201b5..10a901f 100644 --- a/libtommath/bn_mp_cnt_lsb.c +++ b/libtommath/bn_mp_cnt_lsb.c @@ -20,7 +20,7 @@ static const int lnz[16] = { }; /* Counts the number of lsbs which are zero before the first zero bit */ -int mp_cnt_lsb(mp_int *a) +int mp_cnt_lsb(const mp_int *a) { int x; mp_digit q, qq; diff --git a/libtommath/bn_mp_copy.c b/libtommath/bn_mp_copy.c index 84e839e..435729c 100644 --- a/libtommath/bn_mp_copy.c +++ b/libtommath/bn_mp_copy.c @@ -17,7 +17,7 @@ /* copy, b = a */ int -mp_copy (mp_int * a, mp_int * b) +mp_copy (const mp_int * a, mp_int * b) { int res, n; diff --git a/libtommath/bn_mp_count_bits.c b/libtommath/bn_mp_count_bits.c index ff558eb..10ab695 100644 --- a/libtommath/bn_mp_count_bits.c +++ b/libtommath/bn_mp_count_bits.c @@ -17,7 +17,7 @@ /* returns the number of bits in an int */ int -mp_count_bits (mp_int * a) +mp_count_bits (const mp_int * a) { int r; mp_digit q; diff --git a/libtommath/bn_mp_div_2d.c b/libtommath/bn_mp_div_2d.c index 635d374..d76de1a 100644 --- a/libtommath/bn_mp_div_2d.c +++ b/libtommath/bn_mp_div_2d.c @@ -16,7 +16,7 @@ */ /* shift right by a certain bit count (store quotient in c, optional remainder in d) */ -int mp_div_2d (mp_int * a, int b, mp_int * c, mp_int * d) +int mp_div_2d (const mp_int * a, int b, mp_int * c, mp_int * d) { mp_digit D, r, rr; int x, res; diff --git a/libtommath/bn_mp_div_d.c b/libtommath/bn_mp_div_d.c index a5dbc59..d1f2774 100644 --- a/libtommath/bn_mp_div_d.c +++ b/libtommath/bn_mp_div_d.c @@ -19,12 +19,11 @@ static int s_is_power_of_two(mp_digit b, int *p) { int x; - /* fast return if no power of two */ + /* quick out - if (b & (b-1)) isn't zero, b isn't a power of two */ if ((b == 0) || ((b & (b-1)) != 0)) { - return 0; + return 0; } - - for (x = 0; x < DIGIT_BIT; x++) { + for (x = 1; x < DIGIT_BIT; x++) { if (b == (((mp_digit)1)<<x)) { *p = x; return 1; diff --git a/libtommath/bn_mp_init_copy.c b/libtommath/bn_mp_init_copy.c index 37a57ec..2a06e08 100644 --- a/libtommath/bn_mp_init_copy.c +++ b/libtommath/bn_mp_init_copy.c @@ -16,7 +16,7 @@ */ /* creates "a" then copies b into it */ -int mp_init_copy (mp_int * a, mp_int * b) +int mp_init_copy (mp_int * a, const mp_int * b) { int res; diff --git a/libtommath/bn_mp_mod_2d.c b/libtommath/bn_mp_mod_2d.c index 2bb86da..ed427fd 100644 --- a/libtommath/bn_mp_mod_2d.c +++ b/libtommath/bn_mp_mod_2d.c @@ -17,7 +17,7 @@ /* calc a value mod 2**b */ int -mp_mod_2d (mp_int * a, int b, mp_int * c) +mp_mod_2d (const mp_int * a, int b, mp_int * c) { int x, res; diff --git a/libtommath/bn_mp_mul_2d.c b/libtommath/bn_mp_mul_2d.c index c00fd7e..e167fc1 100644 --- a/libtommath/bn_mp_mul_2d.c +++ b/libtommath/bn_mp_mul_2d.c @@ -16,7 +16,7 @@ */ /* shift left by a certain bit count */ -int mp_mul_2d (mp_int * a, int b, mp_int * c) +int mp_mul_2d (const mp_int * a, int b, mp_int * c) { mp_digit d; int res; diff --git a/libtommath/bn_mp_neg.c b/libtommath/bn_mp_neg.c index d03e92e..da50e28 100644 --- a/libtommath/bn_mp_neg.c +++ b/libtommath/bn_mp_neg.c @@ -16,7 +16,7 @@ */ /* b = -a */ -int mp_neg (mp_int * a, mp_int * b) +int mp_neg (const mp_int * a, mp_int * b) { int res; if (a != b) { diff --git a/libtommath/bn_mp_radix_size.c b/libtommath/bn_mp_radix_size.c index cb0c134..9c502c8 100644 --- a/libtommath/bn_mp_radix_size.c +++ b/libtommath/bn_mp_radix_size.c @@ -16,7 +16,7 @@ */ /* returns size of ASCII reprensentation */ -int mp_radix_size (mp_int * a, int radix, int *size) +int mp_radix_size (const mp_int * a, int radix, int *size) { int res, digs; mp_int t; diff --git a/libtommath/bn_mp_sqrt.c b/libtommath/bn_mp_sqrt.c index d3b5d62..775e85f 100644 --- a/libtommath/bn_mp_sqrt.c +++ b/libtommath/bn_mp_sqrt.c @@ -15,11 +15,20 @@ * Tom St Denis, tstdenis82@gmail.com, http://libtom.org */ +#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 +40,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) { diff --git a/libtommath/tommath.h b/libtommath/tommath.h index 80bae69..7d92d97 100644 --- a/libtommath/tommath.h +++ b/libtommath/tommath.h @@ -238,10 +238,10 @@ int mp_init_set (mp_int * a, mp_digit b); int mp_init_set_int (mp_int * a, unsigned long b); /* copy, b = a */ -int mp_copy(mp_int *a, mp_int *b); +int mp_copy(const mp_int *a, mp_int *b); /* inits and copies, a = b */ -int mp_init_copy(mp_int *a, mp_int *b); +int mp_init_copy(mp_int *a, const mp_int *b); /* trim unused digits */ void mp_clamp(mp_int *a); @@ -261,25 +261,25 @@ void mp_rshd(mp_int *a, int b); int mp_lshd(mp_int *a, int b); /* c = a / 2**b, implemented as c = a >> b */ -int mp_div_2d(mp_int *a, int b, mp_int *c, mp_int *d); +int mp_div_2d(const mp_int *a, int b, mp_int *c, mp_int *d); /* b = a/2 */ int mp_div_2(mp_int *a, mp_int *b); /* c = a * 2**b, implemented as c = a << b */ -int mp_mul_2d(mp_int *a, int b, mp_int *c); +int mp_mul_2d(const mp_int *a, int b, mp_int *c); /* b = a*2 */ int mp_mul_2(mp_int *a, mp_int *b); /* c = a mod 2**b */ -int mp_mod_2d(mp_int *a, int b, mp_int *c); +int mp_mod_2d(const mp_int *a, int b, mp_int *c); /* computes a = 2**b */ int mp_2expt(mp_int *a, int b); /* Counts the number of lsbs which are zero before the first zero bit */ -int mp_cnt_lsb(mp_int *a); +int mp_cnt_lsb(const mp_int *a); /* I Love Earth! */ @@ -299,16 +299,16 @@ int mp_and(mp_int *a, mp_int *b, mp_int *c); /* ---> Basic arithmetic <--- */ /* b = -a */ -int mp_neg(mp_int *a, mp_int *b); +int mp_neg(const mp_int *a, mp_int *b); /* b = |a| */ int mp_abs(mp_int *a, mp_int *b); /* compare a to b */ -int mp_cmp(mp_int *a, mp_int *b); +int mp_cmp(const mp_int *a, const mp_int *b); /* compare |a| to |b| */ -int mp_cmp_mag(mp_int *a, mp_int *b); +int mp_cmp_mag(const mp_int *a, const mp_int *b); /* c = a + b */ int mp_add(mp_int *a, mp_int *b, mp_int *c); @@ -331,7 +331,7 @@ int mp_mod(mp_int *a, mp_int *b, mp_int *c); /* ---> single digit functions <--- */ /* compare against a single digit */ -int mp_cmp_d(mp_int *a, mp_digit b); +int mp_cmp_d(const mp_int *a, mp_digit b); /* c = a + b */ int mp_add_d(mp_int *a, mp_digit b, mp_int *c); @@ -455,9 +455,9 @@ int mp_exptmod(mp_int *a, mp_int *b, mp_int *c, mp_int *d); /* number of primes */ #ifdef MP_8BIT - #define PRIME_SIZE 31 +# define PRIME_SIZE 31 #else - #define PRIME_SIZE 256 +# define PRIME_SIZE 256 #endif /* table of first PRIME_SIZE primes */ @@ -524,7 +524,7 @@ int mp_prime_next_prime(mp_int *a, int t, int bbs_style); int mp_prime_random_ex(mp_int *a, int t, int size, int flags, ltm_prime_callback cb, void *dat); /* ---> radix conversion <--- */ -int mp_count_bits(mp_int *a); +int mp_count_bits(const mp_int *a); int mp_unsigned_bin_size(mp_int *a); int mp_read_unsigned_bin(mp_int *a, const unsigned char *b, int c); @@ -539,7 +539,7 @@ int mp_to_signed_bin_n (mp_int * a, unsigned char *b, unsigned long *outlen); int mp_read_radix(mp_int *a, const char *str, int radix); int mp_toradix(mp_int *a, char *str, int radix); int mp_toradix_n(mp_int * a, char *str, int radix, int maxlen); -int mp_radix_size(mp_int *a, int radix, int *size); +int mp_radix_size(const mp_int *a, int radix, int *size); #ifndef LTM_NO_FILE int mp_fread(mp_int *a, int radix, FILE *stream); diff --git a/libtommath/tommath_private.h b/libtommath/tommath_private.h index aeda684..3217ee7 100644 --- a/libtommath/tommath_private.h +++ b/libtommath/tommath_private.h @@ -40,6 +40,7 @@ extern "C" { #endif /* define heap macros */ +#if 0 #ifndef XMALLOC /* default to libc stuff */ #define XMALLOC malloc @@ -53,6 +54,7 @@ extern "C" { extern void *XCALLOC(size_t n, size_t s); extern void XFREE(void *p); #endif +#endif /* lowlevel functions, do not call! */ int s_mp_add(mp_int *a, mp_int *b, mp_int *c); diff --git a/libtommath/updatemakes.sh b/libtommath/updatemakes.sh deleted file mode 100755 index 0f9520e..0000000 --- a/libtommath/updatemakes.sh +++ /dev/null @@ -1,33 +0,0 @@ -#!/bin/bash - -bash genlist.sh > tmplist - -perl filter.pl makefile tmplist -sed -e 's/ *$//' < tmp.delme > makefile -rm -f tmp.delme - -perl filter.pl makefile.icc tmplist -sed -e 's/ *$//' < tmp.delme > makefile.icc -rm -f tmp.delme - -perl filter.pl makefile.shared tmplist -sed -e 's/ *$//' < tmp.delme > makefile.shared -rm -f tmp.delme - -perl filter.pl makefile.cygwin_dll tmplist -sed -e 's/ *$//' < tmp.delme > makefile.cygwin_dll -rm -f tmp.delme - -perl filter.pl makefile.bcc tmplist -sed -e 's/\.o /.obj /g' -e 's/ *$//' < tmp.delme > makefile.bcc -rm -f tmp.delme - -perl filter.pl makefile.msvc tmplist -sed -e 's/\.o /.obj /g' -e 's/ *$//' < tmp.delme > makefile.msvc -rm -f tmp.delme - -rm -f tmplist - -# ref: $Format:%D$ -# git commit: $Format:%H$ -# commit time: $Format:%ai$ |