* sense fast. Here, we do efficient reduction modulo sparse polynomials.
* (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 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
+ * %$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.
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 { \
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;
/* --- Sort out the easy stuff --- */
+ g.r = r;
d = mp_bits(p); assert(d); d--;
r->lim = d/MPW_BITS;
dw = d%MPW_BITS;
INSTR(&g, GFRI_STORE, g.w);
}
+ /* --- 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(&g.iv);
r->iv = xmalloc(r->in * sizeof(gfreduce_instr));
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 f_lsr
#undef f_load
+#undef f_fip
/* --- @gfreduce_destroy@ --- *
*
(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;