X-Git-Url: https://git.distorted.org.uk/~mdw/catacomb/blobdiff_plain/8e663045cd18e14e99affd3f904170dc686fd76c..e99c5bdc78579b1c829a6ee852a68389cba52384:/math/gfreduce.c

diff --git a/math/gfreduce.c b/math/gfreduce.c
index 96e0cedf..582b723f 100644
--- a/math/gfreduce.c
+++ b/math/gfreduce.c
@@ -49,14 +49,15 @@ DA_DECL(instr_v, gfreduce_instr);
  * 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.
@@ -79,9 +80,12 @@ 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 {					\
@@ -97,6 +101,24 @@ struct gen {
 
 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;
@@ -160,6 +182,7 @@ void gfreduce_create(gfreduce *r, mp *p)
 
   /* --- Sort out the easy stuff --- */
 
+  g.r = r;
   d = mp_bits(p); assert(d); d--;
   r->lim = d/MPW_BITS;
   dw = d%MPW_BITS;
@@ -234,9 +257,20 @@ void gfreduce_create(gfreduce *r, mp *p)
     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);
 }
 
@@ -244,6 +278,7 @@ void gfreduce_create(gfreduce *r, mp *p)
 
 #undef f_lsr
 #undef f_load
+#undef f_fip
 
 /* --- @gfreduce_destroy@ --- *
  *
@@ -279,11 +314,15 @@ void gfreduce_dump(gfreduce *r, FILE *fp)
 	  (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@ --- *
@@ -340,7 +379,7 @@ mp *gfreduce_do(gfreduce *r, mp *d, mp *x)
       while (*vl & r->mask) {
 	z = *vl & r->mask;
 	*vl &= ~r->mask;
-	run(r->iv, il, vl, z);
+	run(r->fiv, il, vl, z);
       }
     }
   }
@@ -539,8 +578,6 @@ mp *gfreduce_exp(gfreduce *gr, mp *d, mp *a, mp *e)
 
 #ifdef TEST_RIG
 
-#define MP(x) mp_readstring(MP_NEW, #x, 0, 0)
-
 static int vreduce(dstr *v)
 {
   mp *d = *(mp **)v[0].buf;