Use sliding-window exponentiation.

diff --git a/mpbarrett.c b/mpbarrett.c
index 0af753e..889b268 100644 (file)
--- a/mpbarrett.c
+++ b/mpbarrett.c
@@ -1,6 +1,6 @@
 /* -*-c-*-
  *
- * $Id: mpbarrett.c,v 1.6 2000/10/08 12:03:44 mdw Exp $
+ * $Id: mpbarrett.c,v 1.7 2001/04/19 18:25:26 mdw Exp $
  *
  * Barrett modular reduction
  *
@@ -30,6 +30,9 @@
 /*----- Revision history --------------------------------------------------* 
  *
  * $Log: mpbarrett.c,v $
+ * Revision 1.7  2001/04/19 18:25:26  mdw
+ * Use sliding-window exponentiation.
+ *
  * Revision 1.6  2000/10/08 12:03:44  mdw
  * (mpbarrett_reduce): Cope with negative numbers.
  *
@@ -196,7 +199,12 @@ mp *mpbarrett_reduce(mpbarrett *mb, mp *d, mp *m)
  * Returns:     Result, %$a^e \bmod m$%.
  */
 
-mp *mpbarrett_exp(mpbarrett *mb, mp *d, mp *a, mp *e)
+#define WINSZ 5
+#define TABSZ (1 << (WINSZ - 1))
+
+#define THRESH (((MPW_BITS / WINSZ) << 2) + 1)
+
+static mp *exp_simple(mpbarrett *mb, mp *d, mp *a, mp *e)
 {
   mpscan sc;
   mp *x = MP_ONE;
@@ -255,6 +263,130 @@ exit:
   return (x);
 }
 
+mp *mpbarrett_exp(mpbarrett *mb, mp *d, mp *a, mp *e)
+{
+  mp **tab;
+  mp *a2;
+  mp *spare = (e->f & MP_BURN) ? MP_NEWSEC : MP_NEW;
+  mp *x = MP_ONE;
+  unsigned i, sq = 0;
+  mpscan sc;
+
+  /* --- Do we bother? --- */
+
+  MP_SHRINK(e);
+  if (MP_LEN(e) == 0)
+    goto exit;
+  if (MP_LEN(e) < THRESH) {
+    x->ref--;
+    return (exp_simple(mb, d, a, e));
+  }
+
+  /* --- Do the precomputation --- */
+
+  a2 = mp_sqr(MP_NEW, a);
+  a2 = mpbarrett_reduce(mb, a2, a2);
+  tab = xmalloc(TABSZ * sizeof(mp *));
+  tab[0] = MP_COPY(a);
+  for (i = 1; i < TABSZ; i++) {
+    mp *x = mp_mul(MP_NEW, tab[i - 1], a2);
+    tab[i] = mpbarrett_reduce(mb, x, x);
+  }
+  mp_drop(a2);
+  mp_rscan(&sc, e);
+  
+  /* --- Skip top-end zero bits --- *
+   *
+   * If the initial step worked, there must be a set bit somewhere, so keep
+   * stepping until I find it.
+   */
+
+  MP_RSTEP(&sc);
+  while (!MP_RBIT(&sc))
+    MP_RSTEP(&sc);
+
+  /* --- Now for the main work --- */
+
+  for (;;) {
+    unsigned l = 0;
+    unsigned z = 0;
+
+    /* --- The next bit is set, so read a window index --- *
+     *
+     * Reset @i@ to zero and increment @sq@.  Then, until either I read
+     * @WINSZ@ bits or I run out of bits, scan in a bit: if it's clear, bump
+     * the @z@ counter; if it's set, push a set bit into @i@, shift it over
+     * by @z@ bits, bump @sq@ by @z + 1@ and clear @z@.  By the end of this
+     * palaver, @i@ is an index to the precomputed value in @tab@.
+     */
+
+    i = 0;
+    sq++;
+    for (;;) {
+      l++;
+      if (l >= WINSZ || !MP_RSTEP(&sc))
+       break;
+      if (!MP_RBIT(&sc))
+       z++;
+      else {
+       i = ((i << 1) | 1) << z;
+       sq += z + 1;
+       z = 0;
+      }
+    }
+
+    /* --- Do the squaring --- *
+     *
+     * Remember that @sq@ carries over from the zero-skipping stuff below.
+     */
+
+    while (sq) {
+      mp *y;
+      y = mp_sqr(spare, x);
+      y = mpbarrett_reduce(mb, y, y);
+      spare = x; x = y;
+      sq--;
+    }
+
+    /* --- Do the multiply --- */
+
+    { mp *y = mp_mul(spare, x, tab[i]); spare = x;
+      x = mpbarrett_reduce(mb, y, y); }
+
+    /* --- Now grind along through the rest of the bits --- */
+
+    sq = z;
+    for (;;) {
+      if (!MP_RSTEP(&sc))
+       goto done;
+      if (MP_RBIT(&sc))
+       break;
+      sq++;
+    }
+  }
+
+  /* --- Do a final round of squaring --- */
+
+done:
+  while (sq) {
+    mp *y;
+    y = mp_sqr(spare, x);
+    y = mpbarrett_reduce(mb, y, y);
+    spare = x; x = y;
+    sq--;
+  }
+
+  /* --- Done --- */
+
+  for (i = 0; i < TABSZ; i++)
+    mp_drop(tab[i]);
+  xfree(tab);
+exit:
+  mp_drop(d);
+  mp_drop(spare);
+  return (x);
+}
+
 /*----- Test rig ----------------------------------------------------------*/
 
 #ifdef TEST_RIG