faster-square-root patch

* mpi/mpi.c (s_highest_bit_mp, s_mp_set_bit): New functions.
(mp_sqrt): Rewrite with more efficient algorithm.

1 files changed, 70 insertions, 74 deletions

diff --git a/mpi/mpi.c b/mpi/mpi.c
index ed03059b..8ea4b1ba 100644
--- a/mpi/mpi.c
+++ b/mpi/mpi.c
@@ -156,6 +156,9 @@ static const char *s_dmap_2 =
 mp_err   s_mp_grow(mp_int *mp, mp_size min);   /* increase allocated size */
 mp_err   s_mp_pad(mp_int *mp, mp_size min);    /* left pad with zeroes    */
 
+int      s_highest_bit_mp(mp_int *a);
+mp_err   s_mp_set_bit(mp_int *a, int bit);
+
 void     s_mp_clamp(mp_int *mp);               /* clip leading zeroes     */
 
 void     s_mp_exch(mp_int *a, mp_int *b);      /* swap a and b in place   */
@@ -1533,77 +1536,61 @@ mp_err mp_mod_d(mp_int *a, mp_digit d, mp_digit *c)
 
 /* {{{ mp_sqrt(a, b) */
 
-/*
-  mp_sqrt(a, b)
-
-  Compute the integer square root of a, and store the result in b.
-  Uses an integer-arithmetic version of Newton's iterative linear
-  approximation technique to determine this value; the result has the
-  following two properties:
-
-     b^2 <= a
-     (b+1)^2 >= a
-
-  It is a range error to pass a negative value.
- */
 mp_err mp_sqrt(mp_int *a, mp_int *b)
 {
-  mp_int   x, t;
-  mp_err   res;
+  int mask_shift;
+  mp_int root, guess, *proot = &root, *pguess = &guess;
+  mp_int guess_sqr;
+  mp_err err = MP_MEM;
 
   ARGCHK(a != NULL && b != NULL, MP_BADARG);
 
-  /* Cannot take square root of a negative value */
-  if(SIGN(a) == MP_NEG)
+  if (mp_cmp_z(b) == MP_LT)
     return MP_RANGE;
 
-  /* Special cases for zero and one, trivial     */
-  if(mp_cmp_d(a, 0) == MP_EQ || mp_cmp_d(a, 1) == MP_EQ) 
-    return mp_copy(a, b);
-    
-  /* Initialize the temporaries we'll use below  */
-  if((res = mp_init_size(&t, USED(a))) != MP_OKAY)
-    return res;
-
-  /* Compute an initial guess for the iteration as a itself */
-  if((res = mp_init_copy(&x, a)) != MP_OKAY)
-    goto X;
-
-  for(;;) {
-    /* t = (x * x) - a */
-    mp_copy(&x, &t);      /* can't fail, t is big enough for original x */
-    if((res = mp_sqr(&t, &t)) != MP_OKAY ||
-       (res = mp_sub(&t, a, &t)) != MP_OKAY)
-      goto CLEANUP;
-
-    /* t = t / 2x       */
-    s_mp_mul_2(&x);
-    if((res = mp_div(&t, &x, &t, NULL)) != MP_OKAY)
-      goto CLEANUP;
-    s_mp_div_2(&x);
-
-    /* Terminate the loop, if the quotient is zero */
-    if(mp_cmp_z(&t) == MP_EQ)
-      break;
-
-    /* x = x - t       */
-    if((res = mp_sub(&x, &t, &x)) != MP_OKAY)
-      goto CLEANUP;
-
+  if ((err = mp_init(&root)))
+    goto out;
+  if ((err = mp_init(&guess)))
+    goto cleanup_root;
+  if ((err = mp_init(&guess_sqr)))
+    goto cleanup_guess;
+
+  for (mask_shift = s_highest_bit_mp(a) / 2; mask_shift >= 0; mask_shift--) {
+    mp_int *temp;
+    int cmp;
+
+    if ((err = mp_copy(proot, pguess)))
+	goto cleanup;
+    if ((err = s_mp_set_bit(pguess, mask_shift)))
+	goto cleanup;
+    if ((err = mp_copy(pguess, &guess_sqr)))
+	goto cleanup;
+    if ((err = s_mp_sqr(&guess_sqr)))
+	goto cleanup;
+
+    cmp = s_mp_cmp(&guess_sqr, a);
+
+    if (cmp < 0) {
+	temp = proot;
+	proot = pguess;
+	pguess = temp;
+    } else if (cmp == 0) {
+	proot = pguess;
+	break;
+    }
   }
 
-  /* Copy result to output parameter */
-  mp_sub_d(&x, 1, &x);
-  s_mp_exch(&x, b);
-
- CLEANUP:
-  mp_clear(&x);
- X:
-  mp_clear(&t); 
+  err = mp_copy(proot, b);
 
-  return res;
-
-} /* end mp_sqrt() */
+cleanup:
+  mp_clear(&guess_sqr);
+cleanup_guess:
+  mp_clear(&guess);
+cleanup_root:
+  mp_clear(&root);
+out:
+  return err;
+}
 
 /* }}} */
 
@@ -2552,21 +2539,9 @@ mp_err mp_to_unsigned_bin(mp_int *mp, unsigned char *str)
 
 int    mp_count_bits(mp_int *mp)
 {
-  int      len;
-  mp_digit d;
-
   ARGCHK(mp != NULL, MP_BADARG);
 
-  len = DIGIT_BIT * (USED(mp) - 1);
-  d = DIGIT(mp, USED(mp) - 1);
-
-  while(d != 0) {
-    ++len;
-    d >>= 1;
-  }
-
-  return len;
-  
+  return s_highest_bit_mp(mp);
 } /* end mp_count_bits() */
 
 /* }}} */
@@ -3130,6 +3105,27 @@ static int s_highest_bit(mp_digit n)
   abort();
 }
 
+int s_highest_bit_mp(mp_int *a)
+{
+  int nd1 = USED(a) - 1;
+  return s_highest_bit(DIGIT(a, nd1)) + nd1 * MP_DIGIT_BIT;
+}
+
+mp_err s_mp_set_bit(mp_int *a, int bit)
+{
+  int nd = (bit + MP_DIGIT_BIT) / MP_DIGIT_BIT;
+  int nbit = bit - (nd - 1) * MP_DIGIT_BIT;
+  mp_err res;
+
+  if (nd == 0)
+    return MP_OKAY;
+
+  if ((res = s_mp_pad(a, nd)) != MP_OKAY)
+    return res;
+
+  DIGIT(a, nd - 1) |= ((mp_digit) 1 << nbit);
+  return MP_OKAY;
+}
 
 /* {{{ s_mp_exch(a, b) */