.gitignore: Ignore `auto-version' script.
[u/mdw/catacomb] / mpx.c
diff --git a/mpx.c b/mpx.c
index ef93e3e..29e1ded 100644 (file)
--- a/mpx.c
+++ b/mpx.c
@@ -1,13 +1,13 @@
 /* -*-c-*-
  *
- * $Id: mpx.c,v 1.19 2004/04/03 03:29:40 mdw Exp $
+ * $Id$
  *
  * Low-level multiprecision arithmetic
  *
  * (c) 1999 Straylight/Edgeware
  */
 
-/*----- Licensing notice --------------------------------------------------* 
+/*----- Licensing notice --------------------------------------------------*
  *
  * This file is part of Catacomb.
  *
  * it under the terms of the GNU Library General Public License as
  * published by the Free Software Foundation; either version 2 of the
  * License, or (at your option) any later version.
- * 
+ *
  * Catacomb is distributed in the hope that it will be useful,
  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  * GNU Library General Public License for more details.
- * 
+ *
  * You should have received a copy of the GNU Library General Public
  * License along with Catacomb; if not, write to the Free
  * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
  * MA 02111-1307, USA.
  */
 
-/*----- Revision history --------------------------------------------------* 
- *
- * $Log: mpx.c,v $
- * Revision 1.19  2004/04/03 03:29:40  mdw
- * Fix overrun in @mpx_lsr@.
- *
- * Revision 1.18  2004/04/01 12:50:09  mdw
- * Add cyclic group abstraction, with test code.  Separate off exponentation
- * functions for better static linking.  Fix a buttload of bugs on the way.
- * Generally ensure that negative exponents do inversion correctly.  Add
- * table of standard prime-field subgroups.  (Binary field subgroups are
- * currently unimplemented but easy to add if anyone ever finds a good one.)
- *
- * Revision 1.17  2004/03/27 00:04:46  mdw
- * Implement efficient reduction for pleasant-looking primes.
- *
- * 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.
- *
- * Revision 1.14  2002/10/19 18:55:08  mdw
- * Fix overflows in shift primitives.
- *
- * Revision 1.13  2002/10/19 17:56:50  mdw
- * Fix bit operations.  Test them (a bit) better.
- *
- * Revision 1.12  2002/10/06 22:52:50  mdw
- * Pile of changes for supporting two's complement properly.
- *
- * Revision 1.11  2001/04/03 19:36:05  mdw
- * Add some simple bitwise operations so that Perl can use them.
- *
- * Revision 1.10  2000/10/08 12:06:12  mdw
- * Provide @mpx_ueq@ for rapidly testing equality of two integers.
- *
- * Revision 1.9  2000/06/26 07:52:50  mdw
- * Portability fix for the bug fix.
- *
- * Revision 1.8  2000/06/25 12:59:02  mdw
- * (mpx_udiv): Fix bug in quotient digit estimation.
- *
- * Revision 1.7  1999/12/22 15:49:07  mdw
- * New function for division by a small integer.
- *
- * Revision 1.6  1999/11/20 22:43:44  mdw
- * Integrate testing for MPX routines.
- *
- * 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.
- *
- * Revision 1.1  1999/09/03 08:41:12  mdw
- * Initial import.
- *
- */
-
 /*----- Header files ------------------------------------------------------*/
 
 #include <assert.h>
@@ -485,7 +417,7 @@ void mpx_lsl(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, size_t n)
   }
 
   /* --- Break out word and bit shifts for more sophisticated work --- */
-       
+
   nw = n / MPW_BITS;
   nb = n % MPW_BITS;
 
@@ -530,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);
   }
 
@@ -583,7 +515,7 @@ void mpx_lslc(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, size_t n)
   }
 
   /* --- Break out word and bit shifts for more sophisticated work --- */
-       
+
   nw = n / MPW_BITS;
   nb = n % MPW_BITS;
 
@@ -628,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);
   }
 
@@ -987,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++;
@@ -1018,7 +950,7 @@ void mpx_umul(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl,
   MPX_SHRINK(bv, bvl);
 
   /* --- Deal with a multiply by zero --- */
-  
+
   if (bv == bvl) {
     MPX_ZERO(dv, dvl);
     return;
@@ -1063,9 +995,7 @@ void mpx_umul(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl,
  */
 
 void mpx_umuln(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, mpw m)
-{
-  MPX_UMULN(dv, dvl, av, avl, m);
-}
+  { MPX_UMULN(dv, dvl, av, avl, m); }
 
 /* --- @mpx_umlan@ --- *
  *
@@ -1080,9 +1010,7 @@ void mpx_umuln(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, mpw m)
  */
 
 void mpx_umlan(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, mpw m)
-{
-  MPX_UMLAN(dv, dvl, av, avl, m);
-}
+  { MPX_UMLAN(dv, dvl, av, avl, m); }
 
 /* --- @mpx_usqr@ --- *
  *
@@ -1188,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;
@@ -1402,7 +1330,7 @@ mpw mpx_udivn(mpw *qv, mpw *qvl, const mpw *rv, const mpw *rvl, mpw d)
 } while (0)
 
 #define MAX(x, y) ((x) > (y) ? (x) : (y))
-  
+
 static void dumpbits(const char *msg, const void *pp, size_t sz)
 {
   const octet *p = pp;
@@ -1489,7 +1417,7 @@ static int loadstore(dstr *v)
   if (!ok)
     dumpbits("input data", v->buf, v->len);
 
-  free(m);
+  xfree(m);
   dstr_destroy(&d);
   return (ok);
 }
@@ -1527,7 +1455,7 @@ static int twocl(dstr *v)
     dumpbits("neg", v[1].buf, v[1].len);
   }
 
-  free(m);
+  xfree(m);
   dstr_destroy(&d);
 
   return (ok);
@@ -1566,7 +1494,7 @@ static int twocb(dstr *v)
     dumpbits("neg", v[1].buf, v[1].len);
   }
 
-  free(m);
+  xfree(m);
   dstr_destroy(&d);
 
   return (ok);
@@ -1587,13 +1515,13 @@ static int lsl(dstr *v)
   mpx_lsl(d, dl, a, al, n);
   if (!mpx_ueq(d, dl, c, cl)) {
     fprintf(stderr, "\n*** lsl(%i) failed\n", n);
-    dumpmp("       a", a, al);
+    dumpmp("      a", a, al);
     dumpmp("expected", c, cl);
     dumpmp("  result", d, dl);
     ok = 0;
   }
 
-  free(a); free(c); free(d);
+  xfree(a); xfree(c); xfree(d);
   return (ok);
 }
 
@@ -1612,13 +1540,13 @@ static int lslc(dstr *v)
   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("      a", a, al);
     dumpmp("expected", c, cl);
     dumpmp("  result", d, dl);
     ok = 0;
   }
 
-  free(a); free(c); free(d);
+  xfree(a); xfree(c); xfree(d);
   return (ok);
 }
 
@@ -1637,13 +1565,13 @@ static int lsr(dstr *v)
   mpx_lsr(d, dl, a, al, n);
   if (!mpx_ueq(d, dl, c, cl)) {
     fprintf(stderr, "\n*** lsr(%i) failed\n", n);
-    dumpmp("       a", a, al);
+    dumpmp("      a", a, al);
     dumpmp("expected", c, cl);
     dumpmp("  result", d, dl);
     ok = 0;
   }
 
-  free(a); free(c); free(d);
+  xfree(a); xfree(c); xfree(d);
   return (ok);
 }
 
@@ -1663,14 +1591,14 @@ static int uadd(dstr *v)
   mpx_uadd(d, dl, a, al, b, bl);
   if (!mpx_ueq(d, dl, c, cl)) {
     fprintf(stderr, "\n*** uadd failed\n");
-    dumpmp("       a", a, al);
-    dumpmp("       b", b, bl);
+    dumpmp("      a", a, al);
+    dumpmp("      b", b, bl);
     dumpmp("expected", c, cl);
     dumpmp("  result", d, dl);
     ok = 0;
   }
 
-  free(a); free(b); free(c); free(d);
+  xfree(a); xfree(b); xfree(c); xfree(d);
   return (ok);
 }
 
@@ -1690,14 +1618,14 @@ static int usub(dstr *v)
   mpx_usub(d, dl, a, al, b, bl);
   if (!mpx_ueq(d, dl, c, cl)) {
     fprintf(stderr, "\n*** usub failed\n");
-    dumpmp("       a", a, al);
-    dumpmp("       b", b, bl);
+    dumpmp("      a", a, al);
+    dumpmp("      b", b, bl);
     dumpmp("expected", c, cl);
     dumpmp("  result", d, dl);
     ok = 0;
   }
 
-  free(a); free(b); free(c); free(d);
+  xfree(a); xfree(b); xfree(c); xfree(d);
   return (ok);
 }
 
@@ -1717,14 +1645,14 @@ static int umul(dstr *v)
   mpx_umul(d, dl, a, al, b, bl);
   if (!mpx_ueq(d, dl, c, cl)) {
     fprintf(stderr, "\n*** umul failed\n");
-    dumpmp("       a", a, al);
-    dumpmp("       b", b, bl);
+    dumpmp("      a", a, al);
+    dumpmp("      b", b, bl);
     dumpmp("expected", c, cl);
     dumpmp("  result", d, dl);
     ok = 0;
   }
 
-  free(a); free(b); free(c); free(d);
+  xfree(a); xfree(b); xfree(c); xfree(d);
   return (ok);
 }
 
@@ -1742,13 +1670,13 @@ static int usqr(dstr *v)
   mpx_usqr(d, dl, a, al);
   if (!mpx_ueq(d, dl, c, cl)) {
     fprintf(stderr, "\n*** usqr failed\n");
-    dumpmp("       a", a, al);
+    dumpmp("      a", a, al);
     dumpmp("expected", c, cl);
     dumpmp("  result", d, dl);
     ok = 0;
   }
 
-  free(a); free(c); free(d);
+  xfree(a); xfree(c); xfree(d);
   return (ok);
 }
 
@@ -1781,7 +1709,7 @@ static int udiv(dstr *v)
     ok = 0;
   }
 
-  free(a); free(b); free(r); free(q); free(s); free(qq);
+  xfree(a); xfree(b); xfree(r); xfree(q); xfree(s); xfree(qq);
   return (ok);
 }