X-Git-Url: https://git.distorted.org.uk/u/mdw/catacomb/blobdiff_plain/c8a2f9ef9666e435eb08f46b4059d18182ccd3bc..44c240eef3706e9b1a0451fb4d584ba7caea0cdf:/mpx.c

diff --git a/mpx.c b/mpx.c
index b8bc8bf..80628b8 100644
--- a/mpx.c
+++ b/mpx.c
@@ -1,6 +1,6 @@
 /* -*-c-*-
  *
- * $Id: mpx.c,v 1.2 1999/11/13 01:50:59 mdw Exp $
+ * $Id: mpx.c,v 1.5 1999/11/20 22:23:27 mdw Exp $
  *
  * Low-level multiprecision arithmetic
  *
@@ -30,6 +30,17 @@
 /*----- Revision history --------------------------------------------------* 
  *
  * $Log: mpx.c,v $
+ * Revision 1.5  1999/11/20 22:23:27  mdw
+ * Add function versions of some low-level macros with wider use.
+ *
+ * Revision 1.4  1999/11/17 18:04:09  mdw
+ * Add two's-complement functionality.  Improve mpx_udiv a little by
+ * performing the multiplication of the divisor by q with the subtraction
+ * from r.
+ *
+ * Revision 1.3  1999/11/13 01:57:31  mdw
+ * Remove stray debugging code.
+ *
  * Revision 1.2  1999/11/13 01:50:59  mdw
  * Multiprecision routines finished and tested.
  *
@@ -376,6 +387,30 @@ done:;
 
 /*----- Unsigned arithmetic -----------------------------------------------*/
 
+/* --- @mpx_2c@ --- *
+ *
+ * Arguments:	@mpw *dv, *dvl@ = destination vector
+ *		@const mpw *v, *vl@ = source vector
+ *
+ * Returns:	---
+ *
+ * Use:		Calculates the two's complement of @v@.
+ */
+
+void mpx_2c(mpw *dv, mpw *dvl, const mpw *v, const mpw *vl)
+{
+  mpw c = 0;
+  while (dv < dvl && v < vl)
+    *dv++ = c = MPW(~*v++);
+  if (dv < dvl) {
+    if (c > MPW_MAX / 2)
+      c = MPW(~0);
+    while (dv < dvl)
+      *dv++ = c;
+  }
+  MPX_UADDN(dv, dvl, 1);
+}
+
 /* --- @mpx_ucmp@ --- *
  *
  * Arguments:	@const mpw *av, *avl@ = first argument vector base and limit
@@ -445,6 +480,18 @@ void mpx_uadd(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl,
   }
 }
 
+/* --- @mpx_uaddn@ --- *
+ *
+ * Arguments:	@mpw *dv, *dvl@ = source and destination base and limit
+ *		@mpw n@ = other addend
+ *
+ * Returns:	---
+ *
+ * Use:		Adds a small integer to a multiprecision number.
+ */
+
+void mpx_uaddn(mpw *dv, mpw *dvl, mpw n) { MPX_UADDN(dv, dvl, n); }
+
 /* --- @mpx_usub@ --- *
  *
  * Arguments:	@mpw *dv, *dvl@ = destination vector base and limit
@@ -487,6 +534,18 @@ void mpx_usub(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl,
     *dv++ = c;
 }
 
+/* --- @mpx_usubn@ --- *
+ *
+ * Arguments:	@mpw *dv, *dvl@ = source and destination base and limit
+ *		@n@ = subtrahend
+ *
+ * Returns:	---
+ *
+ * Use:		Subtracts a small integer from a multiprecision number.
+ */
+
+void mpx_usubn(mpw *dv, mpw *dvl, mpw n) { MPX_USUBN(dv, dvl, n); }
+
 /* --- @mpx_umul@ --- *
  *
  * Arguments:	@mpw *dv, *dvl@ = destination vector base and limit
@@ -541,6 +600,41 @@ void mpx_umul(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl,
   }
 }
 
+/* --- @mpx_umuln@ --- *
+ *
+ * Arguments:	@mpw *dv, *dvl@ = destination vector base and limit
+ *		@const mpw *av, *avl@ = multiplicand vector base and limit
+ *		@mpw m@ = multiplier
+ *
+ * Returns:	---
+ *
+ * Use:		Multiplies a multiprecision integer by a single-word value.
+ *		The destination and source may be equal.  The destination
+ *		is completely cleared after use.
+ */
+
+void mpx_umuln(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, mpw m)
+{
+  MPX_UMULN(dv, dvl, av, avl, m);
+}
+
+/* --- @mpx_umlan@ --- *
+ *
+ * Arguments:	@mpw *dv, *dvl@ = destination/accumulator base and limit
+ *		@const mpw *av, *avl@ = multiplicand vector base and limit
+ *		@mpw m@ = multiplier
+ *
+ * Returns:	---
+ *
+ * Use:		Multiplies a multiprecision integer by a single-word value
+ *		and adds the result to an accumulator.
+ */
+
+void mpx_umlan(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, mpw m)
+{
+  MPX_UMLAN(dv, dvl, av, avl, m);
+}
+
 /* --- @mpx_usqr@ --- *
  *
  * Arguments:	@mpw *dv, *dvl@ = destination vector base and limit
@@ -572,7 +666,6 @@ void mpx_usqr(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl)
     /* --- Work out the square at this point in the proceedings --- */
 
     {
-      mpw d = *dvv;
       mpd x = (mpd)a * (mpd)a + *dvv;
       *dvv++ = MPW(x);
       c = MPW(x >> MPW_BITS);
@@ -582,7 +675,6 @@ void mpx_usqr(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl)
 
     avv++;
     while (dvv < dvl && avv < avl) {
-      mpw aa = *avv;
       mpd x = (mpd)a * (mpd)*avv++;
       mpd y = ((x << 1) & MPW_MAX) + c + *dvv;
       c = (x >> (MPW_BITS - 1)) + (y >> MPW_BITS);
@@ -617,8 +709,7 @@ void mpx_usqr(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl)
  *		requiring the dividend to be in the result position but it
  *		does make some sense really.  The remainder must have
  *		headroom for at least two extra words.  The scratch space
- *		must be at least two words larger than twice the size of the
- *		divisor.
+ *		must be at least one word larger than the divisor.
  */
 
 void mpx_udiv(mpw *qv, mpw *qvl, mpw *rv, mpw *rvl,
@@ -661,12 +752,10 @@ void mpx_udiv(mpw *qv, mpw *qvl, mpw *rv, mpw *rvl,
   /* --- Normalize the dividend/remainder to match --- */
 
   if (norm) {
-    mpw *svvl = sv + (dvl - dv) + 1;
     mpx_lsl(rv, rvl, rv, rvl, norm);
-    mpx_lsl(sv, svvl, dv, dvl, norm);
+    mpx_lsl(sv, svl, dv, dvl, norm);
     dv = sv;
-    sv = svvl;
-    dvl = svvl;
+    dvl = svl;
     MPX_SHRINK(dv, dvl);
   }
 
@@ -752,25 +841,42 @@ void mpx_udiv(mpw *qv, mpw *qvl, mpw *rv, mpw *rvl,
       {
 	mpw *svv;
 	const mpw *dvv;
-	mpw c = 0;
+	mpw mc = 0, sc = 0;
 
-	/* --- Calculate the size of the chunk --- */
+	/* --- Calculate the size of the chunk --- *
+	 *
+	 * This does the whole job of calculating @r >> scale - qd@.
+	 */
 
-	for (svv = sv, dvv = dv; dvv < dvl; svv++, dvv++) {
-	  mpd x = (mpd)*dvv * (mpd)q + c;
+	for (svv = rv + scale, dvv = dv;
+	     dvv < dvl && svv < rvl;
+	     svv++, dvv++) {
+	  mpd x = (mpd)*dvv * (mpd)q + mc;
+	  mc = x >> MPW_BITS;
+	  x = (mpd)*svv - MPW(x) - sc;
 	  *svv = MPW(x);
-	  c = x >> MPW_BITS;
+	  if (x >> MPW_BITS)
+	    sc = 1;
+	  else
+	    sc = 0;
+	}
+
+	if (svv < rvl) {
+	  mpd x = (mpd)*svv - mc - sc;
+	  *svv++ = MPW(x);
+	  if (x >> MPW_BITS)
+	    sc = MPW_MAX;
+	  else
+	    sc = 0;
+	  while (svv < rvl)
+	    *svv++ = sc;
 	}
-	if (c)
-	  *svv++ = c;
 
-	/* --- Now make sure that we can cope with the difference --- *
+	/* --- Fix if the quotient was too large --- *
 	 *
-	 * Take advantage of the fact that subtraction works two's-
-	 * complement.
+	 * This doesn't seem to happen very often.
 	 */
 
-	mpx_usub(rv + scale, rvl, rv + scale, rvl, sv, svv);
 	if (rvl[-1] > MPW_MAX / 2) {
 	  mpx_uadd(rv + scale, rvl, rv + scale, rvl, dv, dvl);
 	  q--;