Corrected bugs in tommath installation, improved tommath square root performance, patched around a [clock scan] issue with time zones

1 files changed, 50 insertions, 8 deletions

diff --git a/libtommath/bn_mp_sqrt.c b/libtommath/bn_mp_sqrt.c
index ac8c2b8..72a9ff5 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
  *
@@ -15,11 +16,20 @@
  * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.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
+  double d;
+  mp_digit dig;
+#endif
 
   /* must be positive */
   if (arg->sign == MP_NEG) {
@@ -31,17 +41,49 @@ 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]);
+  }
+  d = sqrt(d);
+  dig = (mp_digit) ldexp(d, -DIGIT_BIT);
+  if (dig) {
+      t1.used = i+2;
+      t1.dp[i+1] = dig;
+      d -= ldexp((double) dig, DIGIT_BIT);
+  } else {
+      t1.used = i+1;
+  }
+  t1.dp[i] = (mp_digit) d;
+
+#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) {
@@ -77,5 +119,5 @@ E2: mp_clear(&t1);
 #endif
 
 /* $Source: /root/tcl/repos-to-convert/tcl/libtommath/bn_mp_sqrt.c,v $ */
-/* $Revision: 1.1.1.2 $ */
-/* $Date: 2005/09/26 16:31:56 $ */
+/* $Revision: 1.2 $ */
+/* $Date: 2005/12/27 17:39:02 $ */