*
* 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^n + p'$% where %$p' = \sum_{0\le i<n} p_i x^i$%,
+ * hopefully with only a few %$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@ --- *
* 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@ */
+#define f_fip 4u /* Final-pass offset is set */
+ instr_v iv; /* Instruction vector */
+ size_t fip; /* Offset for final-pass reduction */
+ size_t w; /* Currently loaded target word */
+ size_t wi; /* Left-shifts for current word */
+ gfreduce *r; /* Reduction context pointer */
+};
+
+#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)
+{
+ /* --- If this is not the low-order word then note final-pass start --- *
+ *
+ * Once we've eliminated the whole high-degree words, there will possibly
+ * remain a few high-degree bits. We can further reduce the subject
+ * polynomial by subtracting an appropriate multiple of %$p'$%, but if we
+ * do this naively we'll end up addressing `low-order' words beyond the
+ * bottom of our input. We solve this problem by storing an alternative
+ * start position for this final pass (which works because we scan bits
+ * right-to-left).
+ */
+
+ if (!(g->f & f_fip) && w < g->r->lim) {
+ g->fip = DA_LEN(&g->iv);
+ g->f |= f_fip;
+ }
+
+ /* --- Actually emit the instruction --- */
+
+ 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 --- */
+ g.r = r;
d = mp_bits(p); assert(d); d--;
r->lim = d/MPW_BITS;
dw = d%MPW_BITS;
}
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;
+ g.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);
+
+ /* --- Wrapping up --- *
+ *
+ * We probably need a final @STORE@, and maybe a sequence of right shifts.
+ */
+
+ if (g.f & f_load) {
+ emit_right_shifts(&g);
+ INSTR(&g, GFRI_STORE, g.w);
}
-#undef INSTR
+ /* --- Copy the instruction vector.
+ *
+ * If we've not set a final-pass offset yet then now would be an excellent
+ * time. Obviously it should be right at the end, because there's nothing
+ * for a final pass to do.
+ */
- r->in = DA_LEN(&iv);
+ 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));
+
+ if (!(g.f & f_fip)) g.fip = DA_LEN(&g.iv);
+ r->fiv = r->iv + g.fip;
+
+ DA_DESTROY(&g.iv);
}
+#undef INSTR
+
+#undef f_lsr
+#undef f_load
+#undef f_fip
+
/* --- @gfreduce_destroy@ --- *
*
* Arguments: @gfreduce *r@ = structure to free
(unsigned long)r->lim, (unsigned long)r->mask);
for (i = 0; i < r->in; i++) {
static const char *opname[] = { "load", "lsl", "lsr", "store" };
+ if (&r->iv[i] == r->fiv)
+ fputs("final:\n", fp);
assert(r->iv[i].op < N(opname));
fprintf(fp, " %s %lu\n",
opname[r->iv[i].op],
(unsigned long)r->iv[i].arg);
}
+ if (&r->iv[i] == r->fiv)
+ fputs("final:\n", fp);
}
/* --- @gfreduce_do@ --- *
while (*vl & r->mask) {
z = *vl & r->mask;
*vl &= ~r->mask;
- run(r->iv, il, vl, z);
+ run(r->fiv, il, vl, z);
}
}
}
#ifdef TEST_RIG
-#define MP(x) mp_readstring(MP_NEW, #x, 0, 0)
-
static int vreduce(dstr *v)
{
mp *d = *(mp **)v[0].buf;