Fix @mp_lsl2c@. Turns out to be surprisingly tricky.

[u/mdw/catacomb] / mpx.c

diff --git a/mpx.c b/mpx.c
index 6375c3e..f40d7d9 100644 (file)
--- a/mpx.c
+++ b/mpx.c
@@ -1,6 +1,6 @@
 /* -*-c-*-
  *
- * $Id: mpx.c,v 1.15 2002/10/20 01:12:31 mdw Exp $
+ * $Id: mpx.c,v 1.16 2003/05/16 09:09:24 mdw Exp $
  *
  * Low-level multiprecision arithmetic
  *
@@ -30,6 +30,9 @@
 /*----- Revision history --------------------------------------------------* 
  *
  * $Log: mpx.c,v $
+ * Revision 1.16  2003/05/16 09:09:24  mdw
+ * Fix @mp_lsl2c@.  Turns out to be surprisingly tricky.
+ *
  * Revision 1.15  2002/10/20 01:12:31  mdw
  * Two's complement I/O fixes.
  *
@@ -525,6 +528,104 @@ void mpx_lsl(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, size_t n)
 done:;
 }
 
+/* --- @mpx_lslc@ --- *
+ *
+ * Arguments:  @mpw *dv, *dvl@ = destination vector base and limit
+ *             @const mpw *av, *avl@ = source vector base and limit
+ *             @size_t n@ = number of bit positions to shift by
+ *
+ * Returns:    ---
+ *
+ * Use:                Performs a logical shift left operation on an integer, only
+ *             it fills in the bits with ones instead of zeroes.
+ */
+
+void mpx_lslc(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, size_t n)
+{
+  size_t nw;
+  unsigned nb;
+
+  /* --- Trivial special case --- */
+
+  if (n == 0)
+    MPX_COPY(dv, dvl, av, avl);
+
+  /* --- Single bit shifting --- */
+
+  else if (n == 1) {
+    mpw w = 1;
+    while (av < avl) {
+      mpw t;
+      if (dv >= dvl)
+       goto done;
+      t = *av++;
+      *dv++ = MPW((t << 1) | w);
+      w = t >> (MPW_BITS - 1);
+    }
+    if (dv >= dvl)
+      goto done;
+    *dv++ = MPW(w);
+    MPX_ZERO(dv, dvl);
+    goto done;
+  }
+
+  /* --- Break out word and bit shifts for more sophisticated work --- */
+       
+  nw = n / MPW_BITS;
+  nb = n % MPW_BITS;
+
+  /* --- Handle a shift by a multiple of the word size --- */
+
+  if (nb == 0) {
+    if (nw >= dvl - dv)
+      MPX_ONE(dv, dvl);
+    else {
+      MPX_COPY(dv + nw, dvl, av, avl);
+      MPX_ONE(dv, dv + nw);
+    }
+  }
+
+  /* --- And finally the difficult case --- *
+   *
+   * This is a little convoluted, because I have to start from the end and
+   * work backwards to avoid overwriting the source, if they're both the same
+   * block of memory.
+   */
+
+  else {
+    mpw w;
+    size_t nr = MPW_BITS - nb;
+    size_t dvn = dvl - dv;
+    size_t avn = avl - av;
+
+    if (dvn <= nw) {
+      MPX_ONE(dv, dvl);
+      goto done;
+    }
+
+    if (dvn > avn + nw) {
+      size_t off = avn + nw + 1;
+      MPX_ZERO(dv + off, dvl);
+      dvl = dv + off;
+      w = 0;
+    } else {
+      avl = av + dvn - nw;
+      w = *--avl << nb;
+    }
+
+    while (avl > av) {
+      mpw t = *--avl;
+      *--dvl = (t >> nr) | w;
+      w = t << nb;
+    }
+
+    *--dvl = (MPW_MAX >> nr) | w;
+    MPX_ONE(dv, dvl);
+  }
+
+done:;
+}
+
 /* --- @mpx_lsr@ --- *
  *
  * Arguments:  @mpw *dv, *dvl@ = destination vector base and limit
@@ -1432,6 +1533,31 @@ static int lsl(dstr *v)
   return (ok);
 }
 
+static int lslc(dstr *v)
+{
+  mpw *a, *al;
+  int n = *(int *)v[1].buf;
+  mpw *c, *cl;
+  mpw *d, *dl;
+  int ok = 1;
+
+  LOAD(a, al, &v[0]);
+  LOAD(c, cl, &v[2]);
+  ALLOC(d, dl, al - a + (n + MPW_BITS - 1) / MPW_BITS);
+
+  mpx_lslc(d, dl, a, al, n);
+  if (!mpx_ueq(d, dl, c, cl)) {
+    fprintf(stderr, "\n*** lslc(%i) failed\n", n);
+    dumpmp("       a", a, al);
+    dumpmp("expected", c, cl);
+    dumpmp("  result", d, dl);
+    ok = 0;
+  }
+
+  free(a); free(c); free(d);
+  return (ok);
+}
+
 static int lsr(dstr *v)
 {
   mpw *a, *al;
@@ -1600,6 +1726,7 @@ static test_chunk defs[] = {
   { "2cl", twocl, { &type_hex, &type_hex, } },
   { "2cb", twocb, { &type_hex, &type_hex, } },
   { "lsl", lsl, { &type_hex, &type_int, &type_hex, 0 } },
+  { "lslc", lslc, { &type_hex, &type_int, &type_hex, 0 } },
   { "lsr", lsr, { &type_hex, &type_int, &type_hex, 0 } },
   { "uadd", uadd, { &type_hex, &type_hex, &type_hex, 0 } },
   { "usub", usub, { &type_hex, &type_hex, &type_hex, 0 } },