author Mark Wooding Fri, 9 Aug 2013 09:30:48 +0000 (10:30 +0100) committer Mark Wooding Fri, 9 Aug 2013 12:51:04 +0000 (13:51 +0100)
There was an unfortunate amount of code duplication in the old code,
particularly around the issuing of LSR instruction sequences.  Gather
the relevant state into a structure, and split the repeated code into
separate functions.

Also add a bunch of commentary explaining the underlying theory.

 math/gfreduce.c patch | blob | blame | history math/t/gfreduce patch | blob | blame | history

index aed3f22..96e0ced 100644 (file)
@@ -47,26 +47,22 @@ DA_DECL(instr_v, gfreduce_instr);
*
* Let's face it, @gfx_div@ sucks.  It works (I hope), but it's not in any
* sense fast.  Here, we do efficient reduction modulo sparse polynomials.
- *
- * Suppose we have a polynomial @X@ we're trying to reduce mod @P@.  If we
- * take the topmost nonzero word of @X@, call it @w@, then we can eliminate
- * it by subtracting off %$w P x^{k}$% for an appropriate value of @k@.  The
- * trick is in observing that if @P@ is sparse we can do this multiplication
- * and subtraction efficiently, just by XORing appropriate shifts of @w@ into
- * @X@.
- *
- * The first tricky bit is in working out when to stop.  I'll use eight-bit
- * words to demonstrate what I'm talking about.
- *
- *  xxxxxx xxxxxxxx xxxxxxxx xxxxxxxx xxxxxxxx xxxxxxxx
- *                 001ppppp pppppppp pppppppp pppppppp
- *                   |<rp>|
- *                   |<------------ bp ------------->|
- *                 |<------------ nw --------------->|
- *
- * The trick of taking whole words off of @X@ stops working when there are
- * only @nw@ words left.  Then we have to mask off the bottom bits of @w@
- * before continuing.
+ * (It works for arbitrary polynomials, but isn't efficient for dense ones.)
+ *
+ * Suppose that %$p(x) = x^n + p'(x) = \sum_{0\le i<n} p_i x^i$%, hopefully
+ * with only a few other %$p_i \ne 0$%.  We're going to compile %$p$% into a
+ * sequence of instructions which can be used to perform reduction modulo
+ * %$p$%.  The important observation is that %$x^n \equiv p' \pmod p$%.
+ *
+ * Suppose we're working with %$w$%-bit words; let %$n = N w + n'$% with
+ * %$0 \le n' < w$%.  Let %$u(x)$% be some arbitrary polynomial.  Write
+ * %$u = z x^k + u'$% with %$\deg u' < k \ge n$%; then a reduction step uses
+ * that %$u \equiv u' + z p' x^{k-n} \pmod p$%: the right hand side has
+ * degree %$\max \{ \deg u', k + \deg p' - n + \deg z \} < \deg u$%, so this
+ * makes progress towards a complete reduction.
+ *
+ * The compiled instruction sequence computes
+ * %$u' + z p' x^{k-n} = u' + \sum_{0\le i<n} z x^{k-n+i}$%.
*/

/* --- @gfreduce_create@ --- *
@@ -79,16 +75,88 @@ DA_DECL(instr_v, gfreduce_instr);
* Use:                Initializes a context structure for reduction.
*/

+struct gen {
+  unsigned f;                          /* Flags */
+#define f_lsr 1u                       /*   Overflow from previous word */
+#define f_load 2u                      /*   Outstanding @LOAD@ */
+  instr_v iv;                          /* Instruction vector */
+  size_t w;                            /* Currently loaded target word */
+  size_t wi;                           /* Left-shifts for current word */
+};
+
+#define INSTR(g_, op_, arg_) do {                                      \
+  struct gen *_g = (g_);                                               \
+  instr_v *_iv = &_g->iv;                                              \
+  size_t _i = DA_LEN(_iv);                                             \
+                                                                       \
+  DA_ENSURE(_iv, 1);                                                   \
+  DA(_iv)[_i].op = (op_);                                              \
+  DA(_iv)[_i].arg = (arg_);                                            \
+  DA_EXTEND(_iv, 1);                                                   \
+} while (0)
+
+static void emit_load(struct gen *g, size_t w)
+{
+  INSTR(g, GFRI_LOAD, w);
+  g->f |= f_load;
+  g->w = w;
+}
+
+static void emit_right_shifts(struct gen *g)
+{
+  gfreduce_instr *ip;
+  size_t i, wl;
+
+  /* --- Close off the current word --- *
+   *
+   * If we shifted into this current word with a nonzero bit offset, then
+   * we'll also need to arrange to perform a sequence of right shifts into
+   * the following word, which we might as well do by scanning the
+   * instruction sequence (which starts at @wi@).
+   *
+   * Either way, we leave a @LOAD@ unmatched if there was one before, in the
+   * hope that callers have an easier time; @g->w@ is updated to reflect the
+   * currently open word.
+   */
+
+  if (!(g->f & f_lsr))
+    return;
+
+  wl = DA_LEN(&g->iv);
+  INSTR(g, GFRI_STORE, g->w);
+  emit_load(g, g->w - 1);
+  for (i = g->wi; i < wl; i++) {
+    ip = &DA(&g->iv)[i];
+    assert(ip->op == GFRI_LSL);
+    if (ip->arg)
+      INSTR(g, GFRI_LSR, MPW_BITS - ip->arg);
+  }
+  g->f &= ~f_lsr;
+}
+
+static void ensure_loaded(struct gen *g, size_t w)
+{
+  if (!(g->f & f_load)) {
+    emit_load(g, w);
+    g->wi = DA_LEN(&g->iv);
+  } else if (w != g->w) {
+    emit_right_shifts(g);
+    if (w != g->w) {
+      INSTR(g, GFRI_STORE, g->w);
+      emit_load(g, w);
+    }
+    g->wi = DA_LEN(&g->iv);
+  }
+}
+
void gfreduce_create(gfreduce *r, mp *p)
{
-  instr_v iv = DA_INIT;
+  struct gen g = { 0, DA_INIT };
unsigned long d;
unsigned dw;
mpscan sc;
unsigned long i;
-  gfreduce_instr *ip;
-  unsigned f = 0;
-  size_t w, ww, wi, wl, ll, bb;
+  size_t w, bb;

/* --- Sort out the easy stuff --- */

@@ -103,79 +171,80 @@ void gfreduce_create(gfreduce *r, mp *p)
}
r->p = mp_copy(p);

-  /* --- Stash a new instruction --- */
-
-#define INSTR(op_, arg_) do {                                          \
-  DA_ENSURE(&iv, 1);                                                   \
-  DA(&iv)[DA_LEN(&iv)].op = (op_);                                     \
-  DA(&iv)[DA_LEN(&iv)].arg = (arg_);                                   \
-  DA_EXTEND(&iv, 1);                                                   \
-} while (0)
-
-#define f_lsr 1u
+  /* --- How this works --- *
+   *
+   * The instruction sequence is run with two ambient parameters: a pointer
+   * (usually) just past the most significant word of the polynomial to be
+   * reduced; and a word %$z$% which is the multiple of %$p'$% we are meant
+   * to add.
+   *
+   * The sequence visits each word of the polynomial at most once.  Suppose
+   * %$u = z x^{w N} + u'$%; our pointer points just past the end of %$u'$%.
+   * Word %$I$% of %$u'$% will be affected by modulus bits %$p_i$% where
+   * %$(N - I - 1) w + 1 \le i \le (N - I + 1) w - 1$%, so %$p_i$% affects
+   * word %$I = \lceil (n - i + 1)/w \rceil$% and (if %$i$% is not a multiple
+   * of %$w$%) also word %$I - 1$%.
+   *
+   * We have four instructions: @LOAD@ reads a specified word of %$u$% into an
+   * accumulator, and @STORE@ stores it back (we'll always store back to the
+   * same word we most recently read, but this isn't a requirement); and
+   * @LSL@ and @LSR@, which XOR in appropriately shifted copies of %$z$% into
+   * the accumulator.  So a typical program will contain sequences of @LSR@
+   * and @LSL@ instructions sandwiched between @LOAD@/@STORE@ pairs.
+   *
+   * We do a single right-to-left pass across %$p$%.
+   */

-  w = (d + MPW_BITS - 1)/MPW_BITS;
-  INSTR(GFRI_LOAD, w);
-  wi = DA_LEN(&iv);
-  f = 0;
-  ll = 0;
bb = MPW_BITS - dw;
+
for (i = 0, mp_scan(&sc, p); mp_step(&sc) && i < d; i++) {
if (!mp_bit(&sc))
continue;
-    ww = (d - i + MPW_BITS - 1)/MPW_BITS;
-    if (ww != w) {
-      wl = DA_LEN(&iv);
-      INSTR(GFRI_STORE, w);
-      if (!ll)
-       ll = DA_LEN(&iv);
-      if (!(f & f_lsr))
-       INSTR(GFRI_LOAD, ww);
-      else {
-       INSTR(GFRI_LOAD, w - 1);
-       for (; wi < wl; wi++) {
-         ip = &DA(&iv)[wi];
-         assert(ip->op == GFRI_LSL);
-         if (ip->arg)
-           INSTR(GFRI_LSR, MPW_BITS - ip->arg);
-       }
-       if (w - 1 != ww) {
-         INSTR(GFRI_STORE, w - 1);
-         INSTR(GFRI_LOAD, ww);
-       }
-       f &= ~f_lsr;
-      }
-      w = ww;
-      wi = DA_LEN(&iv);
-    }
-    INSTR(GFRI_LSL, (bb + i)%MPW_BITS);
+
+    /* --- We've found a set bit, so work out which word it affects --- *
+     *
+     * In general, a bit affects two words: it needs to be shifted left into
+     * one, and shifted right into the next.  We find the former here.
+     */
+
+    w = (d - i + MPW_BITS - 1)/MPW_BITS;
+
+    /* --- Concentrate on the appropriate word --- */
+
+    ensure_loaded(&g, w);
+
+    /* --- Accumulate a new @LSL@ instruction --- *
+     *
+     * If this was a nonzero shift, then we'll need to arrange to do right
+     * shifts into the following word.
+     */
+
+    INSTR(&g, GFRI_LSL, (bb + i)%MPW_BITS);
if ((bb + i)%MPW_BITS)
-      f |= f_lsr;
-  }
-  wl = DA_LEN(&iv);
-  INSTR(GFRI_STORE, w);
-  if (!ll)
-    ll = DA_LEN(&iv);
-  if (f & f_lsr) {
-    INSTR(GFRI_LOAD, w - 1);
-    for (; wi < wl; wi++) {
-      ip = &DA(&iv)[wi];
-      assert(ip->op == GFRI_LSL);
-      if (ip->arg)
-       INSTR(GFRI_LSR, MPW_BITS - ip->arg);
-    }
-    INSTR(GFRI_STORE, w - 1);
+      g.f |= f_lsr;
}

-#undef INSTR
+  /* --- Wrapping up --- *
+   *
+   * We probably need a final @STORE@, and maybe a sequence of right shifts.
+   */

-  r->in = DA_LEN(&iv);
+  if (g.f & f_load) {
+    emit_right_shifts(&g);
+    INSTR(&g, GFRI_STORE, g.w);
+  }
+
+  r->in = DA_LEN(&g.iv);
r->iv = xmalloc(r->in * sizeof(gfreduce_instr));
-  r->liv = r->iv + ll;
-  memcpy(r->iv, DA(&iv), r->in * sizeof(gfreduce_instr));
-  DA_DESTROY(&iv);
+  memcpy(r->iv, DA(&g.iv), r->in * sizeof(gfreduce_instr));
+  DA_DESTROY(&g.iv);
}

+#undef INSTR
+
+#undef f_lsr
+#undef f_load
+
/* --- @gfreduce_destroy@ --- *
*
* Arguments:  @gfreduce *r@ = structure to free
index f712316..ecece92 100644 (file)
@@ -13,6 +13,12 @@ reduce {
0xb2ca471b0867d5fae2e4f27a2d2706da
0xf254423fef93d5d7a76ecf22c656c1352c53257875945d33
0x582f783fc210f72814780e69b0bd29ff;
+
+  # --- Bugs discovered ---
+
+  0x20000000000000000000000000000000000000004000000000000000001
+    0x110414154054140445011511541540514401111414505115044104145451505001151441450000541004550554154040500000411050400055041
+    0x1bf878a39fbee9cf20a2f6f41eadda756518f1669b2c7d8f9234965b6b3;
}

modexp {