X-Git-Url: https://git.distorted.org.uk/u/mdw/catacomb/blobdiff_plain/81578196d5732e443c75768ba9118c581c407cc7..59919ae4b1721ca271c3d3e5955c09d322573821:/mpx.c

diff --git a/mpx.c b/mpx.c
index f40d7d9..ef93e3e 100644
--- a/mpx.c
+++ b/mpx.c
@@ -1,6 +1,6 @@
 /* -*-c-*-
  *
- * $Id: mpx.c,v 1.16 2003/05/16 09:09:24 mdw Exp $
+ * $Id: mpx.c,v 1.19 2004/04/03 03:29:40 mdw Exp $
  *
  * Low-level multiprecision arithmetic
  *
@@ -30,6 +30,19 @@
 /*----- Revision history --------------------------------------------------* 
  *
  * $Log: mpx.c,v $
+ * Revision 1.19  2004/04/03 03:29:40  mdw
+ * Fix overrun in @mpx_lsr@.
+ *
+ * Revision 1.18  2004/04/01 12:50:09  mdw
+ * Add cyclic group abstraction, with test code.  Separate off exponentation
+ * functions for better static linking.  Fix a buttload of bugs on the way.
+ * Generally ensure that negative exponents do inversion correctly.  Add
+ * table of standard prime-field subgroups.  (Binary field subgroups are
+ * currently unimplemented but easy to add if anyone ever finds a good one.)
+ *
+ * Revision 1.17  2004/03/27 00:04:46  mdw
+ * Implement efficient reduction for pleasant-looking primes.
+ *
  * Revision 1.16  2003/05/16 09:09:24  mdw
  * Fix @mp_lsl2c@.  Turns out to be surprisingly tricky.
  *
@@ -650,7 +663,7 @@ void mpx_lsr(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, size_t n)
   /* --- Single bit shifting --- */
 
   else if (n == 1) {
-    mpw w = *av++ >> 1;
+    mpw w = av < avl ? *av++ >> 1 : 0;
     while (av < avl) {
       mpw t;
       if (dv >= dvl)
@@ -877,6 +890,30 @@ void mpx_uadd(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl,
 
 void mpx_uaddn(mpw *dv, mpw *dvl, mpw n) { MPX_UADDN(dv, dvl, n); }
 
+/* --- @mpx_uaddnlsl@ --- *
+ *
+ * Arguments:	@mpw *dv, *dvl@ = destination and first argument vector
+ *		@mpw a@ = second argument
+ *		@unsigned o@ = offset in bits
+ *
+ * Returns:	---
+ *
+ * Use:		Computes %$d + 2^o a$%.  If the result overflows then
+ *		high-order bits are discarded, as usual.  We must have
+ *		@0 < o < MPW_BITS@.
+ */
+
+void mpx_uaddnlsl(mpw *dv, mpw *dvl, mpw a, unsigned o)
+{
+  mpd x = (mpd)a << o;
+
+  while (x && dv < dvl) {
+    x += *dv;
+    *dv++ = MPW(x);
+    x >>= MPW_BITS;
+  }
+}
+
 /* --- @mpx_usub@ --- *
  *
  * Arguments:	@mpw *dv, *dvl@ = destination vector base and limit
@@ -931,6 +968,33 @@ void mpx_usub(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl,
 
 void mpx_usubn(mpw *dv, mpw *dvl, mpw n) { MPX_USUBN(dv, dvl, n); }
 
+/* --- @mpx_uaddnlsl@ --- *
+ *
+ * Arguments:	@mpw *dv, *dvl@ = destination and first argument vector
+ *		@mpw a@ = second argument
+ *		@unsigned o@ = offset in bits
+ *
+ * Returns:	---
+ *
+ * Use:		Computes %$d + 2^o a$%.  If the result overflows then
+ *		high-order bits are discarded, as usual.  We must have
+ *		@0 < o < MPW_BITS@.
+ */
+
+void mpx_usubnlsl(mpw *dv, mpw *dvl, mpw a, unsigned o)
+{
+  mpw b = a >> (MPW_BITS - o);
+  a <<= o;
+
+  if (dv < dvl) {
+    mpd x = (mpd)*dv - (mpd)a;
+    *dv++ = MPW(x);
+    if (x >> MPW_BITS)
+      b++;
+    MPX_USUBN(dv, dvl, b);
+  }
+}
+
 /* --- @mpx_umul@ --- *
  *
  * Arguments:	@mpw *dv, *dvl@ = destination vector base and limit
@@ -1125,7 +1189,7 @@ void mpx_udiv(mpw *qv, mpw *qvl, mpw *rv, mpw *rvl,
 
     d = dvl[-1];
     for (b = MPW_BITS / 2; b; b >>= 1) {
-      if (d < (MPW_MAX >> b)) {
+      if (d <= (MPW_MAX >> b)) {
 	d <<= b;
 	norm += b;
       }