From f45a00c61a25b1c4253bf28296fcb8d21f8eb7e1 Mon Sep 17 00:00:00 2001
From: mdw <mdw>
Date: Wed, 17 Nov 1999 18:04:09 +0000
Subject: [PATCH] Add two's-complement functionality.  Improve mpx_udiv a
 little by performing the multiplication of the divisor by q with the
 subtraction from r.

---
 mpx.c | 79 +++++++++++++++++++++++++++++++++++++++++++++++++++----------------
 1 file changed, 61 insertions(+), 18 deletions(-)

diff --git a/mpx.c b/mpx.c
index b6d765d..e3600b7 100644
--- a/mpx.c
+++ b/mpx.c
@@ -1,6 +1,6 @@
 /* -*-c-*-
  *
- * $Id: mpx.c,v 1.3 1999/11/13 01:57:31 mdw Exp $
+ * $Id: mpx.c,v 1.4 1999/11/17 18:04:09 mdw Exp $
  *
  * Low-level multiprecision arithmetic
  *
@@ -30,6 +30,11 @@
 /*----- Revision history --------------------------------------------------* 
  *
  * $Log: mpx.c,v $
+ * 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.
  *
@@ -379,6 +384,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
@@ -618,8 +647,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,
@@ -662,12 +690,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);
   }
 
@@ -753,25 +779,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--;
-- 
2.11.0