X-Git-Url: https://git.distorted.org.uk/~mdw/catacomb/blobdiff_plain/e85f65dc3662443fa2172bdd45cd97ecb1e1eea6..HEAD:/math/gfreduce.c

diff --git a/math/gfreduce.c b/math/gfreduce.c
index aed3f224..f29c1048 100644
--- a/math/gfreduce.c
+++ b/math/gfreduce.c
@@ -47,26 +47,23 @@ 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^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@ --- *
@@ -79,19 +76,113 @@ 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@ */
+#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;
@@ -103,79 +194,92 @@ 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;
+      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
@@ -193,7 +297,7 @@ void gfreduce_destroy(gfreduce *r)
 
 /* --- @gfreduce_dump@ --- *
  *
- * Arguments:	@gfreduce *r@ = structure to dump
+ * Arguments:	@const gfreduce *r@ = structure to dump
  *		@FILE *fp@ = file to dump on
  *
  * Returns:	---
@@ -201,7 +305,7 @@ void gfreduce_destroy(gfreduce *r)
  * Use:		Dumps a reduction context.
  */
 
-void gfreduce_dump(gfreduce *r, FILE *fp)
+void gfreduce_dump(const gfreduce *r, FILE *fp)
 {
   size_t i;
 
@@ -210,16 +314,20 @@ 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@ --- *
  *
- * Arguments:	@gfreduce *r@ = reduction context
+ * Arguments:	@const gfreduce *r@ = reduction context
  *		@mp *d@ = destination
  *		@mp *x@ = source
  *
@@ -242,7 +350,7 @@ static void run(const gfreduce_instr *i, const gfreduce_instr *il,
   }
 }
 
-mp *gfreduce_do(gfreduce *r, mp *d, mp *x)
+mp *gfreduce_do(const gfreduce *r, mp *d, mp *x)
 {
   mpw *v, *vl;
   const gfreduce_instr *il;
@@ -271,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);
       }
     }
   }
@@ -284,20 +392,26 @@ mp *gfreduce_do(gfreduce *r, mp *d, mp *x)
 
 /* --- @gfreduce_sqrt@ --- *
  *
- * Arguments:	@gfreduce *r@ = pointer to reduction context
+ * Arguments:	@const gfreduce *r@ = pointer to reduction context
  *		@mp *d@ = destination
  *		@mp *x@ = some polynomial
  *
  * Returns:	The square root of @x@ modulo @r->p@, or null.
  */
 
-mp *gfreduce_sqrt(gfreduce *r, mp *d, mp *x)
+mp *gfreduce_sqrt(const gfreduce *r, mp *d, mp *x)
 {
   mp *y = MP_COPY(x);
   mp *z, *spare = MP_NEW;
   unsigned long m = mp_bits(r->p) - 1;
   unsigned long i;
 
+  /* --- This is pretty easy --- *
+   *
+   * Note that %$x = x^{2^m}$%; therefore %$(x^{2^{m-1}})^2 = x^{2^m} = x$%,
+   * so %$x^{2^{m-1}}$% is the square root we seek.
+   */
+
   for (i = 0; i < m - 1; i++) {
     mp *t = gf_sqr(spare, y);
     spare = y;
@@ -316,14 +430,17 @@ mp *gfreduce_sqrt(gfreduce *r, mp *d, mp *x)
 
 /* --- @gfreduce_trace@ --- *
  *
- * Arguments:	@gfreduce *r@ = pointer to reduction context
+ * Arguments:	@const gfreduce *r@ = pointer to reduction context
  *		@mp *x@ = some polynomial
  *
  * Returns:	The trace of @x@. (%$\Tr(x)=x + x^2 + \cdots + x^{2^{m-1}}$%
- *		if %$x \in \gf{2^m}$%).
+ *		if %$x \in \gf{2^m}$%).  Since the trace is invariant under
+ *		the Frobenius automorphism (i.e., %$\Tr(x)^2 = \Tr(x)$%), it
+ *		must be an element of the base field, i.e., %$\gf{2}$%, and
+ *		we only need a single bit to represent it.
  */
 
-int gfreduce_trace(gfreduce *r, mp *x)
+int gfreduce_trace(const gfreduce *r, mp *x)
 {
   mp *y = MP_COPY(x);
   mp *spare = MP_NEW;
@@ -345,7 +462,7 @@ int gfreduce_trace(gfreduce *r, mp *x)
 
 /* --- @gfreduce_halftrace@ --- *
  *
- * Arguments:	@gfreduce *r@ = pointer to reduction context
+ * Arguments:	@const gfreduce *r@ = pointer to reduction context
  *		@mp *d@ = destination
  *		@mp *x@ = some polynomial
  *
@@ -354,7 +471,7 @@ int gfreduce_trace(gfreduce *r, mp *x)
  *		if %$x \in \gf{2^m}$% with %$m$% odd).
  */
 
-mp *gfreduce_halftrace(gfreduce *r, mp *d, mp *x)
+mp *gfreduce_halftrace(const gfreduce *r, mp *d, mp *x)
 {
   mp *y = MP_COPY(x);
   mp *spare = MP_NEW;
@@ -377,27 +494,82 @@ mp *gfreduce_halftrace(gfreduce *r, mp *d, mp *x)
 
 /* --- @gfreduce_quadsolve@ --- *
  *
- * Arguments:	@gfreduce *r@ = pointer to reduction context
+ * Arguments:	@const gfreduce *r@ = pointer to reduction context
  *		@mp *d@ = destination
  *		@mp *x@ = some polynomial
  *
- * Returns:	A polynomial @y@ such that %$y^2 + y = x$%, or null.
+ * Returns:	A polynomial @z@ such that %$z^2 + z = x$%, or null.
+ *
+ * Use:		Solves quadratic equations in a field with characteristic 2.
+ *		Suppose we have an equation %$y^2 + A y + B = 0$% where
+ *		%$A \ne 0$%.  (If %$A = 0$% then %$y = \sqrt{B}$% and you
+ *		want @gfreduce_sqrt@ instead.)  Use this function to solve
+ *		%$z^2 + z = B/A^2$%; then set %$y = A z$%, since
+ *		%$y^2 + y = A^2 z^2 + A^2 z = A^2 (z^2 + z) = B$% as
+ *		required.
+ *
+ *		The two roots are %$z$% and %$z + 1$%; this function always
+ *		returns the one with zero scalar coefficient.
  */
 
-mp *gfreduce_quadsolve(gfreduce *r, mp *d, mp *x)
+mp *gfreduce_quadsolve(const gfreduce *r, mp *d, mp *x)
 {
   unsigned long m = mp_bits(r->p) - 1;
   mp *t;
 
+  /* --- About the solutions --- *
+   *
+   * Factor %$z^2 + z = z (z + 1)$%.  Therefore, if %$z^2 + z = x$% and
+   * %$z' = z + 1$% then %$z'^2 + z' = z^2 + 1 + z + 1 = z^2 + z$%, so
+   * %$z + 1$% is the other solution.
+   *
+   * A solution exists if and only if %$\Tr(x) = 0$%.  To see the `only if'
+   * implication, recall that the trace function is linear, and hence
+   * $%\Tr(z^2 + z) = \Tr(z)^2 + \Tr(z) = \Tr(z) + \Tr(z) = 0$%.  The `if'
+   * direction will be proven using explicit constructions captured in the
+   * code below.
+   */
+
   MP_COPY(x);
-  if (m & 1)
+  if (m & 1) {
+
+    /* --- A short-cut for fields with odd degree ---
+     *
+     * The method below works in all binary fields, but there's a quicker way
+     * which works whenever the degree is odd.  The half-trace is
+     * %$z = \sum_{0\le i\le (m-1)/2} x^{2^{2i}}$%.  Then %$z^2 + z = {}$%
+     * %$\sum_{0\le i\le (m-1)/2} (x^{2^{2i}} + x^{2^{2i+1}}) = {}$%
+     * %$\Tr(x) + x^{2^m} = \Tr(x) + x$%.  This therefore gives us the
+     * solution we want whenever %$\Tr(x) = 0$%.
+     */
+
     d = gfreduce_halftrace(r, d, x);
-  else {
+  } else {
     mp *z, *w, *rho = MP_NEW;
     mp *spare = MP_NEW;
     grand *fr = fibrand_create(0);
     unsigned long i;
 
+    /* --- Unpicking the magic --- *
+     *
+     * Choose %$\rho \inr \gf{2^m}$% with %$\Tr(\rho) = 1$%.  Let
+     * %$z = \sum_{0\le i<m} \rho^{2^i} \sum_{0\le j<i} x^{2^j} = {}$%
+     * %$\sum_{1\le i<m} \rho^{2^i} (x + \sum_{1\le j<i} x^{2^j} = {}$%
+     * %$\rho^2 x + \rho^4 (x + x^2) + \rho^8 (x + x^2 + x^4) + \cdots + {}$%
+     * %$\rho^{2^{m-1}} (x + x^2 + x^{2^{m-2}})$%.  Then %$z^2 = {}$%
+     * %$\sum_{0\le i<m} \rho^{2^{i+1}} \sum_{0\le j<i} x^{2^{j+1}} = {}$%
+     * %$\sum_{1\le i\le m} \rho^{2^i} \sum_{1\le j<i} x^{2^j} = {}$%
+     * %$\sum_{1\le i<m} \rho^{2^i} \sum_{1\le j<i} x^{2^j} + {}$%
+     * %$\rho^{2^m} \sum_{1\le j<m} x^{2^j}$%; and, somewhat miraculously,
+     * %$z^2 + z = \sum_{1\le i<m} \rho^{2^i} x + {}$%
+     * %$\rho \sum_{1\le i<m} x^{2^i} = x (\Tr(\rho) + \rho) + {}$%
+     * %$\rho (\Tr(x) + x) = x \Tr(\rho) + \rho \Tr(x)$%.  Again,
+     * this gives us the root we want whenever %$\Tr(x) = 0$%.
+     *
+     * The loop below calculates %$w = \Tr(\rho)$% and %$z$% simultaneously,
+     * since the same powers of %$\rho$% are wanted in both calculations.
+     */
+
     for (;;) {
       rho = mprand(rho, m, fr, 0);
       z = MP_ZERO;
@@ -422,6 +594,12 @@ mp *gfreduce_quadsolve(gfreduce *r, mp *d, mp *x)
     d = z;
   }
 
+  /* --- Check that we calculated the right answer --- *
+   *
+   * It should be correct; if it's not then maybe the ring we're working in
+   * isn't really a field.
+   */
+
   t = gf_sqr(MP_NEW, d); t = gfreduce_do(r, t, t); t = gf_add(t, t, d);
   if (!MP_EQ(t, x)) {
     MP_DROP(d);
@@ -429,13 +607,20 @@ mp *gfreduce_quadsolve(gfreduce *r, mp *d, mp *x)
   }
   MP_DROP(t);
   MP_DROP(x);
+
+  /* --- Pick a canonical root --- *
+   *
+   * The two roots are %$z$% and %$z + 1$%; pick the one with a zero
+   * scalar coefficient just for consistency's sake.
+   */
+
   if (d) d->v[0] &= ~(mpw)1;
   return (d);
 }
 
 /* --- @gfreduce_exp@ --- *
  *
- * Arguments:	@gfreduce *gr@ = pointer to reduction context
+ * Arguments:	@const gfreduce *gr@ = pointer to reduction context
  *		@mp *d@ = fake destination
  *		@mp *a@ = base
  *		@mp *e@ = exponent
@@ -443,7 +628,7 @@ mp *gfreduce_quadsolve(gfreduce *r, mp *d, mp *x)
  * Returns:	Result, %$a^e \bmod m$%.
  */
 
-mp *gfreduce_exp(gfreduce *gr, mp *d, mp *a, mp *e)
+mp *gfreduce_exp(const gfreduce *gr, mp *d, mp *a, mp *e)
 {
   mp *x = MP_ONE;
   mp *spare = (e->f & MP_BURN) ? MP_NEWSEC : MP_NEW;
@@ -470,8 +655,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;