From: mdw <mdw>
Date: Fri, 10 Dec 1999 23:18:39 +0000 (+0000)
Subject: Change interface for suggested destinations.
X-Git-Url: https://git.distorted.org.uk/u/mdw/catacomb/commitdiff_plain/ef5f48103e83977bda6ef4d7d1aacbb66a629b10

Change interface for suggested destinations.
---

diff --git a/dh-prime.c b/dh-prime.c
index 81dc423..ebe3a33 100644
--- a/dh-prime.c
+++ b/dh-prime.c
@@ -1,6 +1,6 @@
 /* -*-c-*-
  *
- * $Id: dh-prime.c,v 1.1 1999/11/20 22:24:44 mdw Exp $
+ * $Id: dh-prime.c,v 1.2 1999/12/10 23:18:38 mdw Exp $
  *
  * Generate (safe) Diffie-Hellman primes
  *
@@ -30,6 +30,9 @@
 /*----- Revision history --------------------------------------------------* 
  *
  * $Log: dh-prime.c,v $
+ * Revision 1.2  1999/12/10 23:18:38  mdw
+ * Change interface for suggested destinations.
+ *
  * Revision 1.1  1999/11/20 22:24:44  mdw
  * Add Diffie-Hellman support.
  *
@@ -42,7 +45,9 @@
 #include <string.h>
 
 #include "dh.h"
+#include "fibrand.h"
 #include "mp.h"
+#include "mprand.h"
 #include "pgen.h"
 #include "rabin.h"
 
@@ -70,14 +75,14 @@ mp *dh_prime(mp *s, size_t n,
 {
   pgen pq, pp;
   int rc_q, rc_p;
-  mpw bw;
-  mp b;
+  grand *gr = fibrand_create(0);
+  mp *b = MP_NEW;
+  size_t sz = mp_bits(s);
 
   /* --- Initialize prime generators --- */
 
   rc_q = pgen_create(&pq, s);
   rc_p = pgen_muladd(&pp, &pq, 2, 1);
-  mp_build(&b, &bw, &bw + 1);
 
   /* --- Now step along until something crops up --- */
 
@@ -106,12 +111,12 @@ mp *dh_prime(mp *s, size_t n,
      */
 
     for (i = 0; i < 5; i++) {
-      bw = ptab[i];
+      b = mprand(b, sz, gr, 1);
       if (rc_q == PGEN_MAYBE &&
-	  (rc_q = rabin_test(&rq, &b)) == PGEN_COMPOSITE)
+	  (rc_q = rabin_test(&rq, b)) == PGEN_COMPOSITE)
 	break;
       if (rc_p == PGEN_MAYBE &&
-	  (rc_p = rabin_test(&rp, &b)) == PGEN_COMPOSITE)
+	  (rc_p = rabin_test(&rp, b)) == PGEN_COMPOSITE)
 	break;
       if (proc && proc(DHEV_PASS, arg))
 	break;
@@ -146,6 +151,8 @@ mp *dh_prime(mp *s, size_t n,
     mp *p = MP_COPY(pp.m);
     pgen_destroy(&pq);
     pgen_destroy(&pp);
+    mp_drop(b);
+    gr->ops->destroy(gr);
     return (p);
   }
 
@@ -154,6 +161,8 @@ mp *dh_prime(mp *s, size_t n,
 fail:
   pgen_destroy(&pq);
   pgen_destroy(&pp);
+  mp_drop(b);
+  gr->ops->destroy(gr);
   return (0);
 }
 
diff --git a/dsa-gen.c b/dsa-gen.c
index cbddb70..b719a23 100644
--- a/dsa-gen.c
+++ b/dsa-gen.c
@@ -1,6 +1,6 @@
 /* -*-c-*-
  *
- * $Id: dsa-gen.c,v 1.2 1999/11/20 22:23:48 mdw Exp $
+ * $Id: dsa-gen.c,v 1.3 1999/12/10 23:18:38 mdw Exp $
  *
  * Generate DSA shared parameters
  *
@@ -30,6 +30,9 @@
 /*----- Revision history --------------------------------------------------* 
  *
  * $Log: dsa-gen.c,v $
+ * Revision 1.3  1999/12/10 23:18:38  mdw
+ * Change interface for suggested destinations.
+ *
  * Revision 1.2  1999/11/20 22:23:48  mdw
  * Allow event handler to abort the search process.
  *
@@ -44,10 +47,11 @@
 #include <stdlib.h>
 
 #include "dsa.h"
+#include "fibrand.h"
 #include "mp.h"
+#include "mprand.h"
 #include "pgen.h"
 #include "rabin.h"
-#include "rc4.h"
 #include "sha.h"
 
 /*----- Main code ---------------------------------------------------------*/
@@ -57,8 +61,8 @@
  * Arguments:	@dsa_param *dp@ = where to store parameters
  *		@unsigned l@ = bitlength of @p@ in bits
  *		@const void *k@ = pointer to key material
- *		@int (*proc)(int ev, mp *m, void *p)@ = event procedure
  *		@size_t sz@ = size of key material
+ *		@int (*proc)(int ev, mp *m, void *p)@ = event procedure
  *
  * Returns:	Zero if all went well, nonzero if key material was
  *		unsuitable (one of the @DSAEV@ codes).
@@ -83,11 +87,11 @@ int dsa_seed(dsa_param *dp, unsigned l, const void *k, size_t sz,
   mp *q, *p, *g;
   mp *s = mp_loadb(MP_NEW, k, sz);
   octet *sbuf = xmalloc(sz);
-  rc4_ctx rc4;
+  grand *gr;
   int fail = 0;
 
   l /= 8;
-  rc4_init(&rc4, k, sz);
+  gr = fibrand_create(LOAD32(k));
 
   /* --- Generate @q@ --- */
 
@@ -115,6 +119,9 @@ int dsa_seed(dsa_param *dp, unsigned l, const void *k, size_t sz,
 
   /* --- Quick primality check --- */
 
+  if (proc && (fail = proc(DSAEV_FINDQ, 0, arg)) != 0)
+    goto fail_0;
+
   {
     pgen pg;
     int rc = pgen_create(&pg, q);
@@ -138,12 +145,10 @@ int dsa_seed(dsa_param *dp, unsigned l, const void *k, size_t sz,
     rabin r;
     int i;
     mp *g = MP_NEW;
-    octet gbuf[SHA_HASHSZ];
 
     rabin_create(&r, q);
     for (i = 0; i < 18; i++) {
-      rc4_encrypt(&rc4, 0, gbuf, sizeof(gbuf));
-      g = mp_loadb(g, gbuf, sizeof(gbuf));
+      g = mprand(g, 8 * SHA_HASHSZ, gr, 1);
       if (rabin_test(&r, g) == PGEN_COMPOSITE)
 	break;
       if (proc && (fail = proc(DSAEV_PASSQ, q, arg)) != 0)
@@ -156,12 +161,14 @@ int dsa_seed(dsa_param *dp, unsigned l, const void *k, size_t sz,
     if (i < 18) {
       if (proc)
 	fail = proc(DSAEV_FAILQ, q, arg);
-      if (fail)
+      if (!fail)
 	fail = DSAEV_FAILQ;
       goto fail_0;
     }
     if (proc && (fail = proc(DSAEV_GOODQ, q, arg)) != 0)
       goto fail_0;
+    if (proc && (fail = proc(DSAEV_FINDP, 0, arg)) != 0)
+      goto fail_0;
   }
 
   /* --- Now generate @p@ --- *
@@ -241,8 +248,7 @@ int dsa_seed(dsa_param *dp, unsigned l, const void *k, size_t sz,
 	if (proc && (fail = proc(DSAEV_TRYP, p, arg)) != 0)
 	  break;
 	for (i = 0; i < 5; i++) {
-	  rc4_encrypt(&rc4, 0, pbuf, psz - b);
-	  g = mp_loadb(g, pbuf, psz - b);
+	  g = mprand(g, 8 * (psz - b), gr, 1);
 	  if (rabin_test(&r, g) == PGEN_COMPOSITE)
 	    break;
 	  if (proc && (fail = proc(DSAEV_PASSP, p, arg)) != 0)
@@ -263,6 +269,8 @@ int dsa_seed(dsa_param *dp, unsigned l, const void *k, size_t sz,
 
       if (proc)
 	fail = proc(DSAEV_GOODP, p, arg);
+      if (proc && !fail)
+	fail = proc(DSAEV_FINDG, 0, arg);
       break;
     }
 
@@ -284,11 +292,8 @@ int dsa_seed(dsa_param *dp, unsigned l, const void *k, size_t sz,
 
     /* --- Calculate %$p' = (p - 1)/q$% --- */
 
-    {
-      mp *p1 = mp_sub(MP_NEW, p, MP_ONE);
-      mp_div(&pp, 0, p1, q);
-      mp_drop(p1);
-    }
+    pp = mp_sub(MP_NEW, p, MP_ONE);
+    mp_div(&pp, 0, pp, q);
 
     /* --- Now find %$h$% where %$g = h^{p'} \bmod p \ne 1$% --- */
 
@@ -298,7 +303,7 @@ int dsa_seed(dsa_param *dp, unsigned l, const void *k, size_t sz,
       hw = ptab[i];
       if (proc && (fail = proc(DSAEV_TRYH, &h, arg)) != 0)
 	break;
-      g = mpmont_exp(&mm, &h, pp);
+      g = mpmont_exp(&mm, MP_NEW, &h, pp);
       if (MP_CMP(g, !=, MP_ONE))
 	break;
       mp_drop(g);
@@ -306,6 +311,7 @@ int dsa_seed(dsa_param *dp, unsigned l, const void *k, size_t sz,
 	break;
     }
     mp_drop(pp);
+    mpmont_destroy(&mm);
     if (fail)
       goto fail_1;
     if (i >= NPRIME) {
@@ -323,6 +329,7 @@ int dsa_seed(dsa_param *dp, unsigned l, const void *k, size_t sz,
   dp->q = q;
   mp_drop(s);
   free(sbuf);
+  gr->ops->destroy(gr);
   return (0);
 
   /* --- Tidy up after failure --- */
@@ -335,6 +342,7 @@ fail_0:
   mp_drop(q);
   mp_drop(s);
   free(sbuf);
+  gr->ops->destroy(gr);
   return (-1);
 }
 
@@ -389,6 +397,7 @@ static int verify(dstr *v)
   if (!rc) {
     mp_drop(dp.q); mp_drop(dp.p); mp_drop(dp.g);
   }
+  assert(mparena_count(MPARENA_GLOBAL) == 0);
   return (ok);
 }
 
diff --git a/dsa-sign.c b/dsa-sign.c
index 8e5c997..613a07e 100644
--- a/dsa-sign.c
+++ b/dsa-sign.c
@@ -1,6 +1,6 @@
 /* -*-c-*-
  *
- * $Id: dsa-sign.c,v 1.1 1999/11/19 19:28:00 mdw Exp $
+ * $Id: dsa-sign.c,v 1.2 1999/12/10 23:18:38 mdw Exp $
  *
  * DSA signing operation
  *
@@ -30,6 +30,9 @@
 /*----- Revision history --------------------------------------------------* 
  *
  * $Log: dsa-sign.c,v $
+ * Revision 1.2  1999/12/10 23:18:38  mdw
+ * Change interface for suggested destinations.
+ *
  * Revision 1.1  1999/11/19 19:28:00  mdw
  * Implementation of the Digital Signature Algorithm.
  *
@@ -39,6 +42,7 @@
 
 #include "dsa.h"
 #include "mp.h"
+#include "mpbarrett.h"
 #include "mpmont.h"
 
 /*----- Main code ---------------------------------------------------------*/
@@ -46,9 +50,9 @@
 /* --- @dsa_mksig@ --- *
  *
  * Arguments:	@const dsa_param *dp@ = pointer to DSA parameters
- *		@const mp *a@ = secret signing key
- *		@const mp *m@ = message to be signed
- *		@const mp *k@ = random data
+ *		@mp *a@ = secret signing key
+ *		@mp *m@ = message to be signed
+ *		@mp *k@ = random data
  *		@mp **rr, **ss@ = where to put output parameters
  *
  * Returns:	---
@@ -56,45 +60,38 @@
  * Use:		Computes a DSA signature of a message.
  */
 
-void dsa_mksig(const dsa_param *dp, const mp *a, const mp *m, const mp *k,
-	       mp **rr, mp **ss)
+void dsa_mksig(const dsa_param *dp, mp *a, mp *m, mp *k, mp **rr, mp **ss)
 {
-  mpmont pm, qm;
-  mp *k1, *r;
-  mp *rrr, *ar;
-
-  /* --- Create the Montgomery contexts --- */
-
-  mpmont_create(&pm, dp->p);
-  mpmont_create(&qm, dp->q);
+  mpmont pm;
+  mpbarrett qb;
+  mp *k1 = MP_NEW, *r;
+  mp *ar;
 
   /* --- Compute %$r = (g^k \bmod p) \bmod q$% --- */
 
-  r = mpmont_exp(&pm, dp->g, k);
+  mpmont_create(&pm, dp->p);
+  r = mpmont_exp(&pm, MP_NEW, dp->g, k);
+  mpmont_destroy(&pm);
   mp_div(0, &r, r, dp->q);
-  *rr = r;
 
   /* --- Compute %$k^{-1} \bmod q$% --- */
 
-  mp_gcd(0, 0, &k1, dp->q, (mp *)k);
+  mp_gcd(0, 0, &k1, dp->q, k);
 
   /* --- Now for %$k^{-1}(m + ar)$% --- */
 
-  rrr = mpmont_mul(&qm, MP_NEW, r, qm.r2);
-  ar = mpmont_mul(&qm, MP_NEW, a, rrr);
+  mpbarrett_create(&qb, dp->q);
+  ar = mp_mul(MP_NEW, a, r);
   ar = mp_add(ar, ar, m);
-  if (MP_CMP(ar, >=, dp->q))
-    ar = mp_sub(ar, ar, dp->q);
-  rrr = mpmont_mul(&qm, rrr, ar, qm.r2);
-  ar = mpmont_mul(&qm, ar, rrr, k1);
+  ar = mpbarrett_reduce(&qb, ar, ar);
+  ar = mp_mul(ar, ar, k1);
+  ar = mpbarrett_reduce(&qb, ar, ar);
+  mpbarrett_destroy(&qb);
+  MP_DROP(k1);
+  if (*rr) MP_DROP(*rr);
+  if (*ss) MP_DROP(*ss);
+  *rr = r;
   *ss = ar;
-
-  /* --- Tidy things up a little --- */
-
-  mp_drop(rrr);
-  mp_drop(k1);
-  mpmont_destroy(&pm);
-  mpmont_destroy(&qm);
 }
 
 /* --- @dsa_sign@ --- *
@@ -122,7 +119,7 @@ void dsa_sign(dsa_param *dp, mp *a,
 {
   mp *mm = mp_loadb(MP_NEW, m, msz);
   mp *km = mp_loadb(MP_NEW, k, ksz);
-  mp *rm, *sm;
+  mp *rm = MP_NEW, *sm = MP_NEW;
   dsa_mksig(dp, a, mm, km, &rm, &sm);
   mp_storeb(rm, r, rsz);
   mp_storeb(sm, s, ssz);
@@ -191,6 +188,7 @@ static int verify(dstr *v)
   mp_drop(dp.q);
   mp_drop(dp.g);
   mp_drop(x);
+  assert(mparena_count(MPARENA_GLOBAL) == 0);
   return (ok);
 }
 
diff --git a/dsa-verify.c b/dsa-verify.c
index 4b931bf..fee0dd3 100644
--- a/dsa-verify.c
+++ b/dsa-verify.c
@@ -1,6 +1,6 @@
 /* -*-c-*-
  *
- * $Id: dsa-verify.c,v 1.2 1999/11/23 00:20:04 mdw Exp $
+ * $Id: dsa-verify.c,v 1.3 1999/12/10 23:18:38 mdw Exp $
  *
  * DSA signature verification
  *
@@ -30,6 +30,9 @@
 /*----- Revision history --------------------------------------------------* 
  *
  * $Log: dsa-verify.c,v $
+ * Revision 1.3  1999/12/10 23:18:38  mdw
+ * Change interface for suggested destinations.
+ *
  * Revision 1.2  1999/11/23 00:20:04  mdw
  * Remove stray debugging code.
  *
@@ -49,17 +52,16 @@
 /* --- @dsa_vrfy@ --- *
  *
  * Arguments:	@const dsa_param *dp@ = pointer to DSA parameters
- *		@const mp *y@ = public verification key
- *		@const mp *m@ = message which was signed
- *		@const mp *r, *s@ = the signature
+ *		@mp *y@ = public verification key
+ *		@mp *m@ = message which was signed
+ *		@mp *r, *s@ = the signature
  *
  * Returns:	Zero if the signature is a forgery, nonzero if it's valid.
  *
  * Use:		Verifies a DSA digital signature.
  */
 
-int dsa_vrfy(const dsa_param *dp, const mp *y,
-	     const mp *m, const mp *r, const mp *s)
+int dsa_vrfy(const dsa_param *dp, mp *y, mp *m, mp *r, mp *s)
 {
   mpmont pm, qm;
   mp *w;
@@ -81,8 +83,8 @@ int dsa_vrfy(const dsa_param *dp, const mp *y,
   /* --- Compute %$w = s^{-1} \bmod q$% --- */
 
   {
-    mp *z;
-    mp_gcd(0, 0, &z, dp->q, (mp *)s);
+    mp *z = MP_NEW;
+    mp_gcd(0, 0, &z, dp->q, s);
     w = mpmont_mul(&qm, MP_NEW, z, qm.r2);
     mp_drop(z);
   }
@@ -97,8 +99,8 @@ int dsa_vrfy(const dsa_param *dp, const mp *y,
   /* --- Do the exponentiation and take residue mod @q@ --- */
 
   f[0].base = dp->g;
-  f[1].base = (mp *)y;
-  w = mpmont_mexp(&pm, f, 2);
+  f[1].base = y;
+  w = mpmont_mexp(&pm, MP_NEW, f, 2);
   mp_div(0, &w, w, dp->q);
   ok = MP_CMP(w, ==, r);
 
@@ -114,7 +116,7 @@ int dsa_vrfy(const dsa_param *dp, const mp *y,
 /* --- @dsa_verify@ --- *
  *
  * Arguments:	@const dsa_param *dp@ = pointer to DSA parameters
- *		@const mp *y@ = public verification key
+ *		@mp *y@ = public verification key
  *		@const void *m@ = pointer to message block
  *		@size_t msz@ = size of message block
  *		@const void *r@ = pointer to @r@ signature half
@@ -127,7 +129,7 @@ int dsa_vrfy(const dsa_param *dp, const mp *y,
  * Use:		Verifies a DSA digital signature.
  */
 
-int dsa_verify(const dsa_param *dp, const mp *y,
+int dsa_verify(const dsa_param *dp, mp *y,
 	       const void *m, size_t msz,
 	       const void *r, size_t rsz,
 	       const void *s, size_t ssz)
@@ -191,6 +193,7 @@ static int verify(int good, dstr *v)
   mp_drop(dp.q);
   mp_drop(dp.g);
   mp_drop(y);
+  assert(mparena_count(MPARENA_GLOBAL) == 0);
   return (ok);
 }
 
diff --git a/mp-arith.c b/mp-arith.c
index d8381ee..6f00656 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.2 1999/12/10 23:18:39 mdw Exp $
  *
  * Basic arithmetic on multiprecision integers
  *
@@ -30,6 +30,9 @@
 /*----- Revision history --------------------------------------------------* 
  *
  * $Log: mp-arith.c,v $
+ * 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 +42,10 @@
 
 #include "mp.h"
 
+/*----- Macros ------------------------------------------------------------*/
+
+#define MAX(x, y) ((x) >= (y) ? (x) : (y))
+
 /*----- Main code ---------------------------------------------------------*/
 
 /* --- @mp_2c@ --- *
@@ -84,13 +91,13 @@ 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);
   mpx_lsl(d->v, d->vl, a->v, a->vl, n);
@@ -102,13 +109,13 @@ 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));
   mpx_lsr(d->v, d->vl, a->v, a->vl, n);
@@ -138,19 +145,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_MODIFY(d, MAX(MP_LEN(a), MP_LEN(b)) + 1);
   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 +169,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_MODIFY(d, MAX(MP_LEN(a), MP_LEN(b)) + 1);
   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 +195,67 @@ 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;
+  a = MP_COPY(a);
+  b = MP_COPY(b);
+
   MP_MODIFY(d, MP_LEN(a) + MP_LEN(b));
-  mpx_umul(d->v, d->vl, a->v, a->vl, b->v, b->vl);
+  if (MP_LEN(a) <= KARATSUBA_CUTOFF || MP_LEN(b) <= KARATSUBA_CUTOFF)
+    mpx_umul(d->v, d->vl, a->v, a->vl, b->v, b->vl);
+  else {
+    size_t m = MAX(MP_LEN(a), MP_LEN(b)) * 2 + KARATSUBA_SLOP;
+    mpw *s;
+    m += 32;
+    s = MP_ALLOC(m);
+    mpx_kmul(d->v, d->vl, a->v, a->vl, b->v, b->vl, s, s + m);
+    MP_FREE(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_MODIFY(d, 2 * m);
+  if (m > KARATSUBA_CUTOFF) {
+    mpw *s;
+    m = 2 * (m + 1) + 32;
+    s = MP_ALLOC(m);
+    mpx_kmul(d->v, d->vl, a->v, a->vl, a->v, a->vl, s, s + m);
+    MP_FREE(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,7 +271,7 @@ 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;
@@ -267,12 +296,13 @@ void mp_div(mp **qq, mp **rr, const mp *a, const mp *b)
     if (MP_LEN(b) > rq)
       rq = MP_LEN(b);
 
+    b = MP_COPY(b);
     if (r == a) {
-      MP_SPLIT(r);
+      MP_SPLIT(a);
+      a = r = MP_COPY(a);
       MP_ENSURE(r, MP_LEN(r) + 2);
     } else {
-      if (r == b)
-	r = MP_NEW;
+      a = MP_COPY(a);
       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));
@@ -281,8 +311,6 @@ void mp_div(mp **qq, mp **rr, const mp *a, const mp *b)
 
   /* --- Fix up the quotient too --- */
 
-  if (q == a || q == b)
-    q = MP_NEW;
   MP_MODIFY(q, MP_LEN(a));
 
   /* --- Perform the calculation --- */
@@ -298,8 +326,8 @@ void mp_div(mp **qq, mp **rr, const mp *a, const mp *b)
 
   q->f = ((a->f | b->f) & MP_BURN) | ((a->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);
@@ -326,6 +354,8 @@ void mp_div(mp **qq, mp **rr, const mp *a, const mp *b)
     *rr = r;
   }
 
+  MP_DROP(a);
+  MP_DROP(b);
   MP_FREE(sv);
 }
 
@@ -348,7 +378,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 +389,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 +399,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 +407,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 +429,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);
 }
 
diff --git a/mp-gcd.c b/mp-gcd.c
index 19192f5..c772512 100644
--- a/mp-gcd.c
+++ b/mp-gcd.c
@@ -1,6 +1,6 @@
 /* -*-c-*-
  *
- * $Id: mp-gcd.c,v 1.2 1999/11/22 20:49:56 mdw Exp $
+ * $Id: mp-gcd.c,v 1.3 1999/12/10 23:18:39 mdw Exp $
  *
  * Extended GCD calculation
  *
@@ -30,6 +30,9 @@
 /*----- Revision history --------------------------------------------------* 
  *
  * $Log: mp-gcd.c,v $
+ * Revision 1.3  1999/12/10 23:18:39  mdw
+ * Change interface for suggested destinations.
+ *
  * Revision 1.2  1999/11/22 20:49:56  mdw
  * Fix bug which failed to favour `x' when `y' wasn't wanted and the two
  * arguments needed swapping.
@@ -54,9 +57,7 @@
  *
  * Use:		Calculates @gcd(a, b)@, and two numbers @x@ and @y@ such that
  *		@ax + by = gcd(a, b)@.  This is useful for computing modular
- *		inverses.  Neither @a@ nor @b@ may be zero.  Note that,
- *		unlike @mp_div@ for example, it is not possible to specify
- *		explicit destinations -- new MPs are always allocated.
+ *		inverses.
  */
 
 void mp_gcd(mp **gcd, mp **xx, mp **yy, mp *a, mp *b)
@@ -65,14 +66,73 @@ void mp_gcd(mp **gcd, mp **xx, mp **yy, mp *a, mp *b)
   mp *x = MP_ZERO, *y = MP_ONE;
   mp *u, *v;
   size_t shift = 0;
-  int ext = xx || yy;
-  int swap = 0;
+  unsigned f = 0;
+
+  enum {
+    f_swap = 1u,
+    f_aneg = 2u,
+    f_bneg = 4u,
+    f_ext = 8u
+  };
+
+  /* --- Sort out some initial flags --- */
 
-  /* --- Ensure that @a@ is larger than @b@ --- */
+  if (xx || yy)
+    f |= f_ext;
 
-  if (MP_CMP(a, <, b)) {
+  if (a->f & MP_NEG)
+    f |= f_aneg;
+  if (b->f & MP_NEG)
+    f |= f_bneg;
+
+  /* --- Ensure that @a@ is larger than @b@ --- *
+   *
+   * Use absolute values here!
+   */
+
+  if (MPX_UCMP(a->v, a->vl, <, b->v, b->vl)) {
     { mp *t = a; a = b; b = t; }
-    swap = 1;
+    f |= f_swap;
+  }
+
+  /* --- Check for zeroness --- */
+
+  if (MP_CMP(b, ==, MP_ZERO)) {
+
+    /* --- Store %$|a|$% as the GCD --- */
+
+    if (gcd) {
+      if (*gcd) MP_DROP(*gcd);
+      a = MP_COPY(a);
+      if (a->f & MP_NEG) {
+	MP_SPLIT(a);
+	a->f &= ~MP_NEG;
+	f |= f_aneg;
+      }
+      *gcd = a;
+    }
+
+    /* --- Store %$1$% and %$0$% in the appropriate bins --- */
+
+    if (f & f_ext) {
+      if (f & f_swap) {
+	mp **t = xx; xx = yy; yy = t;
+      }
+      if (xx) {
+	if (*xx) MP_DROP(*xx);
+	if (MP_CMP(a, ==, MP_ZERO))
+	  *xx = MP_ZERO;
+	else if (f & f_aneg)
+	  *xx = MP_MONE;
+	else
+	  *xx = MP_ONE;
+      }
+      if (yy) {
+	if (*yy) MP_DROP(*yy);
+	*yy = MP_ZERO;
+      }
+    }
+    return;
   }
 
   /* --- Take a reference to the arguments --- */
@@ -82,13 +142,14 @@ void mp_gcd(mp **gcd, mp **xx, mp **yy, mp *a, mp *b)
 
   /* --- Make sure @a@ and @b@ are not both even --- */
 
+  MP_SPLIT(a); a->f &= ~MP_NEG;
+  MP_SPLIT(b); b->f &= ~MP_NEG;
+
   if (((a->v[0] | b->v[0]) & 1) == 0) {
     mpscan asc, bsc;
 
     /* --- Break off my copies --- */
 
-    MP_SPLIT(a);
-    MP_SPLIT(b);
     MP_SCAN(&asc, a);
     MP_SCAN(&bsc, b);
 
@@ -130,7 +191,7 @@ void mp_gcd(mp **gcd, mp **xx, mp **yy, mp *a, mp *b)
 	MP_STEP(&xsc); MP_STEP(&ysc);
 	if (MP_BIT(&sc))
 	  break;
-	if (ext && (MP_BIT(&xsc) | MP_BIT(&ysc))) {
+	if ((f & f_ext) && (MP_BIT(&xsc) | MP_BIT(&ysc))) {
 	  if (n) {
 	    X = mp_lsr(X, X, n);
 	    Y = mp_lsr(Y, Y, n);
@@ -147,7 +208,7 @@ void mp_gcd(mp **gcd, mp **xx, mp **yy, mp *a, mp *b)
 
       if (nn) {
 	u = mp_lsr(u, u, nn);
-	if (ext && n) {
+	if ((f & f_ext) && n) {
 	  X = mp_lsr(X, X, n);
 	  Y = mp_lsr(Y, Y, n);
 	}
@@ -167,7 +228,7 @@ void mp_gcd(mp **gcd, mp **xx, mp **yy, mp *a, mp *b)
 	MP_STEP(&xsc); MP_STEP(&ysc);
 	if (MP_BIT(&sc))
 	  break;
-	if (ext && (MP_BIT(&xsc) | MP_BIT(&ysc))) {
+	if ((f & f_ext) && (MP_BIT(&xsc) | MP_BIT(&ysc))) {
 	  if (n) {
 	    x = mp_lsr(x, x, n);
 	    y = mp_lsr(y, y, n);
@@ -184,7 +245,7 @@ void mp_gcd(mp **gcd, mp **xx, mp **yy, mp *a, mp *b)
 
       if (nn) {
 	v = mp_lsr(v, v, nn);
-	if (ext && n) {
+	if ((f & f_ext) && n) {
 	  x = mp_lsr(x, x, n);
 	  y = mp_lsr(y, y, n);
 	}
@@ -195,13 +256,13 @@ void mp_gcd(mp **gcd, mp **xx, mp **yy, mp *a, mp *b)
 
     if (MP_CMP(u, >=, v)) {
       u = mp_sub(u, u, v);
-      if (ext) {
+      if (f & f_ext) {
 	X = mp_sub(X, X, x);
 	Y = mp_sub(Y, Y, y);
       }
     } else {
       v = mp_sub(v, v, u);
-      if (ext) {
+      if (f & f_ext) {
 	x = mp_sub(x, x, X);
 	y = mp_sub(y, y, Y);
       }
@@ -213,22 +274,34 @@ void mp_gcd(mp **gcd, mp **xx, mp **yy, mp *a, mp *b)
 
   /* --- Write the results out --- */
 
-  if (gcd)
-    *gcd = mp_lsl(v, v, shift);
-  else
+  if (!gcd)
     MP_DROP(v);
+  else {
+    if (*gcd) MP_DROP(*gcd);
+    *gcd = mp_lsl(v, v, shift);
+  }
 
   /* --- Perform a little normalization --- *
    *
    * Ensure that the coefficient returned is positive, if there is only one.
-   * If there are two, favour @y@.
+   * If there are two, favour @y@.  Of course, if the original arguments were
+   * negative then I'll need to twiddle their signs as well.
    */
 
-  if (ext) {
-    if (swap) {
+  if (f & f_ext) {
+
+    /* --- If @a@ and @b@ got swapped, swap the coefficients back --- */
+
+    if (f & f_swap) {
       mp *t = x; x = y; y = t;
       t = a; a = b; b = t;
     }
+
+    /* --- Sort out the signs --- *
+     *
+     * Note that %$ax + by = a(x - b) + b(y + a)$%.
+     */
+
     if (yy) {
       if (y->f & MP_NEG) {
 	y = mp_add(y, y, a);
@@ -237,8 +310,28 @@ void mp_gcd(mp **gcd, mp **xx, mp **yy, mp *a, mp *b)
     } else if (x->f & MP_NEG)
       x = mp_add(x, x, b);
 
-    if (xx) *xx = x; else MP_DROP(x);
-    if (yy) *yy = y; else MP_DROP(y);
+    /* --- Twiddle the signs --- */
+
+    if (f & f_aneg)
+      x->f ^= MP_NEG;
+    if (f & f_bneg)
+      y->f ^= MP_NEG;
+
+    /* --- Store the results --- */
+
+    if (!xx)
+      MP_DROP(x);
+    else {
+      if (*xx) MP_DROP(*xx);
+      *xx = x;
+    }
+
+    if (!yy)
+      MP_DROP(y);
+    else {
+      if (*yy) MP_DROP(*yy);
+      *yy = y;
+    }
   }
 
   MP_DROP(u);
@@ -259,7 +352,7 @@ static int gcd(dstr *v)
   mp *x = *(mp **)v[3].buf;
   mp *y = *(mp **)v[4].buf;
 
-  mp *gg, *xx, *yy;
+  mp *gg = MP_NEW, *xx = MP_NEW, *yy = MP_NEW;
   mp_gcd(&gg, &xx, &yy, a, b);
   if (MP_CMP(x, !=, xx)) {
     fputs("\n*** mp_gcd(x) failed", stderr);
@@ -301,6 +394,7 @@ static int gcd(dstr *v)
   }
   MP_DROP(a); MP_DROP(b); MP_DROP(g); MP_DROP(x); MP_DROP(y);
   MP_DROP(gg); MP_DROP(xx); MP_DROP(yy);
+  assert(mparena_count(MPARENA_GLOBAL) == 0);
   return (ok);
 }
 
diff --git a/mpmont-mexp.c b/mpmont-mexp.c
index 1d1266f..b541f8c 100644
--- a/mpmont-mexp.c
+++ b/mpmont-mexp.c
@@ -1,6 +1,6 @@
 /* -*-c-*-
  *
- * $Id: mpmont-mexp.c,v 1.2 1999/11/21 11:35:10 mdw Exp $
+ * $Id: mpmont-mexp.c,v 1.3 1999/12/10 23:18:39 mdw Exp $
  *
  * Multiplle simultaneous exponentiations
  *
@@ -30,6 +30,9 @@
 /*----- Revision history --------------------------------------------------* 
  *
  * $Log: mpmont-mexp.c,v $
+ * Revision 1.3  1999/12/10 23:18:39  mdw
+ * Change interface for suggested destinations.
+ *
  * Revision 1.2  1999/11/21 11:35:10  mdw
  * Performance improvement: use @mp_sqr@ and @mpmont_reduce@ instead of
  * @mpmont_mul@ for squaring in exponentiation.
@@ -49,6 +52,7 @@
 /* --- @mpmont_mexpr@ --- *
  *
  * Arguments:	@mpmont *mm@ = pointer to Montgomery reduction context
+ *		@mp *d@ = fake destination
  *		@mpmont_factor *f@ = pointer to array of factors
  *		@size_t n@ = number of factors supplied
  *
@@ -64,7 +68,7 @@ typedef struct scan {
   mpw w;
 } scan;
 
-mp *mpmont_mexpr(mpmont *mm, mpmont_factor *f, size_t n)
+mp *mpmont_mexpr(mpmont *mm, mp *d, mpmont_factor *f, size_t n)
 {
   size_t vn = 1 << n;
   mp **v = xmalloc(vn * sizeof(mp *));
@@ -186,12 +190,16 @@ mp *mpmont_mexpr(mpmont *mm, mpmont_factor *f, size_t n)
     free(s);
   }
 
+  if (d != MP_NEW)
+    MP_DROP(d);
+
   return (a);
 }
 
 /* --- @mpmont_mexp@ --- *
  *
  * Arguments:	@mpmont *mm@ = pointer to Montgomery reduction context
+ *		@mp *d@ = fake destination
  *		@mpmont_factor *f@ = pointer to array of factors
  *		@size_t n@ = number of factors supplied
  *
@@ -200,9 +208,9 @@ mp *mpmont_mexpr(mpmont *mm, mpmont_factor *f, size_t n)
  * Use:		Convenient interface over @mpmont_mexpr@.
  */
 
-mp *mpmont_mexp(mpmont *mm, mpmont_factor *f, size_t n)
+mp *mpmont_mexp(mpmont *mm, mp *d, mpmont_factor *f, size_t n)
 {
-  mp *d = mpmont_mexpr(mm, f, n);
+  d = mpmont_mexpr(mm, d, f, n);
   d = mpmont_reduce(mm, d, d);
   return (d);
 }
@@ -230,7 +238,7 @@ static int verify(size_t n, dstr *v)
 
   rr = *(mp **)v[j].buf;
   mpmont_create(&mm, m);
-  r = mpmont_mexp(&mm, f, n);
+  r = mpmont_mexp(&mm, MP_NEW, f, n);
   if (MP_CMP(r, !=, rr)) {
     fputs("\n*** mexp failed\n", stderr);
     fputs("m = ", stderr); mp_writefile(m, stderr, 10);
@@ -254,6 +262,7 @@ static int verify(size_t n, dstr *v)
   MP_DROP(r);
   MP_DROP(rr);
   mpmont_destroy(&mm);
+  assert(mparena_count(MPARENA_GLOBAL) == 0);
   return (ok);
 }
 
diff --git a/mpmont.c b/mpmont.c
index 7522dba..d061934 100644
--- a/mpmont.c
+++ b/mpmont.c
@@ -1,6 +1,6 @@
 /* -*-c-*-
  *
- * $Id: mpmont.c,v 1.5 1999/11/22 13:58:40 mdw Exp $
+ * $Id: mpmont.c,v 1.6 1999/12/10 23:18:39 mdw Exp $
  *
  * Montgomery reduction
  *
@@ -30,6 +30,9 @@
 /*----- Revision history --------------------------------------------------* 
  *
  * $Log: mpmont.c,v $
+ * Revision 1.6  1999/12/10 23:18:39  mdw
+ * Change interface for suggested destinations.
+ *
  * Revision 1.5  1999/11/22 13:58:40  mdw
  * Add an option to disable Montgomery reduction, so that performance
  * comparisons can be done.
@@ -77,6 +80,7 @@
  * Returns:	---
  *
  * Use:		Initializes a Montgomery reduction context ready for use.
+ *		The argument @m@ must be a positive odd integer.
  */
 
 #ifdef MPMONT_DISABLE
@@ -93,6 +97,12 @@ void mpmont_create(mpmont *mm, mp *m)
 
 void mpmont_create(mpmont *mm, mp *m)
 {
+  /* --- Validate the arguments --- */
+
+  assert(((void)"Montgomery modulus must be positive",
+	  (m->f & MP_NEG) == 0));
+  assert(((void)"Montgomery modulus must be odd", m->v[0] & 1));
+
   /* --- Take a copy of the modulus --- */
 
   mp_shrink(m);
@@ -106,7 +116,7 @@ void mpmont_create(mpmont *mm, mp *m)
   {
     mpw av[2] = { 0, 1 };
     mp a, b;
-    mp *i;
+    mp *i = MP_NEW;
     mpw mi;
 
     mp_build(&a, av, av + 2);
@@ -158,14 +168,14 @@ void mpmont_destroy(mpmont *mm)
  *
  * Arguments:	@mpmont *mm@ = pointer to Montgomery reduction context
  *		@mp *d@ = destination
- *		@const mp *a@ = source, assumed positive
+ *		@mp *a@ = source, assumed positive
  *
  * Returns:	Result, %$a R^{-1} \bmod m$%.
  */
 
 #ifdef MPMONT_DISABLE
 
-mp *mpmont_reduce(mpmont *mm, mp *d, const mp *a)
+mp *mpmont_reduce(mpmont *mm, mp *d, mp *a)
 {
   mp_div(0, &d, a, mm->m);
   return (d);
@@ -173,10 +183,10 @@ mp *mpmont_reduce(mpmont *mm, mp *d, const mp *a)
 
 #else
 
-mp *mpmont_reduce(mpmont *mm, mp *d, const mp *a)
+mp *mpmont_reduce(mpmont *mm, mp *d, mp *a)
 {
   mpw *dv, *dvl;
-  const mpw *mv, *mvl;
+  mpw *mv, *mvl;
   size_t n;
 
   /* --- Initial conditioning of the arguments --- */
@@ -219,14 +229,14 @@ mp *mpmont_reduce(mpmont *mm, mp *d, const mp *a)
  *
  * Arguments:	@mpmont *mm@ = pointer to Montgomery reduction context
  *		@mp *d@ = destination
- *		@const mp *a, *b@ = sources, assumed positive
+ *		@mp *a, *b@ = sources, assumed positive
  *
  * Returns:	Result, %$a b R^{-1} \bmod m$%.
  */
 
 #ifdef MPMONT_DISABLE
 
-mp *mpmont_mul(mpmont *mm, mp *d, const mp *a, const mp *b)
+mp *mpmont_mul(mpmont *mm, mp *d, mp *a, mp *b)
 {
   d = mp_mul(d, a, b);
   mp_div(0, &d, d, mm->m);
@@ -235,59 +245,69 @@ mp *mpmont_mul(mpmont *mm, mp *d, const mp *a, const mp *b)
 
 #else
 
-mp *mpmont_mul(mpmont *mm, mp *d, const mp *a, const mp *b)
+mp *mpmont_mul(mpmont *mm, mp *d, mp *a, mp *b)
 {
-  mpw *dv, *dvl;
-  const mpw *av, *avl;
-  const mpw *bv, *bvl;
-  const mpw *mv, *mvl;
-  mpw y;
-  size_t n, i;
-
-  /* --- Initial conditioning of the arguments --- */
+  if (MP_LEN(a) > KARATSUBA_CUTOFF && MP_LEN(b) > KARATSUBA_CUTOFF) {
+    d = mp_mul(d, a, b);
+    d = mpmont_reduce(mm, d, d);
+  } else {
+    mpw *dv, *dvl;
+    mpw *av, *avl;
+    mpw *bv, *bvl;
+    mpw *mv, *mvl;
+    mpw y;
+    size_t n, i;
+
+    /* --- Initial conditioning of the arguments --- */
+
+    if (MP_LEN(a) > MP_LEN(b)) {
+      mp *t = a; a = b; b = t;
+    }
+    n = MP_LEN(mm->m);
 
-  if (MP_LEN(a) > MP_LEN(b)) {
-    const mp *t = a; a = b; b = t;
-  }
-  n = MP_LEN(mm->m);
-    
-  MP_MODIFY(d, 2 * n + 1);
-  dv = d->v; dvl = d->vl;
-  MPX_ZERO(dv, dvl);
-  av = a->v; avl = a->vl;
-  bv = b->v; bvl = b->vl;
-  mv = mm->m->v; mvl = mm->m->vl;
-  y = *bv;
+    a = MP_COPY(a);
+    b = MP_COPY(b);
+    MP_MODIFY(d, 2 * n + 1);
+    dv = d->v; dvl = d->vl;
+    MPX_ZERO(dv, dvl);
+    av = a->v; avl = a->vl;
+    bv = b->v; bvl = b->vl;
+    mv = mm->m->v; mvl = mm->m->vl;
+    y = *bv;
+
+    /* --- Montgomery multiplication phase --- */
+
+    i = 0;
+    while (i < n && av < avl) {
+      mpw x = *av++;
+      mpw u = MPW((*dv + x * y) * mm->mi);
+      MPX_UMLAN(dv, dvl, bv, bvl, x);
+      MPX_UMLAN(dv, dvl, mv, mvl, u);
+      dv++;
+      i++;
+    }
 
-  /* --- Montgomery multiplication phase --- */
+    /* --- Simpler Montgomery reduction phase --- */
 
-  i = 0;
-  while (i < n && av < avl) {
-    mpw x = *av++;
-    mpw u = MPW((*dv + x * y) * mm->mi);
-    MPX_UMLAN(dv, dvl, bv, bvl, x);
-    MPX_UMLAN(dv, dvl, mv, mvl, u);
-    dv++;
-    i++;
-  }
+    while (i < n) {
+      mpw u = MPW(*dv * mm->mi);
+      MPX_UMLAN(dv, dvl, mv, mvl, u);
+      dv++;
+      i++;
+    }
 
-  /* --- Simpler Montgomery reduction phase --- */
+    /* --- Done --- */
 
-  while (i < n) {
-    mpw u = MPW(*dv * mm->mi);
-    MPX_UMLAN(dv, dvl, mv, mvl, u);
-    dv++;
-    i++;
+    memmove(d->v, dv, MPWS(dvl - dv));
+    d->vl -= dv - d->v;
+    MP_SHRINK(d);
+    d->f = (a->f | b->f) & MP_BURN;
+    if (MP_CMP(d, >=, mm->m))
+      d = mp_sub(d, d, mm->m);
+    MP_DROP(a);
+    MP_DROP(b);
   }
 
-  /* --- Done --- */
-
-  memmove(d->v, dv, MPWS(dvl - dv));
-  d->vl -= dv - d->v;
-  MP_SHRINK(d);
-  d->f = (a->f | b->f) & MP_BURN;
-  if (MP_CMP(d, >=, mm->m))
-    d = mp_sub(d, d, mm->m);
   return (d);
 }
 
@@ -296,17 +316,18 @@ mp *mpmont_mul(mpmont *mm, mp *d, const mp *a, const mp *b)
 /* --- @mpmont_expr@ --- *
  *
  * Arguments:	@mpmont *mm@ = pointer to Montgomery reduction context
- *		@const mp *a@ = base
- *		@const mp *e@ = exponent
+ *		@mp *d@ = fake destination
+ *		@mp *a@ = base
+ *		@mp *e@ = exponent
  *
  * Returns:	Result, %$a^e R \bmod m$%.
  */
 
-mp *mpmont_expr(mpmont *mm, const mp *a, const mp *e)
+mp *mpmont_expr(mpmont *mm, mp *d, mp *a, mp *e)
 {
   mpscan sc;
   mp *ar = mpmont_mul(mm, MP_NEW, a, mm->r2);
-  mp *d = MP_COPY(mm->r);
+  mp *x = MP_COPY(mm->r);
   mp *spare = MP_NEW;
 
   mp_scan(&sc, e);
@@ -322,8 +343,8 @@ mp *mpmont_expr(mpmont *mm, const mp *a, const mp *e)
 	  spare = ar; ar = dd;
 	  sq--;
 	}
-	dd = mpmont_mul(mm, spare, d, ar);
-	spare = d; d = dd;
+	dd = mpmont_mul(mm, spare, x, ar);
+	spare = x; x = dd;
       }
       sq++;
       if (!MP_STEP(&sc))
@@ -333,21 +354,24 @@ mp *mpmont_expr(mpmont *mm, const mp *a, const mp *e)
   MP_DROP(ar);
   if (spare != MP_NEW)
     MP_DROP(spare);
-  return (d);
+  if (d != MP_NEW)
+    MP_DROP(d);
+  return (x);
 }
 
 /* --- @mpmont_exp@ --- *
  *
  * Arguments:	@mpmont *mm@ = pointer to Montgomery reduction context
- *		@const mp *a@ = base
- *		@const mp *e@ = exponent
+ *		@mp *d@ = fake destination
+ *		@mp *a@ = base
+ *		@mp *e@ = exponent
  *
  * Returns:	Result, %$a^e \bmod m$%.
  */
 
-mp *mpmont_exp(mpmont *mm, const mp *a, const mp *e)
+mp *mpmont_exp(mpmont *mm, mp *d, mp *a, mp *e)
 {
-  mp *d = mpmont_expr(mm, a, e);
+  d = mpmont_expr(mm, d, a, e);
   d = mpmont_reduce(mm, d, d);
   return (d);
 }
@@ -399,6 +423,7 @@ static int tcreate(dstr *v)
   MP_DROP(r);
   MP_DROP(r2);
   mpmont_destroy(&mm);
+  assert(mparena_count(MPARENA_GLOBAL) == 0);
   return (ok);
 }
 
@@ -456,6 +481,7 @@ static int tmul(dstr *v)
   MP_DROP(b);
   MP_DROP(r);
   mpmont_destroy(&mm);
+  assert(mparena_count(MPARENA_GLOBAL) == 0);
   return ok;
 }
 
@@ -471,7 +497,7 @@ static int texp(dstr *v)
   mpmont mm;
   mpmont_create(&mm, m);
 
-  mr = mpmont_exp(&mm, a, b);
+  mr = mpmont_exp(&mm, MP_NEW, a, b);
 
   if (MP_CMP(mr, !=, r)) {
     fputs("\n*** montgomery modexp failed", stderr);
@@ -490,14 +516,15 @@ static int texp(dstr *v)
   MP_DROP(r);
   MP_DROP(mr);
   mpmont_destroy(&mm);
+  assert(mparena_count(MPARENA_GLOBAL) == 0);
   return ok;
 }
 
 
 static test_chunk tests[] = {
-  { "create", tcreate, { &type_mp, &type_mp, &type_mp, &type_mp } },
-  { "mul", tmul, { &type_mp, &type_mp, &type_mp, &type_mp } },
-  { "exp", texp, { &type_mp, &type_mp, &type_mp, &type_mp } },
+  { "create", tcreate, { &type_mp, &type_mp, &type_mp, &type_mp, 0 } },
+  { "mul", tmul, { &type_mp, &type_mp, &type_mp, &type_mp, 0 } },
+  { "exp", texp, { &type_mp, &type_mp, &type_mp, &type_mp, 0 } },
   { 0, 0, { 0 } },
 };