Add two's-complement functionality. Improve mpx_udiv a little by

performing the multiplication of the divisor by q with the subtraction
from r.

diff --git a/mpx.c b/mpx.c
index b6d765d..e3600b7 100644 (file)
--- 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--;