Fix various assumptions about mpw sizes.
authorMark Wooding <mdw@distorted.org.uk>
Tue, 4 Apr 2006 16:17:45 +0000 (17:17 +0100)
committerMark Wooding <mdw@distorted.org.uk>
Tue, 16 Jan 2007 22:20:06 +0000 (22:20 +0000)
  * configure, mptypes: New configure switches force mpw type to either
    sane but small (16/32 bits) or cussid (19/38 bits).  This found a
    bunch of exciting bugs...

  * gfreduce, mpreduce: If MPW_BITS is not a power of two, modular
    reduction of a `negative' unsigned value does the wrong thing.

  * mpx_lsl and friends: Shifting ops weren't masking high-order bits
    correctly when writing the output.  Apply MPW().

  * mpx_usubnlsl: More failure to elide high-order junk bits.

  * mptypes, mpx_udiv, mpx_bits, mp_odd: The binary search is neato, but
    starts in the wrong place if MPW_BITS is not a power of two.  Have
    mptypes  compute MPW_P2 as the largest power of two less than
    MPW_BITS.

configure.in
gfreduce.c
mp-arith.c
mpreduce.c
mptypes.c
mpx.c
mpx.h
tests/mpreduce
tests/mpx

index f35a8df..aef91e0 100644 (file)
@@ -84,6 +84,22 @@ dnl --- Support for the passphrase pixie ---
 mdw_CHECK_MANYLIBS(socket, socket)
 AC_CHECK_FUNCS(mlock)
 
+dnl --- Debugging things ---
+
+AC_ARG_ENABLE([mpw], 
+              [AS_HELP_STRING([--enable-mpw], 
+                              [force small-width mp digits])],
+              [case "$enableval" in
+                 y*|t*|short)
+                   AC_DEFINE([FORCE_MPW_SHORT], [1],
+                             [Define to force small-width mp digits.])
+                   ;;
+                 cussid)
+                   AC_DEFINE([FORCE_MPW_CUSSID], [1],
+                             [Define to force strange-width mp digits.])
+                   ;;
+               esac])
+
 dnl --- Done ---
 
 mdw_MLIB(2.0.0)
index bdf3579..01ca678 100644 (file)
@@ -90,7 +90,7 @@ void gfreduce_create(gfreduce *r, mp *p)
   unsigned long i;
   gfreduce_instr *ip;
   unsigned f = 0;
-  size_t w, ww, wi, wl, ll;
+  size_t w, ww, wi, wl, ll, bb;
 
   /* --- Sort out the easy stuff --- */
 
@@ -121,6 +121,7 @@ void gfreduce_create(gfreduce *r, mp *p)
   wi = DA_LEN(&iv);
   f = 0;
   ll = 0;
+  bb = MPW_BITS - dw;
   for (i = 0, mp_scan(&sc, p); mp_step(&sc) && i < d; i++) {
     if (!mp_bit(&sc))
       continue;
@@ -149,8 +150,8 @@ void gfreduce_create(gfreduce *r, mp *p)
       w = ww;
       wi = DA_LEN(&iv);
     }
-    INSTR(GFRI_LSL, (MPW_BITS + i - d)%MPW_BITS);
-    if ((MPW_BITS + i - d)%MPW_BITS)
+    INSTR(GFRI_LSL, (bb + i)%MPW_BITS);
+    if ((bb + i)%MPW_BITS)
       f |= f_lsr;
   }
   wl = DA_LEN(&iv);
index 3af1447..9cb5178 100644 (file)
@@ -652,16 +652,16 @@ mp *mp_odd(mp *d, mp *m, size_t *s)
     ss = 0;
   else {
     mpw x = *v;
-    mpw mask = MPW_MAX;
-    unsigned z = MPW_BITS / 2;
+    unsigned z = MPW_P2;
+    mpw mask = ((mpw)1 << z) - 1;
 
     while (z) {
-      mask >>= z;
       if (!(x & mask)) {
        x >>= z;
        ss += z;
       }
       z >>= 1;
+      mask >>= z;
     }
   }
 
index 1c39564..bc41f60 100644 (file)
@@ -60,7 +60,7 @@ int mpreduce_create(mpreduce *r, mp *p)
   instr_v iv = DA_INIT;
   unsigned long d, i;
   unsigned op;
-  size_t w, b;
+  size_t w, b, bb;
 
   /* --- Fill in the easy stuff --- */
 
@@ -105,6 +105,7 @@ int mpreduce_create(mpreduce *r, mp *p)
   st = Z;
 #endif
 
+  bb = MPW_BITS - (d + 1)%MPW_BITS;
   for (i = 0, mp_scan(&sc, p); i < d && mp_step(&sc); i++) {
     switch (st | mp_bit(&sc)) {
       case  Z | 1: st = Z1; break;
@@ -115,7 +116,7 @@ int mpreduce_create(mpreduce *r, mp *p)
       case X0 | 0: st =  Z; op = MPRI_SUB; goto instr;
       instr:
        w = (d - i)/MPW_BITS + 1;
-       b = (MPW_BITS + i - d - 1)%MPW_BITS;
+       b = (bb + i)%MPW_BITS;
        INSTR(op | !!b, w, b);
     }
   }
@@ -154,6 +155,10 @@ int mpreduce_create(mpreduce *r, mp *p)
     }
   }
   DA_DESTROY(&iv);
+
+#ifdef DEBUG
+  mpreduce_dump(r, stdout);
+#endif
   return (0);
 }
 
index 613e214..381f3dc 100644 (file)
--- a/mptypes.c
+++ b/mptypes.c
@@ -30,6 +30,8 @@
 /*----- Header files ------------------------------------------------------*/
 
 #define _GNU_SOURCE
+#include "config.h"
+
 #include <stdio.h>
 #include <limits.h>
 #if __STDC_VERSION__ >= 199900l
@@ -105,6 +107,7 @@ int main(int argc, char *argv[])
   itype *i;
   itype *largest, *mpw, *mpd;
   const static char *extstr = "CATACOMB_MPTYPES_EXTENSION ";
+  unsigned p2;
 
   /* --- Find the bitcounts --- */
 
@@ -123,6 +126,17 @@ int main(int argc, char *argv[])
    * which is twice as big as that one.
    */
 
+#if defined(FORCE_MPW_CUSSID)
+  largest = mpd = &tytab[3];
+  mpw = &tytab[2];
+  mpw->bits = 19; mpw->max = 0x7ffff;
+  mpd->bits = 38; mpd->max = 0x3fffffffffll;
+#elif defined(FORCE_MPW_SHORT)
+  largest = mpd = &tytab[2];
+  mpw = &tytab[1];
+  mpw->bits = 16; mpw->max = 0xffff;
+  mpd->bits = 32; mpd->max = 0xffffffff;
+#else
   largest = tytab;
   for (i = tytab; i->name; i++) {
     if (i->bits > largest->bits)
@@ -145,6 +159,8 @@ int main(int argc, char *argv[])
     d.bits = w.bits * 2; d.max = ~(~((umax)0) << d.bits);
     mpw = &w; mpd = &d;
   }    
+#endif
+  for (p2 = 1; (p2 << 1) < mpw->bits; p2 <<= 1);
 
   /* --- Output time --- */
 
@@ -176,6 +192,7 @@ int main(int argc, char *argv[])
   printf("\
 %stypedef %s mpw;\n\
 #define MPW_BITS %u\n\
+#define MPW_P2 %u\n\
 #define MPW_MAX %s%" P_UMAX "%s\n\
 \n\
 %stypedef %s mpd;\n\
@@ -185,7 +202,7 @@ int main(int argc, char *argv[])
 #endif\n\
 ",
   mpw->flags & f_ext ? extstr : "", mpw->name,
-  mpw->bits,
+  mpw->bits, p2,
   mpw->flags & f_ext ? extstr : "", mpw->max, mpw->suff,
   mpd->flags & f_ext ? extstr : "", mpd->name,
   mpd->bits,
diff --git a/mpx.c b/mpx.c
index c6a23be..01264b1 100644 (file)
--- a/mpx.c
+++ b/mpx.c
@@ -462,11 +462,11 @@ void mpx_lsl(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, size_t n)
 
     while (avl > av) {
       mpw t = *--avl;
-      *--dvl = (t >> nr) | w;
+      *--dvl = MPW((t >> nr) | w);
       w = t << nb;
     }
 
-    *--dvl = w;
+    *--dvl = MPW(w);
     MPX_ZERO(dv, dvl);
   }
 
@@ -560,11 +560,11 @@ void mpx_lslc(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, size_t n)
 
     while (avl > av) {
       mpw t = *--avl;
-      *--dvl = (t >> nr) | w;
+      *--dvl = MPW((t >> nr) | w);
       w = t << nb;
     }
 
-    *--dvl = (MPW_MAX >> nr) | w;
+    *--dvl = MPW((MPW_MAX >> nr) | w);
     MPX_ONE(dv, dvl);
   }
 
@@ -919,7 +919,7 @@ void mpx_usubnlsl(mpw *dv, mpw *dvl, mpw a, unsigned o)
   a <<= o;
 
   if (dv < dvl) {
-    mpd x = (mpd)*dv - (mpd)a;
+    mpd x = (mpd)*dv - MPW(a);
     *dv++ = MPW(x);
     if (x >> MPW_BITS)
       b++;
@@ -1116,7 +1116,7 @@ void mpx_udiv(mpw *qv, mpw *qvl, mpw *rv, mpw *rvl,
     unsigned b;
 
     d = dvl[-1];
-    for (b = MPW_BITS / 2; b; b >>= 1) {
+    for (b = MPW_P2; b; b >>= 1) {
       if (d <= (MPW_MAX >> b)) {
        d <<= b;
        norm += b;
diff --git a/mpx.h b/mpx.h
index f79cffd..15253ed 100644 (file)
--- a/mpx.h
+++ b/mpx.h
@@ -89,7 +89,7 @@
   else {                                                               \
     unsigned long _b = MPW_BITS * (_vl - _v - 1) + 1;                  \
     mpw _w = _vl[-1];                                                  \
-    unsigned _k = MPW_BITS / 2;                                                \
+    unsigned _k = MPW_P2;                                              \
     while (_k) {                                                       \
       if (_w >> _k) {                                                  \
        _w >>= _k;                                                      \
index faa84d2..1216f2e 100644 (file)
@@ -6,6 +6,12 @@ reduce {
   0xc000 0x16cb3 0xacb3;
   0x8000 0x345545 0x5545;
 
+  0xfffef 0x100000 0x11;
+
+  0x1ffffffe 0x26fc6567 0x6fc6569;
+  0x3ffffffe 0x45445dc0 0x5445dc2;
+  0x7ffffffe 0xd4827a70 0x54827a72;
+
   0x72e2c37447f8bca34c4a39b130ea8e5c9a7d8b54564aa88ea773
   0x367aa8f5ba9ac4e8e2ea198b8af2c3b3081deab392ffc05715783b245a62a6fa
   0x08e8c03ebf398c63d71d8fd7ca4ece12367a8dde180ca650afb6;
index 7f6ff5e..1ebc869 100644 (file)
--- a/tests/mpx
+++ b/tests/mpx
@@ -885,4 +885,6 @@ udiv {
   ffffffffffffffffc90fdaa22168c234c4c6628b80dc1cd129024e088a67cc74020bbea63b139b22514a08798e3404ddef9519b3cd3a431b302b0a6df25f14374fe1356d6d51c245e485b576625e7ec6f44c42e9a63a3620ffffffffffffffff
   7fffffffffffffffe487ed5110b4611a62633145c06e0e68948127044533e63a0105df531d89cd9128a5043cc71a026ef7ca8cd9e69d218d98158536f92f8a1ba7f09ab6b6a8e122f242dabb312f3f637a262174d31d1b107fffffffffffffff
   02 01;
+
+  26737e 0ffffc 02 067386;
 }