X-Git-Url: https://git.distorted.org.uk/u/mdw/catacomb/blobdiff_plain/d3409d5ecf2492cff862616de72a580d1a8e8dc0..d34decd2b2b88240cf4ca68a2a5feb7bf36de6e7:/mp-arith.c

diff --git a/mp-arith.c b/mp-arith.c
index d8381ee..fe6a45d 100644
--- a/mp-arith.c
+++ b/mp-arith.c
@@ -1,6 +1,6 @@
 /* -*-c-*-
  *
- * $Id: mp-arith.c,v 1.1 1999/11/17 18:02:16 mdw Exp $
+ * $Id: mp-arith.c,v 1.6 2000/06/17 11:45:09 mdw Exp $
  *
  * Basic arithmetic on multiprecision integers
  *
@@ -30,6 +30,23 @@
 /*----- Revision history --------------------------------------------------* 
  *
  * $Log: mp-arith.c,v $
+ * Revision 1.6  2000/06/17 11:45:09  mdw
+ * Major memory management overhaul.  Added arena support.  Use the secure
+ * arena for secret integers.  Replace and improve the MP management macros
+ * (e.g., replace MP_MODIFY by MP_DEST).
+ *
+ * Revision 1.5  1999/12/22 15:54:41  mdw
+ * Adjust Karatsuba parameters.  Calculate destination size better.
+ *
+ * Revision 1.4  1999/12/13 15:35:16  mdw
+ * Slightly different rules on memory allocation.
+ *
+ * Revision 1.3  1999/12/11 10:57:43  mdw
+ * Karatsuba squaring algorithm.
+ *
+ * Revision 1.2  1999/12/10 23:18:39  mdw
+ * Change interface for suggested destinations.
+ *
  * Revision 1.1  1999/11/17 18:02:16  mdw
  * New multiprecision integer arithmetic suite.
  *
@@ -39,6 +56,10 @@
 
 #include "mp.h"
 
+/*----- Macros ------------------------------------------------------------*/
+
+#define MAX(x, y) ((x) >= (y) ? (x) : (y))
+
 /*----- Main code ---------------------------------------------------------*/
 
 /* --- @mp_2c@ --- *
@@ -53,7 +74,7 @@ mp *mp_2c(mp *d, mp *a)
   if (!(a->f & MP_NEG))
     return (MP_COPY(a));
 
-  MP_MODIFY(d, MP_LEN(a));
+  MP_DEST(d, MP_LEN(a), a->f);
   mpx_2c(d->v, d->vl, a->v, a->vl);
   d->f = a->f & MP_BURN;
   MP_SHRINK(d);
@@ -74,7 +95,7 @@ mp *mp_sm(mp *d, mp *a)
   if (!MP_LEN(a) || a->vl[-1] < MPW_MAX / 2)
     return (MP_COPY(a));
 
-  MP_MODIFY(d, MP_LEN(a));
+  MP_DEST(d, MP_LEN(a), a->f);
   mpx_2c(d->v, d->vl, a->v, a->vl);
   d->f = (a->f & (MP_BURN | MP_NEG)) ^ MP_NEG;
   MP_SHRINK(d);
@@ -84,15 +105,15 @@ mp *mp_sm(mp *d, mp *a)
 /* --- @mp_lsl@ --- *
  *
  * Arguments:	@mp *d@ = destination
- *		@const mp *a@ = source
+ *		@mp *a@ = source
  *		@size_t n@ = number of bits to move
  *
  * Returns:	Result, @a@ shifted left by @n@.
  */
 
-mp *mp_lsl(mp *d, const mp *a, size_t n)
+mp *mp_lsl(mp *d, mp *a, size_t n)
 {
-  MP_MODIFY(d, MP_LEN(a) + (n + MPW_BITS - 1) / MPW_BITS);
+  MP_DEST(d, MP_LEN(a) + (n + MPW_BITS - 1) / MPW_BITS, a->f);
   mpx_lsl(d->v, d->vl, a->v, a->vl, n);
   d->f = a->f & (MP_NEG | MP_BURN);
   MP_SHRINK(d);
@@ -102,15 +123,15 @@ mp *mp_lsl(mp *d, const mp *a, size_t n)
 /* --- @mp_lsr@ --- *
  *
  * Arguments:	@mp *d@ = destination
- *		@const mp *a@ = source
+ *		@mp *a@ = source
  *		@size_t n@ = number of bits to move
  *
  * Returns:	Result, @a@ shifted left by @n@.
  */
 
-mp *mp_lsr(mp *d, const mp *a, size_t n)
+mp *mp_lsr(mp *d, mp *a, size_t n)
 {
-  MP_MODIFY(d, MP_LEN(a));
+  MP_DEST(d, MP_LEN(a), a->f);
   mpx_lsr(d->v, d->vl, a->v, a->vl, n);
   d->f = a->f & (MP_NEG | MP_BURN);
   MP_SHRINK(d);
@@ -138,19 +159,19 @@ int mp_cmp(const mp *a, const mp *b)
 /* --- @mp_add@ --- *
  *
  * Arguments:	@mp *d@ = destination
- *		@const mp *a, *b@ = sources
+ *		@mp *a, *b@ = sources
  *
  * Returns:	Result, @a@ added to @b@.
  */
 
-mp *mp_add(mp *d, const mp *a, const mp *b)
+mp *mp_add(mp *d, mp *a, mp *b)
 {
-  MP_MODIFY(d, (MP_LEN(a) > MP_LEN(b) ? MP_LEN(a) : MP_LEN(b)) + 1);
+  MP_DEST(d, MAX(MP_LEN(a), MP_LEN(b)) + 1, a->f | b->f);
   if (!((a->f ^ b->f) & MP_NEG))
     mpx_uadd(d->v, d->vl, a->v, a->vl, b->v, b->vl);
   else {
     if (MPX_UCMP(a->v, a->vl, <, b->v, b->vl)) {
-      const mp *t = a; a = b; b = t;
+      mp *t = a; a = b; b = t;
     }
     mpx_usub(d->v, d->vl, a->v, a->vl, b->v, b->vl);
   }
@@ -162,20 +183,20 @@ mp *mp_add(mp *d, const mp *a, const mp *b)
 /* --- @mp_sub@ --- *
  *
  * Arguments:	@mp *d@ = destination
- *		@const mp *a, *b@ = sources
+ *		@mp *a, *b@ = sources
  *
  * Returns:	Result, @b@ subtracted from @a@.
  */
 
-mp *mp_sub(mp *d, const mp *a, const mp *b)
+mp *mp_sub(mp *d, mp *a, mp *b)
 {
   unsigned sgn = 0;
-  MP_MODIFY(d, (MP_LEN(a) > MP_LEN(b) ? MP_LEN(a) : MP_LEN(b)) + 1);
+  MP_DEST(d, MAX(MP_LEN(a), MP_LEN(b)) + 1, a->f | b->f);
   if ((a->f ^ b->f) & MP_NEG)
     mpx_uadd(d->v, d->vl, a->v, a->vl, b->v, b->vl);
   else {
     if (MPX_UCMP(a->v, a->vl, <, b->v, b->vl)) {
-      const mp *t = a; a = b; b = t;
+      mp *t = a; a = b; b = t;
       sgn = MP_NEG;
     }
     mpx_usub(d->v, d->vl, a->v, a->vl, b->v, b->vl);
@@ -188,45 +209,68 @@ mp *mp_sub(mp *d, const mp *a, const mp *b)
 /* --- @mp_mul@ --- *
  *
  * Arguments:	@mp *d@ = destination
- *		@const mp *a, *b@ = sources
+ *		@mp *a, *b@ = sources
  *
  * Returns:	Result, @a@ multiplied by @b@.
  */
 
-mp *mp_mul(mp *d, const mp *a, const mp *b)
+mp *mp_mul(mp *d, mp *a, mp *b)
 {
-  if (d == a || d == b)
-    d = MP_NEW;
-  MP_MODIFY(d, MP_LEN(a) + MP_LEN(b));
-  mpx_umul(d->v, d->vl, a->v, a->vl, b->v, b->vl);
+  a = MP_COPY(a);
+  b = MP_COPY(b);
+
+  if (MP_LEN(a) <= KARATSUBA_CUTOFF || MP_LEN(b) <= KARATSUBA_CUTOFF) {
+    MP_DEST(d, MP_LEN(a) + MP_LEN(b), a->f | b->f | MP_UNDEF);
+    mpx_umul(d->v, d->vl, a->v, a->vl, b->v, b->vl);
+  } else {
+    size_t m = 2 * MAX(MP_LEN(a), MP_LEN(b)) + 2;
+    mpw *s;
+    MP_DEST(d, m, a->f | b->f | MP_UNDEF);
+    m += KARATSUBA_SLOP;
+    s = mpalloc(d->a, m);
+    mpx_kmul(d->v, d->vl, a->v, a->vl, b->v, b->vl, s, s + m);
+    mpfree(d->a, s);
+  }
+
   d->f = ((a->f | b->f) & MP_BURN) | ((a->f ^ b->f) & MP_NEG);
   MP_SHRINK(d);
+  MP_DROP(a);
+  MP_DROP(b);
   return (d);
 }
 
 /* --- @mp_sqr@ --- *
  *
  * Arguments:	@mp *d@ = destination
- *		@const mp *a@ = source
+ *		@mp *a@ = source
  *
  * Returns:	Result, @a@ squared.
  */
 
-mp *mp_sqr(mp *d, const mp *a)
+mp *mp_sqr(mp *d, mp *a)
 {
-  if (d == a)
-    d = MP_NEW;
-  MP_MODIFY(d, 2 * MP_LEN(a));
-  mpx_usqr(d->v, d->vl, a->v, a->vl);
+  size_t m = MP_LEN(a);
+
+  a = MP_COPY(a);
+  MP_DEST(d, 2 * m + 2, a->f | MP_UNDEF);
+  if (m > KARATSUBA_CUTOFF) {
+    mpw *s;
+    m = 2 * (m + 1) + KARATSUBA_SLOP;
+    s = mpalloc(d->a, m);
+    mpx_ksqr(d->v, d->vl, a->v, a->vl, s, s + m);
+    mpfree(d->a, s);
+  } else 
+    mpx_usqr(d->v, d->vl, a->v, a->vl);
   d->f = a->f & MP_BURN;
   MP_SHRINK(d);
+  MP_DROP(a);
   return (d);
 }
 
 /* --- @mp_div@ --- *
  *
  * Arguments:	@mp **qq, **rr@ = destination, quotient and remainder
- *		@const mp *a, *b@ = sources
+ *		@mp *a, *b@ = sources
  *
  * Use:		Calculates the quotient and remainder when @a@ is divided by
  *		@b@.  The destinations @*qq@ and @*rr@ must be distinct.
@@ -242,48 +286,38 @@ mp *mp_sqr(mp *d, const mp *a)
  *		straightforward.
  */
 
-void mp_div(mp **qq, mp **rr, const mp *a, const mp *b)
+void mp_div(mp **qq, mp **rr, mp *a, mp *b)
  {
   mp *r = rr ? *rr : MP_NEW;
   mp *q = qq ? *qq : MP_NEW;
   mpw *sv, *svl;
 
-  /* --- Set up some temporary workspace --- */
-
-  {
-    size_t rq = MP_LEN(b) + 1;
-    sv = MP_ALLOC(rq);
-    svl = sv + rq;
-  }
-
   /* --- Set the remainder up right --- *
    *
    * Just in case the divisor is larger, be able to cope with this.  It's not
    * important in @mpx_udiv@, but it is here because of the sign correction.
    */
 
-  {
-    size_t rq = MP_LEN(a) + 2;
-    if (MP_LEN(b) > rq)
-      rq = MP_LEN(b);
-
-    if (r == a) {
-      MP_SPLIT(r);
-      MP_ENSURE(r, MP_LEN(r) + 2);
-    } else {
-      if (r == b)
-	r = MP_NEW;
-      MP_MODIFY(r, MP_LEN(a) + 2);
-      memcpy(r->v, a->v, MPWS(MP_LEN(a)));
-      memset(r->v + MP_LEN(a), 0, MPWS(2));
-    }
-  }
+  b = MP_COPY(b);
+  a = MP_COPY(a);
+  if (r)
+    MP_DROP(r);
+  r = a;
+  MP_DEST(r, MP_LEN(a) + 2, a->f | b->f);
 
   /* --- Fix up the quotient too --- */
 
-  if (q == a || q == b)
-    q = MP_NEW;
-  MP_MODIFY(q, MP_LEN(a));
+  r = MP_COPY(r);
+  MP_DEST(q, MP_LEN(r), r->f | MP_UNDEF);
+  MP_DROP(r);
+
+  /* --- Set up some temporary workspace --- */
+
+  {
+    size_t rq = MP_LEN(b) + 1;
+    sv = mpalloc(r->a, rq);
+    svl = sv + rq;
+  }
 
   /* --- Perform the calculation --- */
 
@@ -296,10 +330,10 @@ void mp_div(mp **qq, mp **rr, const mp *a, const mp *b)
    * remainder from @b@.
    */
 
-  q->f = ((a->f | b->f) & MP_BURN) | ((a->f ^ b->f) & MP_NEG);
+  q->f = ((r->f | b->f) & MP_BURN) | ((r->f ^ b->f) & MP_NEG);
   if (q->f & MP_NEG) {
-    mpw *v = r->v;
-    while (v < r->vl) {
+    mpw *v;
+    for (v = r->v; v < r->vl; v++) {
       if (*v) {
 	MPX_UADDN(q->v, q->vl, 1);
 	mpx_usub(r->v, r->vl, b->v, b->vl, r->v, r->vl);
@@ -308,10 +342,13 @@ void mp_div(mp **qq, mp **rr, const mp *a, const mp *b)
     }
   }
 
-  r->f = ((a->f | b->f) & MP_BURN) | (b->f & MP_NEG);
+  r->f = ((r->f | b->f) & MP_BURN) | (b->f & MP_NEG);
 
   /* --- Store the return values --- */
 
+  mpfree(r->a, sv);
+  MP_DROP(b);
+
   if (!qq)
     MP_DROP(q);
   else {
@@ -325,8 +362,6 @@ void mp_div(mp **qq, mp **rr, const mp *a, const mp *b)
     MP_SHRINK(r);
     *rr = r;
   }
-
-  MP_FREE(sv);
 }
 
 /*----- Test rig ----------------------------------------------------------*/
@@ -348,7 +383,7 @@ static int verify(const char *op, mp *expect, mp *result, mp *a, mp *b)
 }
 
 #define RIG(name, op)							\
-  static int t ## name(dstr *v)						\
+  static int t##name(dstr *v)						\
   {									\
     mp *a = *(mp **)v[0].buf;						\
     mpw n = *(int *)v[1].buf;						\
@@ -359,6 +394,7 @@ static int verify(const char *op, mp *expect, mp *result, mp *a, mp *b)
     mp_build(&b, &n, &n + 1);						\
     ok = verify(#name, r, c, a, &b);					\
     mp_drop(a); mp_drop(c); mp_drop(r);					\
+    assert(mparena_count(MPARENA_GLOBAL) == 0);				\
     return (ok);							\
   }
 
@@ -368,7 +404,7 @@ RIG(lsr, mp_lsr)
 #undef RIG
 
 #define RIG(name, op)							\
-  static int t ## name(dstr *v)						\
+  static int t##name(dstr *v)						\
   {									\
     mp *a = *(mp **)v[0].buf;						\
     mp *b = *(mp **)v[1].buf;						\
@@ -376,6 +412,7 @@ RIG(lsr, mp_lsr)
     mp *c = op(MP_NEW, a, b);						\
     int ok = verify(#name, r, c, a, b);					\
     mp_drop(a); mp_drop(b); mp_drop(c); mp_drop(r);			\
+    assert(mparena_count(MPARENA_GLOBAL) == 0);				\
     return (ok);							\
   }
 
@@ -397,6 +434,7 @@ static int tdiv(dstr *v)
   ok &= verify("div(quotient)", q, c, a, b);
   ok &= verify("div(remainder)", r, d, a, b);
   mp_drop(a); mp_drop(b); mp_drop(c); mp_drop(d); mp_drop(r); mp_drop(q);
+  assert(mparena_count(MPARENA_GLOBAL) == 0);
   return (ok);
 }