X-Git-Url: https://git.distorted.org.uk/u/mdw/catacomb/blobdiff_plain/51f5bbe0bfe8a344c3cc1f5a9fddf885fa75d057..HEAD:/math/gfreduce.c

diff --git a/math/gfreduce.c b/math/gfreduce.c
index 96e0ced..45591dd 100644
--- a/math/gfreduce.c
+++ b/math/gfreduce.c
@@ -79,9 +79,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 +100,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 +181,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 +256,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 +277,7 @@ void gfreduce_create(gfreduce *r, mp *p)
 
 #undef f_lsr
 #undef f_load
+#undef f_fip
 
 /* --- @gfreduce_destroy@ --- *
  *
@@ -279,11 +313,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 +378,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);
       }
     }
   }