X-Git-Url: https://git.distorted.org.uk/~mdw/catacomb/blobdiff_plain/192053f16dc688ff9d73d42b52bf35f39948c0e9..HEAD:/math/mpreduce.c

diff --git a/math/mpreduce.c b/math/mpreduce.c
index cc74bc2b..0fd4b458 100644
--- a/math/mpreduce.c
+++ b/math/mpreduce.c
@@ -47,7 +47,7 @@ DA_DECL(instr_v, mpreduce_instr);
  *
  * Suppose %$x = x' + z 2^k$%, where %$k \ge n$%; then
  * %$x \equiv x' + d z 2^{k-n} \pmod p$%.  We can use this to trim the
- * representation of %$x$%; each time, we reduce %$x$% by a mutliple of
+ * representation of %$x$%; each time, we reduce %$x$% by a multiple of
  * %$2^{k-n} p$%.  We can do this in two passes: firstly by taking whole
  * words off the top, and then (if necessary) by trimming the top word.
  * Finally, if %$p \le x < 2^n$% then %$0 \le x - p < p$% and we're done.
@@ -81,7 +81,10 @@ DA_DECL(instr_v, mpreduce_instr);
  * Arguments:	@gfreduce *r@ = structure to fill in
  *		@mp *x@ = an integer
  *
- * Returns:	Zero if successful; nonzero on failure.
+ * Returns:	Zero if successful; nonzero on failure.  The current
+ *		algorithm always succeeds when given positive @x@.  Earlier
+ *		versions used to fail on particular kinds of integers, but
+ *		this is guaranteed not to happen any more.
  *
  * Use:		Initializes a context structure for reduction.
  */
@@ -196,18 +199,27 @@ int mpreduce_create(mpreduce *r, mp *p)
     }
   }
 
-  /* --- This doesn't always work --- *
+  /* --- Fix up wrong-sided decompositions --- *
    *
-   * If %$d \ge 2^{n-1}$% then the above recurrence will output a subtraction
-   * as the final instruction, which may sometimes underflow.  (It interprets
-   * such numbers as being in the form %$2^{n-1} + d$%.)  This is clearly
-   * bad, so detect the situation and fail gracefully.
+   * At this point, we haven't actually finished up the state machine
+   * properly.  We stopped scanning just after bit %$n - 1$% -- the most
+   * significant one, which we know in advance must be set (since @x@ is
+   * strictly positive).  Therefore we are either in state @X@ or @Z1@.  In
+   * the former case, we have nothing to do.  In the latter, there are two
+   * subcases to deal with.  If there are no other instructions, then @x@ is
+   * a perfect power of two, and %$d = 0$%, so again there is nothing to do.
+   *
+   * In the remaining case, we have decomposed @x@ as %$2^{n-1} + d$%, for
+   * some positive %$d%, which is unfortunate: if we're asked to reduce
+   * %$2^n$%, say, we'll end up with %$-d$% (or would do, if we weren't
+   * sticking to unsigned arithmetic for good performance).  So instead, we
+   * rewrite this as %$2^n - 2^{n-1} + d$% and everything will be good.
    */
 
-  if (DA_LEN(&iv) && (DA(&iv)[DA_LEN(&iv) - 1].op & ~1u) == MPRI_SUB) {
-    mp_drop(r->p);
-    DA_DESTROY(&iv);
-    return (-1);
+  if (st == Z1 && DA_LEN(&iv)) {
+    w = 1;
+    b = (bb + d)%MPW_BITS;
+    INSTR(MPRI_ADD | !!b, w, b);
   }
 
 #undef INSTR
@@ -283,7 +295,7 @@ void mpreduce_destroy(mpreduce *r)
 
 /* --- @mpreduce_dump@ --- *
  *
- * Arguments:	@mpreduce *r@ = structure to dump
+ * Arguments:	@const mpreduce *r@ = structure to dump
  *		@FILE *fp@ = file to dump on
  *
  * Returns:	---
@@ -291,7 +303,7 @@ void mpreduce_destroy(mpreduce *r)
  * Use:		Dumps a reduction context.
  */
 
-void mpreduce_dump(mpreduce *r, FILE *fp)
+void mpreduce_dump(const mpreduce *r, FILE *fp)
 {
   size_t i;
   static const char *opname[] = { "add", "addshift", "sub", "subshift" };
@@ -319,7 +331,7 @@ void mpreduce_dump(mpreduce *r, FILE *fp)
 
 /* --- @mpreduce_do@ --- *
  *
- * Arguments:	@mpreduce *r@ = reduction context
+ * Arguments:	@const mpreduce *r@ = reduction context
  *		@mp *d@ = destination
  *		@mp *x@ = source
  *
@@ -341,7 +353,7 @@ static void run(const mpreduce_instr *i, const mpreduce_instr *il,
   }
 }
 
-mp *mpreduce_do(mpreduce *r, mp *d, mp *x)
+mp *mpreduce_do(const mpreduce *r, mp *d, mp *x)
 {
   mpw *v, *vl;
   const mpreduce_instr *il;
@@ -398,7 +410,7 @@ mp *mpreduce_do(mpreduce *r, mp *d, mp *x)
 
 /* --- @mpreduce_exp@ --- *
  *
- * Arguments:	@mpreduce *mr@ = pointer to reduction context
+ * Arguments:	@const mpreduce *mr@ = pointer to reduction context
  *		@mp *d@ = fake destination
  *		@mp *a@ = base
  *		@mp *e@ = exponent
@@ -406,7 +418,7 @@ mp *mpreduce_do(mpreduce *r, mp *d, mp *x)
  * Returns:	Result, %$a^e \bmod m$%.
  */
 
-mp *mpreduce_exp(mpreduce *mr, mp *d, mp *a, mp *e)
+mp *mpreduce_exp(const mpreduce *mr, mp *d, mp *a, mp *e)
 {
   mp *x = MP_ONE;
   mp *spare = (e->f & MP_BURN) ? MP_NEWSEC : MP_NEW;
@@ -431,11 +443,8 @@ mp *mpreduce_exp(mpreduce *mr, mp *d, mp *a, mp *e)
 
 /*----- Test rig ----------------------------------------------------------*/
 
-
 #ifdef TEST_RIG
 
-#define MP(x) mp_readstring(MP_NEW, #x, 0, 0)
-
 static int vreduce(dstr *v)
 {
   mp *d = *(mp **)v[0].buf;