X-Git-Url: https://git.distorted.org.uk/u/mdw/catacomb/blobdiff_plain/f46efa79cd2bb9adc81541f1218965f85a6b2eef..80be023065ced106a4078a36371c135a60d2bd6c:/mpreduce.c

diff --git a/mpreduce.c b/mpreduce.c
index 857549a..0f20ed8 100644
--- a/mpreduce.c
+++ b/mpreduce.c
@@ -1,6 +1,6 @@
 /* -*-c-*-
  *
- * $Id: mpreduce.c,v 1.1 2004/03/27 00:04:46 mdw Exp $
+ * $Id$
  *
  * Efficient reduction modulo nice primes
  *
@@ -27,14 +27,6 @@
  * MA 02111-1307, USA.
  */
 
-/*----- Revision history --------------------------------------------------* 
- *
- * $Log: mpreduce.c,v $
- * Revision 1.1  2004/03/27 00:04:46  mdw
- * Implement efficient reduction for pleasant-looking primes.
- *
- */
-
 /*----- Header files ------------------------------------------------------*/
 
 #include <mLib/darray.h>
@@ -55,12 +47,12 @@ DA_DECL(instr_v, mpreduce_instr);
  * Arguments:	@gfreduce *r@ = structure to fill in
  *		@mp *x@ = an integer
  *
- * Returns:	---
+ * Returns:	Zero if successful; nonzero on failure.
  *
  * Use:		Initializes a context structure for reduction.
  */
 
-void mpreduce_create(mpreduce *r, mp *p)
+int mpreduce_create(mpreduce *r, mp *p)
 {
   mpscan sc;
   enum { Z = 0, Z1 = 2, X = 4, X0 = 6 };
@@ -72,7 +64,8 @@ void mpreduce_create(mpreduce *r, mp *p)
 
   /* --- Fill in the easy stuff --- */
 
-  assert(MP_ISPOS(p));
+  if (!MP_POSP(p))
+    return (-1);
   d = mp_bits(p);
   r->lim = d/MPW_BITS;
   r->s = d%MPW_BITS;
@@ -126,9 +119,9 @@ void mpreduce_create(mpreduce *r, mp *p)
     }
   }
   if ((DA(&iv)[DA_LEN(&iv) - 1].op & ~1u) == MPRI_SUB) {
-    fprintf(stderr,
-	    "mpreduce can't handle moduli of the form 2^m + x\n");
-    abort();
+    mp_drop(r->p);
+    DA_DESTROY(&iv);
+    return (-1);
   }
 
 #undef INSTR
@@ -157,7 +150,8 @@ void mpreduce_create(mpreduce *r, mp *p)
       r->iv[i + r->in].argy = b;
     }
   }
-  DA_DESTROY(&iv);  
+  DA_DESTROY(&iv);
+  return (0);
 }
 
 /* --- @mpreduce_destroy@ --- *
@@ -261,7 +255,7 @@ mp *mpreduce_do(mpreduce *r, mp *d, mp *x)
 
   /* --- If source is negative, divide --- */
 
-  if (MP_ISNEG(x)) {
+  if (MP_NEGP(x)) {
     mp_div(0, &d, x, r->p);
     return (d);
   }
@@ -291,9 +285,9 @@ mp *mpreduce_do(mpreduce *r, mp *d, mp *x)
 	*vl = 0;
 	run(r->iv, il, vl, z);
 #ifdef DEBUG
-  MP_PRINTX("x", x);
-  mp_div(0, &_rr, x, r->p);
-  assert(MP_EQ(_r, _rr));
+	MP_PRINTX("x", x);
+	mp_div(0, &_rr, x, r->p);
+	assert(MP_EQ(_r, _rr));
 #endif
       }
     }
@@ -303,9 +297,9 @@ mp *mpreduce_do(mpreduce *r, mp *d, mp *x)
 	*vl &= ((1 << r->s) - 1);
 	run(r->iv + r->in, il + r->in, vl, z);
 #ifdef DEBUG
-  MP_PRINTX("x", x);
-  mp_div(0, &_rr, x, r->p);
-  assert(MP_EQ(_r, _rr));
+	MP_PRINTX("x", x);
+	mp_div(0, &_rr, x, r->p);
+	assert(MP_EQ(_r, _rr));
 #endif
       }
     }
@@ -343,12 +337,18 @@ mp *mpreduce_exp(mpreduce *mr, mp *d, mp *a, mp *e)
   mp *spare = (e->f & MP_BURN) ? MP_NEWSEC : MP_NEW;
 
   MP_SHRINK(e);
-  if (!MP_LEN(e))
+  MP_COPY(a);
+  if (MP_ZEROP(e))
     ;
-  else if (MP_LEN(e) < EXP_THRESH)
-    EXP_SIMPLE(x, a, e);
-  else
-    EXP_WINDOW(x, a, e);
+  else {
+    if (MP_NEGP(e))
+      a = mp_modinv(a, a, mr->p);
+    if (MP_LEN(e) < EXP_THRESH)
+      EXP_SIMPLE(x, a, e);
+    else
+      EXP_WINDOW(x, a, e);
+  }
+  mp_drop(a);
   mp_drop(d);
   mp_drop(spare);
   return (x);