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$