From: mdw <mdw>
Date: Sun, 21 Mar 2004 22:52:06 +0000 (+0000)
Subject: Merge and close elliptic curve branch.
X-Git-Url: https://git.distorted.org.uk/u/mdw/catacomb/commitdiff_plain/c3caa2face1cda7002eb58245ad75865bf437455

Merge and close elliptic curve branch.
---

c3caa2face1cda7002eb58245ad75865bf437455
diff --cc BRANCHES
index 2f32e82,2f32e82..5310815
--- a/BRANCHES
+++ b/BRANCHES
@@@ -5,4 -5,4 +5,4 @@@ For a branch FOO, we have
  	FOO-merge-N	Nth branch merge point
  
  `ec' -- elliptic curve work
--	No merges
++	ec-merge-1	Closed.
diff --cc Makefile.m4
index 814fd34,4a95044..9be48d5
--- a/Makefile.m4
+++ b/Makefile.m4
@@@ -1,6 -1,6 +1,6 @@@
 -## -*-makefile-*-
 +## -*-m4-*-
  ##
- ## $Id: Makefile.m4,v 1.66 2004/03/21 22:43:50 mdw Exp $
 -## $Id: Makefile.m4,v 1.60.2.2 2004/03/21 22:39:46 mdw Exp $
++## $Id: Makefile.m4,v 1.67 2004/03/21 22:52:06 mdw Exp $
  ##
  ## Makefile for Catacomb
  ##
@@@ -29,23 -29,11 +29,32 @@@
  ##----- Revision history ----------------------------------------------------
  ##
  ## $Log: Makefile.m4,v $
 -## Revision 1.60.2.2  2004/03/21 22:39:46  mdw
 -## Elliptic curves on binary fields work.
++## Revision 1.67  2004/03/21 22:52:06  mdw
++## Merge and close elliptic curve branch.
+ ##
 -## Revision 1.60.2.1  2003/06/10 13:43:53  mdw
 -## Simple (non-projective) curves over prime fields now seem to work.
++##   Revision 1.60.2.2  2004/03/21 22:39:46  mdw
++##   Elliptic curves on binary fields work.
++##
++##   Revision 1.60.2.1  2003/06/10 13:43:53  mdw
++##   Simple (non-projective) curves over prime fields now seem to work.
++##
 +## Revision 1.66  2004/03/21 22:43:50  mdw
 +## New hash variant SHA224.
 +##
 +## Revision 1.65  2003/11/29 23:39:36  mdw
 +## Debianization.
 +##
 +## Revision 1.64  2003/11/10 22:18:30  mdw
 +## Build fixes.
 +##
 +## Revision 1.63  2003/10/17 16:30:46  mdw
 +## Report errors if key files don't exist!
 +##
 +## Revision 1.62  2003/10/12 15:02:09  mdw
 +## Reliability fixes.
 +##
 +## Revision 1.61  2003/10/11 21:02:33  mdw
 +## Import buf stuff from tripe.
  ##
  ## Revision 1.60  2003/05/16 01:12:37  mdw
  ## Ship `rc2-tab.h' and `skipjack-tab.h'.
@@@ -353,11 -341,14 +364,14 @@@ define(`MP_SOURCES'
  	exp.c mpcrt.c mpmul.c mprand.c \
  	mpbarrett.c mpbarrett-mexp.c mpbarrett-exp.h \
  	mpmont.c mpmont-mexp.c mpmont-exp.h \
 -	rho.c \
 +	rho.c buf.c \
- 	GF_SOURCES PGEN_SOURCES')
+ 	GF_SOURCES PGEN_SOURCES EC_SOURCES')
  
  define(`GF_SOURCES',
- 	`gfx.c gfx-kmul.c gfx-sqr.c')
+ 	`gfx.c gfx-kmul.c gfx-sqr.c gf-arith.c gf-gcd.c gfreduce.c')
+ 
+ define(`EC_SOURCES',
+ 	`field.c f-prime.c f-binpoly.c ec.c ec-prime.c ec-bin.c')
  
  define(`PGEN_SOURCES',
  	`pfilt.c rabin.c \
diff --cc calc/ec2.cal
index 0000000,0d7ceb9..3e89034
mode 000000,100644..100644
--- a/calc/ec2.cal
+++ b/calc/ec2.cal
@@@ -1,0 -1,180 +1,183 @@@
+ /* -*-apcalc-*-
+  *
 - * $Id: ec2.cal,v 1.1.2.1 2004/03/21 22:39:46 mdw Exp $
++ * $Id: ec2.cal,v 1.2 2004/03/21 22:52:06 mdw Exp $
+  *
+  * Testbed for elliptic curve arithmetic over binary fields
+  *
+  * (c) 2004 Straylight/Edgeware
+  */
+ 
+ /*----- Licensing notice --------------------------------------------------* 
+  *
+  * This file is part of Catacomb.
+  *
+  * Catacomb is free software; you can redistribute it and/or modify
+  * 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: ec2.cal,v $
++ * Revision 1.2  2004/03/21 22:52:06  mdw
++ * Merge and close elliptic curve branch.
++ *
+  * Revision 1.1.2.1  2004/03/21 22:39:46  mdw
+  * Elliptic curves on binary fields work.
+  *
+  * Revision 1.1.4.2  2004/03/20 00:13:31  mdw
+  * Projective coordinates for prime curves
+  *
+  * Revision 1.1.4.1  2003/06/10 13:43:53  mdw
+  * Simple (non-projective) curves over prime fields now seem to work.
+  *
+  * Revision 1.1  2000/10/08 16:01:37  mdw
+  * Prototypes of various bits of code.
+  *
+  */
+ 
+ /*----- Object types ------------------------------------------------------*/
+ 
+ obj ec2_curve { a, b, p };
+ obj ec2_pt { x, y, e };
+ obj ecpp_pt { x, y, z, e };
+ 
+ /*----- Main code ---------------------------------------------------------*/
+ 
+ define ec2_curve(a, b, p)
+ {
+   local obj ec2_curve e;
+   e.a = a;
+   e.b = b;
+   e.p = p;
+   return (e);
+ }
+ 
+ define ec2_pt(x, y, e)
+ {
+   local obj ec2_pt p;
+   p.x = x % e.p;
+   p.y = y % e.p;
+   p.e = e;
+   return (p);
+ }
+ 
+ define ec2_pt_print(a)
+ {
+   print "(" : a.x : ", " : a.y : ")" :;
+ }
+ 
+ define ec2_pt_add(a, b)
+ {
+   local e, alpha;
+   local obj ec2_pt d;
+ 
+   print "> ecadd: ", a, b;
+   if (a == 0)
+     d = b;
+   else if (b == 0)
+     d = a;
+   else if (!istype(a, b))
+     quit "bad type arguments to ec2_pt_add";
+   else if (a.e != b.e)
+     quit "points from different curves in ec2_pt_add";
+   else {
+     e = a.e;
+     if (a.x != b.x) {
+       alpha = ((a.y + b.y) * gf_inv(a.x + b.x, e.p)) % e.p;
+       d.x = (e.a + alpha^2 + alpha + a.x + b.x) % e.p;
+     } else if (a.y != b.y || a.x == gf(0))
+       return 0;
+     else {
+       alpha = a.x + a.y * gf_inv(a.x, e.p) % e.p;
+       d.x = (e.a + alpha^2 + alpha) % e.p;
+     }
+     d.y = ((a.x + d.x) * alpha + d.x + a.y) % e.p;
+     d.e = e;
+   }
+ 
+   print "< ecadd: ", d;
+   return (d);
+ }
+ 
+ define ec2_pt_dbl(a)
+ {
+   local e, alpha;
+   local obj ec2_pt d;
+   print "> ecdbl: ", a;
+   if (istype(a, 1))
+     return (0);
+   e = a.e;
+   alpha = a.x + a.y * gf_inv(a.x, e.p) % e.p;
+   d.x = (e.a + alpha^2 + alpha) % e.p;
+   d.y = ((a.x + d.x) * alpha + d.x + a.y) % e.p;
+   d.e = e;
+   print "< ecdbl: ", d;
+   return (d);
+ }
+ 
+ define ec2_pt_sub(a, b) { return ec2_pt_add(a, ec2_pt_neg(b)); }
+ 
+ define ec2_pt_neg(a)
+ {
+   local obj ec2_pt d;
+   d.x = a.x;
+   d.y = a.x + a.y;
+   d.e = a.e;
+   return (d);
+ }
+ 
+ define ec2_pt_check(a)
+ {
+   local e;
+ 
+   e = a.e;
+   if ((a.y^2 + a.x * a.y) % e.p != (a.x^3 + e.a * a.x^2 + e.b) % e.p)
+     quit "bad curve point";
+ }
+ 
+ define ec2_pt_mul(a, b)
+ {
+   local p, n;
+   local d;
+ 
+   if (istype(a, 1)) {
+     n = a;
+     p = b;
+   } else if (istype(b, 1)) {
+     n = b;
+     p = a;
+   } else
+     return (newerror("bad arguments to ec2_pt_mul"));
+ 
+   d = 0;
+   while (n) {
+     if (n & 1)
+       d += p;
+     n >>= 1;
+     p = ec2_pt_dbl(p);
+   }
+   return (d);
+ }
+ 
+ /*----- FIPS186-2 standard curves -----------------------------------------*/
+ 
+ b163 = ec2_curve(gf(1),gf(0x20a601907b8c953ca1481eb10512f78744a3205fd),
+ 		 gf(0x800000000000000000000000000000000000000c9));
+ b163_r = 5846006549323611672814742442876390689256843201587;
+ b163_g = ec2_pt(0x3f0eba16286a2d57ea0991168d4994637e8343e36,
+ 		0x0d51fbc6c71a0094fa2cdd545b11c5c0c797324f1, b163);
+ 
+ /*----- That's all, folks -------------------------------------------------*/
+ 
diff --cc calc/ecp.cal
index 04971aa,43ac1b3..7c560c5
--- a/calc/ecp.cal
+++ b/calc/ecp.cal
@@@ -1,6 -1,6 +1,6 @@@
  /* -*-apcalc-*-
   *
-  * $Id: ecp.cal,v 1.1 2000/10/08 16:01:37 mdw Exp $
 - * $Id: ecp.cal,v 1.1.4.2 2004/03/20 00:13:31 mdw Exp $
++ * $Id: ecp.cal,v 1.2 2004/03/21 22:52:06 mdw Exp $
   *
   * Testbed for elliptic curve arithmetic over prime fields
   *
@@@ -30,6 -30,12 +30,15 @@@
  /*----- Revision history --------------------------------------------------* 
   *
   * $Log: ecp.cal,v $
++ * Revision 1.2  2004/03/21 22:52:06  mdw
++ * Merge and close elliptic curve branch.
++ *
+  * Revision 1.1.4.2  2004/03/20 00:13:31  mdw
+  * Projective coordinates for prime curves
+  *
+  * Revision 1.1.4.1  2003/06/10 13:43:53  mdw
+  * Simple (non-projective) curves over prime fields now seem to work.
+  *
   * Revision 1.1  2000/10/08 16:01:37  mdw
   * Prototypes of various bits of code.
   *
diff --cc calc/gfx.cal
index 8d8fd00,5b19cb3..446061e
--- a/calc/gfx.cal
+++ b/calc/gfx.cal
@@@ -1,6 -1,6 +1,6 @@@
  /* -*-apcalc-*-
   *
-  * $Id: gfx.cal,v 1.1 2000/10/08 16:01:37 mdw Exp $
 - * $Id: gfx.cal,v 1.1.4.1 2004/03/21 22:39:46 mdw Exp $
++ * $Id: gfx.cal,v 1.2 2004/03/21 22:52:06 mdw Exp $
   *
   * Testbed for %$\gf{2}$% poltnomial arithmetic
   *
@@@ -30,6 -30,9 +30,12 @@@
  /*----- Revision history --------------------------------------------------* 
   *
   * $Log: gfx.cal,v $
++ * Revision 1.2  2004/03/21 22:52:06  mdw
++ * Merge and close elliptic curve branch.
++ *
+  * Revision 1.1.4.1  2004/03/21 22:39:46  mdw
+  * Elliptic curves on binary fields work.
+  *
   * Revision 1.1  2000/10/08 16:01:37  mdw
   * Prototypes of various bits of code.
   *
diff --cc configure.in
index 555e0ef,599cc28..7e2245d
--- a/configure.in
+++ b/configure.in
@@@ -1,6 -1,6 +1,6 @@@
 -dnl -*-fundamental-*-
 +dnl -*-m4-*-
  dnl
- dnl $Id: configure.in,v 1.26 2003/11/29 23:39:36 mdw Exp $
 -dnl $Id: configure.in,v 1.24.2.1 2003/06/10 13:43:53 mdw Exp $
++dnl $Id: configure.in,v 1.27 2004/03/21 22:52:06 mdw Exp $
  dnl
  dnl Autoconfiguration for Catacomb
  dnl
@@@ -29,11 -29,8 +29,17 @@@ dnl MA 02111-1307, USA
  dnl ----- Revision history --------------------------------------------------
  dnl
  dnl $Log: configure.in,v $
 -dnl Revision 1.24.2.1  2003/06/10 13:43:53  mdw
 -dnl Simple (non-projective) curves over prime fields now seem to work.
++dnl Revision 1.27  2004/03/21 22:52:06  mdw
++dnl Merge and close elliptic curve branch.
++dnl
++dnl   Revision 1.24.2.1  2003/06/10 13:43:53  mdw
++dnl   Simple (non-projective) curves over prime fields now seem to work.
++dnl
 +dnl Revision 1.26  2003/11/29 23:39:36  mdw
 +dnl Debianization.
 +dnl
 +dnl Revision 1.25  2003/10/11 21:02:33  mdw
 +dnl Import buf stuff from tripe.
  dnl
  dnl Revision 1.24  2003/05/16 00:30:28  mdw
  dnl Version bump.
@@@ -84,7 -81,7 +90,7 @@@ dn
  dnl --- Boring boilerplate ---
  
  AC_INIT(blkc.h)
- mdw_INIT_LIB(catacomb, Catacomb, 2.0.1)
 -mdw_INIT_LIB(catacomb, Catacomb, 2.1.0ec1)
++mdw_INIT_LIB(catacomb, Catacomb, 2.1.0)
  AM_CONFIG_HEADER(config.h)
  
  dnl --- Make sure I can compile and build libraries ---
diff --cc debian/changelog
index 5717a50,0000000..ce65b69
mode 100644,000000..100644
--- a/debian/changelog
+++ b/debian/changelog
@@@ -1,6 -1,0 +1,15 @@@
++catacomb (2.1.0) experimental; urgency=low
++
++  * Added support for elliptic curves, on both prime and binary fields
++    (polynomial basis only).  No actual crypto, but there's enough already
++    to do ECDH and stuff on well-known curves  Testing is currently a bit
++    patchy.
++
++ -- Mark Wooding <mdw@nsict.org>  Sun, 21 Mar 2004 22:47:56 +0000
++
 +catacomb (2.0.1) experimental; urgency=low
 +
 +  * Debianization!
 +  * (pixie): Don't report uninteresting errors when accepting connections.
 +
 + -- Mark Wooding <mdw@nsict.org>  Thu, 11 Dec 2003 10:47:59 +0000
diff --cc ec-bin.c
index 0000000,4f79c3d..3e85e65
mode 000000,100644..100644
--- a/ec-bin.c
+++ b/ec-bin.c
@@@ -1,0 -1,428 +1,431 @@@
+ /* -*-c-*-
+  *
 - * $Id: ec-bin.c,v 1.1.2.1 2004/03/21 22:39:46 mdw Exp $
++ * $Id: ec-bin.c,v 1.2 2004/03/21 22:52:06 mdw Exp $
+  *
+  * Arithmetic for elliptic curves over binary fields
+  *
+  * (c) 2004 Straylight/Edgeware
+  */
+ 
+ /*----- Licensing notice --------------------------------------------------* 
+  *
+  * This file is part of Catacomb.
+  *
+  * Catacomb is free software; you can redistribute it and/or modify
+  * 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: ec-bin.c,v $
++ * Revision 1.2  2004/03/21 22:52:06  mdw
++ * Merge and close elliptic curve branch.
++ *
+  * Revision 1.1.2.1  2004/03/21 22:39:46  mdw
+  * Elliptic curves on binary fields work.
+  *
+  */
+ 
+ /*----- Header files ------------------------------------------------------*/
+ 
+ #include <mLib/sub.h>
+ 
+ #include "ec.h"
+ 
+ /*----- Data structures ---------------------------------------------------*/
+ 
+ typedef struct ecctx {
+   ec_curve c;
+   mp *a, *b;
+   mp *bb;
+ } ecctx;
+ 
+ /*----- Main code ---------------------------------------------------------*/
+ 
+ static const ec_ops ec_binops, ec_binprojops;
+ 
+ static ec *ecneg(ec_curve *c, ec *d, const ec *p)
+ {
+   EC_COPY(d, p);
+   if (d->x)
+     d->y = F_ADD(c->f, d->y, d->y, d->x);
+   return (d);
+ }
+ 
+ static ec *ecprojneg(ec_curve *c, ec *d, const ec *p)
+ {
+   EC_COPY(d, p);
+   if (d->x) {
+     mp *t = F_MUL(c->f, MP_NEW, d->x, d->z);
+     d->y = F_ADD(c->f, d->y, d->y, t);
+     MP_DROP(t);
+   }
+   return (d);
+ }
+ 
+ static ec *ecfind(ec_curve *c, ec *d, mp *x)
+ {
+   /* write me */
+   return (0);
+ }
+ 
+ static ec *ecdbl(ec_curve *c, ec *d, const ec *a)
+ {
+   if (EC_ATINF(a) || F_ZEROP(c->f, a->x))
+     EC_SETINF(d);
+   else {
+     field *f = c->f;
+     ecctx *cc = (ecctx *)c;
+     mp *lambda;
+     mp *dx, *dy;
+ 
+     dx = F_INV(f, MP_NEW, a->x);	/* %$x^{-1}$% */
+     dy = F_MUL(f, MP_NEW, dx, a->y);	/* %$y/x$% */
+     lambda = F_ADD(f, dy, dy, a->x);	/* %$\lambda = x + y/x$% */
+ 
+     dx = F_SQR(f, dx, lambda);		/* %$\lambda^2$% */
+     dx = F_ADD(f, dx, dx, lambda);	/* %$\lambda^2 + \lambda$% */
+     dx = F_ADD(f, dx, dx, cc->a);     /* %$x' = a + \lambda^2 + \lambda$% */
+ 
+     dy = F_ADD(f, MP_NEW, a->x, dx);	/* %$ x + x' $% */
+     dy = F_MUL(f, dy, dy, lambda);	/* %$ (x + x') \lambda$% */
+     dy = F_ADD(f, dy, dy, a->y);	/* %$ (x + x') \lambda + y$% */
+     dy = F_ADD(f, dy, dy, dx);	    /* %$ y' = (x + x') \lambda + y + x'$% */
+ 
+     EC_DESTROY(d);
+     d->x = dx;
+     d->y = dy;
+     d->z = 0;
+     MP_DROP(lambda);
+   }
+   return (d);
+ }
+ 
+ static ec *ecprojdbl(ec_curve *c, ec *d, const ec *a)
+ {
+   if (EC_ATINF(a) || F_ZEROP(c->f, a->x))
+     EC_SETINF(d);
+   else {
+     field *f = c->f;
+     ecctx *cc = (ecctx *)c;
+     mp *dx, *dy, *dz, *u, *v;
+ 
+     dy = F_SQR(f, MP_NEW, a->z);	/* %$z^2$% */
+     dx = F_MUL(f, MP_NEW, dy, cc->bb);	/* %$c z^2$% */
+     dx = F_ADD(f, dx, dx, a->x);	/* %$x + c z^2$% */
+     dz = F_SQR(f, MP_NEW, dx);		/* %$(x + c z^2)^2$% */
+     dx = F_SQR(f, dx, dz);		/* %$x' = (x + c z^2)^4$% */
+ 
+     dz = F_MUL(f, dz, dy, a->x);	/* %$z' = x z^2$% */
+ 
+     dy = F_SQR(f, dy, a->x);		/* %$x^2$% */
+     u = F_MUL(f, MP_NEW, a->y, a->z);	/* %$y z$% */
+     u = F_ADD(f, u, u, dz);		/* %$z' + y z$% */
+     u = F_ADD(f, u, u, dy);		/* %$u = z' + x^2 + y z$% */
+ 
+     v = F_SQR(f, MP_NEW, dy);		/* %$x^4$% */
+     dy = F_MUL(f, dy, v, dz);		/* %$x^4 z'$% */
+     v = F_MUL(f, v, u, dx);		/* %$u x'$% */
+     dy = F_ADD(f, dy, dy, v);		/* %$y' = x^4 z' + u x'$% */
+ 
+     EC_DESTROY(d);
+     d->x = dx;
+     d->y = dy;
+     d->z = dz;
+     MP_DROP(u);
+     MP_DROP(v);
+     assert(!(d->x->f & MP_DESTROYED));
+     assert(!(d->y->f & MP_DESTROYED));
+     assert(!(d->z->f & MP_DESTROYED));
+   }
+   return (d);
+ }
+ 
+ static ec *ecadd(ec_curve *c, ec *d, const ec *a, const ec *b)
+ {
+   if (a == b)
+     ecdbl(c, d, a);
+   else if (EC_ATINF(a))
+     EC_COPY(d, b);
+   else if (EC_ATINF(b))
+     EC_COPY(d, a);
+   else {
+     field *f = c->f;
+     ecctx *cc = (ecctx *)c;
+     mp *lambda;
+     mp *dx, *dy;
+ 
+     if (!MP_EQ(a->x, b->x)) {
+       dx = F_ADD(f, MP_NEW, a->x, b->x); /* %$x_0 + x_1$% */
+       dy = F_INV(f, MP_NEW, dx);	/* %$(x_0 + x_1)^{-1}$% */
+       dx = F_ADD(f, dx, a->y, b->y);	/* %$y_0 + y_1$% */
+       lambda = F_MUL(f, MP_NEW, dy, dx);
+ 				  /* %$\lambda = (y_0 + y_1)/(x_0 + x_1)$% */
+ 
+       dx = F_SQR(f, dx, lambda);	/* %$\lambda^2$% */
+       dx = F_ADD(f, dx, dx, lambda);	/* %$\lambda^2 + \lambda$% */
+       dx = F_ADD(f, dx, dx, cc->a);	/* %$a + \lambda^2 + \lambda$% */
+       dx = F_ADD(f, dx, dx, a->x);    /* %$a + \lambda^2 + \lambda + x_0$% */
+       dx = F_ADD(f, dx, dx, b->x);
+                            /* %$x' = a + \lambda^2 + \lambda + x_0 + x_1$% */
+     } else if (!MP_EQ(a->y, b->y) || F_ZEROP(f, a->x)) {
+       EC_SETINF(d);
+       return (d);
+     } else {
+       dx = F_INV(f, MP_NEW, a->x);	/* %$x^{-1}$% */
+       dy = F_MUL(f, MP_NEW, dx, a->y);	/* %$y/x$% */
+       lambda = F_ADD(f, dy, dy, a->x);	/* %$\lambda = x + y/x$% */
+ 
+       dx = F_SQR(f, dx, lambda);	/* %$\lambda^2$% */
+       dx = F_ADD(f, dx, dx, lambda);	/* %$\lambda^2 + \lambda$% */
+       dx = F_ADD(f, dx, dx, cc->a);    /* %$x' = a + \lambda^2 + \lambda$% */
+       dy = MP_NEW;
+     }
+       
+     dy = F_ADD(f, dy, a->x, dx);	/* %$ x + x' $% */
+     dy = F_MUL(f, dy, dy, lambda);	/* %$ (x + x') \lambda$% */
+     dy = F_ADD(f, dy, dy, a->y);	/* %$ (x + x') \lambda + y$% */
+     dy = F_ADD(f, dy, dy, dx);	    /* %$ y' = (x + x') \lambda + y + x'$% */
+ 
+     EC_DESTROY(d);
+     d->x = dx;
+     d->y = dy;
+     d->z = 0;
+     MP_DROP(lambda);
+   }
+   return (d);
+ }
+ 
+ static ec *ecprojadd(ec_curve *c, ec *d, const ec *a, const ec *b)
+ {
+   if (a == b)
+     c->ops->dbl(c, d, a);
+   else if (EC_ATINF(a))
+     EC_COPY(d, b);
+   else if (EC_ATINF(b))
+     EC_COPY(d, a);
+   else {
+     field *f = c->f;
+     ecctx *cc = (ecctx *)c;
+     mp *dx, *dy, *dz, *u, *uu, *v, *t, *s, *ss, *r, *w, *l;
+ 
+     dz = F_SQR(f, MP_NEW, b->z);	/* %$z_1^2$% */
+     u = F_MUL(f, MP_NEW, dz, a->x);	/* %$u_0 = x_0 z_1^2$% */
+     t = F_MUL(f, MP_NEW, dz, b->z);	/* %$z_1^3$% */
+     s = F_MUL(f, MP_NEW, t, a->y);	/* %$s_0 = y_0 z_1^3$% */
+ 
+     dz = F_SQR(f, dz, a->z);		/* %$z_0^2$% */
+     uu = F_MUL(f, MP_NEW, dz, b->x);	/* %$u_1 = x_1 z_0^2$% */
+     t = F_MUL(f, t, dz, a->z);		/* %$z_0^3$% */
+     ss = F_MUL(f, MP_NEW, t, b->y);	/* %$s_1 = y_1 z_0^3$% */
+ 
+     w = F_ADD(f, u, u, uu);		/* %$r = u_0 + u_1$% */
+     r = F_ADD(f, s, s, ss);		/* %$w = s_0 + s_1$% */
+     if (F_ZEROP(f, w)) {
+       MP_DROP(w);
+       MP_DROP(uu);
+       MP_DROP(ss);
+       MP_DROP(t);
+       MP_DROP(dz);
+       if (F_ZEROP(f, r)) {
+ 	MP_DROP(r);
+ 	return (c->ops->dbl(c, d, a));
+       } else {
+ 	MP_DROP(r);
+ 	EC_SETINF(d);
+ 	return (d);
+       }
+     }
+ 
+     l = F_MUL(f, t, a->z, w);		/* %$l = z_0 w$% */
+ 
+     dz = F_MUL(f, dz, l, b->z);		/* %$z' = l z_1$% */
+ 
+     ss = F_MUL(f, ss, r, b->x);		/* %$r x_1$% */
+     t = F_MUL(f, uu, l, b->y);		/* %$l y_1$% */
+     v = F_ADD(f, ss, ss, t);		/* %$v = r x_1 + l y_1$% */
+ 
+     t = F_ADD(f, t, r, dz);		/* %$t = r + z'$% */
+ 
+     uu = F_SQR(f, MP_NEW, dz);		/* %$z'^2$% */
+     dx = F_MUL(f, MP_NEW, uu, cc->a);	/* %$a z'^2$% */
+     uu = F_MUL(f, uu, t, r);		/* %$t r$% */
+     dx = F_ADD(f, dx, dx, uu);		/* %$a z'^2 + t r$% */
+     r = F_SQR(f, r, w);			/* %$w^2$% */
+     uu = F_MUL(f, uu, r, w);		/* %$w^3$% */
+     dx = F_ADD(f, dx, dx, uu);		/* %$x' = a z'^2 + t r + w^3$% */
+ 
+     r = F_SQR(f, r, l);			/* %$l^2$% */
+     dy = F_MUL(f, uu, v, r);		/* %$v l^2$% */
+     l = F_MUL(f, l, t, dx);		/* %$t x'$% */
+     dy = F_ADD(f, dy, dy, l);		/* %$y' = t x' + v l^2$% */
+ 
+     EC_DESTROY(d);
+     d->x = dx;
+     d->y = dy;
+     d->z = dz;
+     MP_DROP(l);
+     MP_DROP(r);
+     MP_DROP(w);
+     MP_DROP(t);
+     MP_DROP(v);
+   }
+   return (d);
+ }
+ 
+ static int eccheck(ec_curve *c, const ec *p)
+ {
+   ecctx *cc = (ecctx *)c;
+   field *f = c->f;
+   int rc;
+   mp *u, *v;
+ 
+   v = F_SQR(f, MP_NEW, p->x);
+   u = F_MUL(f, MP_NEW, v, p->x);
+   v = F_MUL(f, v, v, cc->a);
+   u = F_ADD(f, u, u, v);
+   u = F_ADD(f, u, u, cc->b);
+   v = F_MUL(f, v, p->x, p->y);
+   u = F_ADD(f, u, u, v);
+   v = F_SQR(f, v, p->y);
+   u = F_ADD(f, u, u, v);
+   rc = F_ZEROP(f, u);
+   mp_drop(u);
+   mp_drop(v);
+   return (rc);
+ }
+ 
+ static int ecprojcheck(ec_curve *c, const ec *p)
+ {
+   ec t = EC_INIT;
+   int rc;
+   
+   c->ops->fix(c, &t, p);
+   rc = eccheck(c, &t);
+   EC_DESTROY(&t);
+   return (rc);
+ }
+ 
+ static void ecdestroy(ec_curve *c)
+ {
+   ecctx *cc = (ecctx *)c;
+   MP_DROP(cc->a);
+   MP_DROP(cc->b);
+   if (cc->bb) MP_DROP(cc->bb);
+   DESTROY(cc);
+ }
+ 
+ /* --- @ec_bin@, @ec_binproj@ --- *
+  *
+  * Arguments:	@field *f@ = the underlying field for this elliptic curve
+  *		@mp *a, *b@ = the coefficients for this curve
+  *
+  * Returns:	A pointer to the curve.
+  *
+  * Use:		Creates a curve structure for an elliptic curve defined over
+  *		a binary field.  The @binproj@ variant uses projective
+  *		coordinates, which can be a win.
+  */
+ 
+ ec_curve *ec_bin(field *f, mp *a, mp *b)
+ {
+   ecctx *cc = CREATE(ecctx);
+   cc->c.ops = &ec_binops;
+   cc->c.f = f;
+   cc->a = F_IN(f, MP_NEW, a);
+   cc->b = F_IN(f, MP_NEW, b);
+   cc->bb = 0;
+   return (&cc->c);
+ }
+ 
+ ec_curve *ec_binproj(field *f, mp *a, mp *b)
+ {
+   ecctx *cc = CREATE(ecctx);
+   cc->c.ops = &ec_binprojops;
+   cc->c.f = f;
+   cc->a = F_IN(f, MP_NEW, a);
+   cc->b = F_IN(f, MP_NEW, b);
+   cc->bb = F_SQRT(f, MP_NEW, b);
+   cc->bb = F_SQRT(f, cc->bb, cc->bb);
+   return (&cc->c);
+ }
+ 
+ static const ec_ops ec_binops = {
+   ecdestroy, ec_idin, ec_idout, ec_idfix,
+   0, ecneg, ecadd, ec_stdsub, ecdbl, eccheck
+ };
+ 
+ static const ec_ops ec_binprojops = {
+   ecdestroy, ec_projin, ec_projout, ec_projfix,
+   0, ecprojneg, ecprojadd, ec_stdsub, ecprojdbl, ecprojcheck
+ };
+ 
+ /*----- Test rig ----------------------------------------------------------*/
+ 
+ #ifdef TEST_RIG
+ 
+ #define MP(x) mp_readstring(MP_NEW, #x, 0, 0)
+ 
+ int main(int argc, char *argv[])
+ {
+   field *f;
+   ec_curve *c;
+   ec g = EC_INIT, d = EC_INIT;
+   mp *p, *a, *b, *r;
+   int i, n = argc == 1 ? 1 : atoi(argv[1]);
+ 
+   printf("ec-bin: ");
+   fflush(stdout);
+   a = MP(1);
+   b = MP(0x066647ede6c332c7f8c0923bb58213b333b20e9ce4281fe115f7d8f90ad);
+   p = MP(0x20000000000000000000000000000000000000004000000000000000001);
+   r =
+   MP(6901746346790563787434755862277025555839812737345013555379383634485462);
+ 
+   f = field_binpoly(p);
+   c = ec_binproj(f, a, b);
+   
+   g.x = MP(0x0fac9dfcbac8313bb2139f1bb755fef65bc391f8b36f8f8eb7371fd558b);
+   g.y = MP(0x1006a08a41903350678e58528bebf8a0beff867a7ca36716f7e01f81052);
+ 
+   for (i = 0; i < n; i++) { 
+     ec_mul(c, &d, &g, r);
+     if (EC_ATINF(&d)) {
+       fprintf(stderr, "zero too early\n");
+       return (1);
+     }
+     ec_add(c, &d, &d, &g);
+     if (!EC_ATINF(&d)) {
+       fprintf(stderr, "didn't reach zero\n");
+       MP_EPRINTX("d.x", d.x);
+       MP_EPRINTX("d.y", d.y);
+       MP_EPRINTX("d.z", d.y);
+       return (1);
+     }
+     ec_destroy(&d);
+   }
+ 
+   ec_destroy(&g);
+   ec_destroycurve(c);
+   F_DESTROY(f);
+   MP_DROP(p); MP_DROP(a); MP_DROP(b); MP_DROP(r);
+   assert(!mparena_count(&mparena_global));
+   printf("ok\n");
+   return (0);
+ }
+ 
+ #endif
+ 
+ /*----- That's all, folks -------------------------------------------------*/
diff --cc ec-exp.h
index 9ba7bc5,bb4e08a..0daf717
--- a/ec-exp.h
+++ b/ec-exp.h
@@@ -1,6 -1,6 +1,6 @@@
  /* -*-c-*-
   *
-  * $Id: ec-exp.h,v 1.2 2003/05/15 23:25:59 mdw Exp $
 - * $Id: ec-exp.h,v 1.2.4.1 2004/03/20 00:13:31 mdw Exp $
++ * $Id: ec-exp.h,v 1.3 2004/03/21 22:52:06 mdw Exp $
   *
   * Exponentiation operations for elliptic curves
   *
@@@ -30,6 -30,9 +30,12 @@@
  /*----- Revision history --------------------------------------------------* 
   *
   * $Log: ec-exp.h,v $
++ * Revision 1.3  2004/03/21 22:52:06  mdw
++ * Merge and close elliptic curve branch.
++ *
+  * Revision 1.2.4.1  2004/03/20 00:13:31  mdw
+  * Projective coordinates for prime curves
+  *
   * Revision 1.2  2003/05/15 23:25:59  mdw
   * Make elliptic curve stuff build.
   *
diff --cc ec-prime.c
index 4611855,40f487e..bdc6368
--- a/ec-prime.c
+++ b/ec-prime.c
@@@ -1,6 -1,6 +1,6 @@@
  /* -*-c-*-
   *
-  * $Id: ec-prime.c,v 1.3 2003/05/15 23:25:59 mdw Exp $
 - * $Id: ec-prime.c,v 1.3.4.3 2004/03/21 22:39:46 mdw Exp $
++ * $Id: ec-prime.c,v 1.4 2004/03/21 22:52:06 mdw Exp $
   *
   * Elliptic curves over prime fields
   *
@@@ -30,6 -30,15 +30,18 @@@
  /*----- Revision history --------------------------------------------------* 
   *
   * $Log: ec-prime.c,v $
++ * Revision 1.4  2004/03/21 22:52:06  mdw
++ * Merge and close elliptic curve branch.
++ *
+  * Revision 1.3.4.3  2004/03/21 22:39:46  mdw
+  * Elliptic curves on binary fields work.
+  *
+  * Revision 1.3.4.2  2004/03/20 00:13:31  mdw
+  * Projective coordinates for prime curves
+  *
+  * Revision 1.3.4.1  2003/06/10 13:43:53  mdw
+  * Simple (non-projective) curves over prime fields now seem to work.
+  *
   * Revision 1.3  2003/05/15 23:25:59  mdw
   * Make elliptic curve stuff build.
   *
diff --cc ec.c
index 47c01a9,a2b229f..c95333f
--- a/ec.c
+++ b/ec.c
@@@ -1,6 -1,6 +1,6 @@@
  /* -*-c-*-
   *
-  * $Id: ec.c,v 1.4 2003/05/15 23:25:59 mdw Exp $
 - * $Id: ec.c,v 1.4.4.2 2004/03/20 00:13:31 mdw Exp $
++ * $Id: ec.c,v 1.5 2004/03/21 22:52:06 mdw Exp $
   *
   * Elliptic curve definitions
   *
@@@ -30,6 -30,12 +30,15 @@@
  /*----- Revision history --------------------------------------------------* 
   *
   * $Log: ec.c,v $
++ * Revision 1.5  2004/03/21 22:52:06  mdw
++ * Merge and close elliptic curve branch.
++ *
+  * Revision 1.4.4.2  2004/03/20 00:13:31  mdw
+  * Projective coordinates for prime curves
+  *
+  * Revision 1.4.4.1  2003/06/10 13:43:53  mdw
+  * Simple (non-projective) curves over prime fields now seem to work.
+  *
   * Revision 1.4  2003/05/15 23:25:59  mdw
   * Make elliptic curve stuff build.
   *
diff --cc ec.h
index 105838c,bb5e201..07f1468
--- a/ec.h
+++ b/ec.h
@@@ -1,6 -1,6 +1,6 @@@
  /* -*-c-*-
   *
-  * $Id: ec.h,v 1.4 2003/05/15 23:25:59 mdw Exp $
 - * $Id: ec.h,v 1.4.4.3 2004/03/21 22:39:46 mdw Exp $
++ * $Id: ec.h,v 1.5 2004/03/21 22:52:06 mdw Exp $
   *
   * Elliptic curve definitions
   *
@@@ -30,6 -30,15 +30,18 @@@
  /*----- Revision history --------------------------------------------------* 
   *
   * $Log: ec.h,v $
++ * Revision 1.5  2004/03/21 22:52:06  mdw
++ * Merge and close elliptic curve branch.
++ *
+  * Revision 1.4.4.3  2004/03/21 22:39:46  mdw
+  * Elliptic curves on binary fields work.
+  *
+  * Revision 1.4.4.2  2004/03/20 00:13:31  mdw
+  * Projective coordinates for prime curves
+  *
+  * Revision 1.4.4.1  2003/06/10 13:43:53  mdw
+  * Simple (non-projective) curves over prime fields now seem to work.
+  *
   * Revision 1.4  2003/05/15 23:25:59  mdw
   * Make elliptic curve stuff build.
   *
diff --cc exp.h
index 6cfdfd8,fc9e3a9..6bfd686
--- a/exp.h
+++ b/exp.h
@@@ -1,6 -1,6 +1,6 @@@
  /* -*-c-*-
   *
-  * $Id: exp.h,v 1.1 2001/06/16 13:00:59 mdw Exp $
 - * $Id: exp.h,v 1.1.4.1 2004/03/20 00:13:31 mdw Exp $
++ * $Id: exp.h,v 1.2 2004/03/21 22:52:06 mdw Exp $
   *
   * Generalized exponentiation
   *
@@@ -30,6 -30,9 +30,12 @@@
  /*----- Revision history --------------------------------------------------* 
   *
   * $Log: exp.h,v $
++ * Revision 1.2  2004/03/21 22:52:06  mdw
++ * Merge and close elliptic curve branch.
++ *
+  * Revision 1.1.4.1  2004/03/20 00:13:31  mdw
+  * Projective coordinates for prime curves
+  *
   * Revision 1.1  2001/06/16 13:00:59  mdw
   * New generic exponentation code.  Includes sliding-window simultaneous
   * exponentiation.
diff --cc f-binpoly.c
index 0000000,509efc4..02e683d
mode 000000,100644..100644
--- a/f-binpoly.c
+++ b/f-binpoly.c
@@@ -1,0 -1,142 +1,145 @@@
+ /* -*-c-*-
+  *
 - * $Id: f-binpoly.c,v 1.1.2.1 2004/03/21 22:39:46 mdw Exp $
++ * $Id: f-binpoly.c,v 1.2 2004/03/21 22:52:06 mdw Exp $
+  *
+  * Binary fields with polynomial basis representation
+  *
+  * (c) 2004 Straylight/Edgeware
+  */
+ 
+ /*----- Licensing notice --------------------------------------------------* 
+  *
+  * This file is part of Catacomb.
+  *
+  * Catacomb is free software; you can redistribute it and/or modify
+  * 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: f-binpoly.c,v $
++ * Revision 1.2  2004/03/21 22:52:06  mdw
++ * Merge and close elliptic curve branch.
++ *
+  * Revision 1.1.2.1  2004/03/21 22:39:46  mdw
+  * Elliptic curves on binary fields work.
+  *
+  */
+ 
+ /*----- Header files ------------------------------------------------------*/
+ 
+ #include <mLib/sub.h>
+ 
+ #include "field.h"
+ #include "gf.h"
+ #include "gfreduce.h"
+ 
+ /*----- Data structures ---------------------------------------------------*/
+ 
+ typedef struct fctx {
+   field f;
+   gfreduce r;
+ } fctx;
+ 
+ /*----- Main code ---------------------------------------------------------*/
+ 
+ /* --- Field operations --- */
+ 
+ static void fdestroy(field *ff)
+ {
+   fctx *f = (fctx *)ff;
+   gfreduce_destroy(&f->r);
+   DESTROY(f);
+ }
+ 
+ static int fzerop(field *ff, mp *x)
+ {
+   return (!MP_LEN(x));
+ }
+ 
+ static mp *fadd(field *ff, mp *d, mp *x, mp *y)
+ {
+   return (gf_add(d, x, y));
+ }
+ 
+ static mp *fmul(field *ff, mp *d, mp *x, mp *y)
+ {
+   fctx *f = (fctx *)ff;
+   d = gf_mul(d, x, y);
+   return (gfreduce_do(&f->r, d, d));
+ }
+ 
+ static mp *fsqr(field *ff, mp *d, mp *x)
+ {
+   fctx *f = (fctx *)ff;
+   d = gf_sqr(d, x);
+   return (gfreduce_do(&f->r, d, d));
+ }
+ 
+ static mp *finv(field *ff, mp *d, mp *x)
+ {
+   fctx *f = (fctx *)ff;
+   gf_gcd(0, 0, &d, f->r.p, x);
+   return (d);
+ }
+ 
+ static mp *freduce(field *ff, mp *d, mp *x)
+ {
+   fctx *f = (fctx *)ff;
+   return (gfreduce_do(&f->r, d, x));
+ }
+ 
+ static mp *fsqrt(field *ff, mp *d, mp *x)
+ {
+   fctx *f = (fctx *)ff;
+   return (gfreduce_sqrt(&f->r, d, x));
+ }
+ 
+ static mp *fquadsolve(field *ff, mp *d, mp *x)
+ {
+   fctx *f = (fctx *)ff;
+   return (gfreduce_quadsolve(&f->r, d, x));
+ }
+ 
+ /* --- Field operations table --- */
+ 
+ static field_ops fops = {
+   fdestroy,
+   freduce, field_id,
+   fzerop, field_id, fadd, fadd, fmul, fsqr, finv, freduce, fsqrt,
+   fquadsolve,
+   0, 0, 0, 0
+ };
+ 
+ /* --- @field_binpoly@ --- *
+  *
+  * Arguments:	@mp *p@ = the reduction polynomial
+  *
+  * Returns:	A pointer to the field.
+  *
+  * Use:		Creates a field structure for a binary field mod @p@.
+  */
+ 
+ field *field_binpoly(mp *p)
+ {
+   fctx *f = CREATE(fctx);
+   f->f.ops = &fops;
+   f->f.zero = MP_ZERO;
+   f->f.one = MP_ONE;
+   gfreduce_create(&f->r, p);
+   return (&f->f);
+ }
+ 
+ /*----- That's all, folks -------------------------------------------------*/
diff --cc f-prime.c
index f4bf3a7,7215ec8..8454902
--- a/f-prime.c
+++ b/f-prime.c
@@@ -1,6 -1,6 +1,6 @@@
  /* -*-c-*-
   *
-  * $Id: f-prime.c,v 1.3 2003/05/15 23:25:59 mdw Exp $
 - * $Id: f-prime.c,v 1.3.4.3 2004/03/21 22:39:46 mdw Exp $
++ * $Id: f-prime.c,v 1.4 2004/03/21 22:52:06 mdw Exp $
   *
   * Prime fields with Montgomery arithmetic
   *
@@@ -30,6 -30,15 +30,18 @@@
  /*----- Revision history --------------------------------------------------* 
   *
   * $Log: f-prime.c,v $
++ * Revision 1.4  2004/03/21 22:52:06  mdw
++ * Merge and close elliptic curve branch.
++ *
+  * Revision 1.3.4.3  2004/03/21 22:39:46  mdw
+  * Elliptic curves on binary fields work.
+  *
+  * Revision 1.3.4.2  2004/03/20 00:13:31  mdw
+  * Projective coordinates for prime curves
+  *
+  * Revision 1.3.4.1  2003/06/10 13:43:53  mdw
+  * Simple (non-projective) curves over prime fields now seem to work.
+  *
   * Revision 1.3  2003/05/15 23:25:59  mdw
   * Make elliptic curve stuff build.
   *
diff --cc field.c
index c7e7559,0f02fa7..c310b93
--- a/field.c
+++ b/field.c
@@@ -1,8 -1,9 +1,9 @@@
  /* -*-c-*-
   *
-  * $Id: field.c,v 1.1 2001/05/07 17:30:13 mdw Exp $
 - * $Id: field.c,v 1.1.4.1 2003/06/10 13:43:53 mdw Exp $
++ * $Id: field.c,v 1.2 2004/03/21 22:52:06 mdw Exp $
+  *
+  * Abstract field operations
   *
-  * [Abstract field operations *
   * (c) 2001 Straylight/Edgeware
   */
  
@@@ -29,6 -30,9 +30,12 @@@
  /*----- Revision history --------------------------------------------------* 
   *
   * $Log: field.c,v $
++ * Revision 1.2  2004/03/21 22:52:06  mdw
++ * Merge and close elliptic curve branch.
++ *
+  * Revision 1.1.4.1  2003/06/10 13:43:53  mdw
+  * Simple (non-projective) curves over prime fields now seem to work.
+  *
   * Revision 1.1  2001/05/07 17:30:13  mdw
   * Add an internal-representation no-op function.
   *
diff --cc field.h
index 5a70649,ea019c5..dd674c9
--- a/field.h
+++ b/field.h
@@@ -1,6 -1,6 +1,6 @@@
  /* -*-c-*-
   *
-  * $Id: field.h,v 1.3 2002/01/13 13:48:44 mdw Exp $
 - * $Id: field.h,v 1.3.4.2 2004/03/21 22:39:46 mdw Exp $
++ * $Id: field.h,v 1.4 2004/03/21 22:52:06 mdw Exp $
   *
   * Definitions for field arithmetic
   *
@@@ -30,6 -30,12 +30,15 @@@
  /*----- Revision history --------------------------------------------------* 
   *
   * $Log: field.h,v $
++ * Revision 1.4  2004/03/21 22:52:06  mdw
++ * Merge and close elliptic curve branch.
++ *
+  * Revision 1.3.4.2  2004/03/21 22:39:46  mdw
+  * Elliptic curves on binary fields work.
+  *
+  * Revision 1.3.4.1  2004/03/20 00:13:31  mdw
+  * Projective coordinates for prime curves
+  *
   * Revision 1.3  2002/01/13 13:48:44  mdw
   * Further progress.
   *
diff --cc gf-arith.c
index 0000000,6838e44..debfba1
mode 000000,100644..100644
--- a/gf-arith.c
+++ b/gf-arith.c
@@@ -1,0 -1,263 +1,266 @@@
+ /* -*-c-*-
+  *
 - * $Id: gf-arith.c,v 1.1.2.1 2004/03/21 22:39:46 mdw Exp $
++ * $Id: gf-arith.c,v 1.2 2004/03/21 22:52:06 mdw Exp $
+  *
+  * Basic arithmetic on binary polynomials
+  *
+  * (c) 2004 Straylight/Edgeware
+  */
+ 
+ /*----- Licensing notice --------------------------------------------------* 
+  *
+  * This file is part of Catacomb.
+  *
+  * Catacomb is free software; you can redistribute it and/or modify
+  * 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: gf-arith.c,v $
++ * Revision 1.2  2004/03/21 22:52:06  mdw
++ * Merge and close elliptic curve branch.
++ *
+  * Revision 1.1.2.1  2004/03/21 22:39:46  mdw
+  * Elliptic curves on binary fields work.
+  *
+  */
+ 
+ /*----- Header files ------------------------------------------------------*/
+ 
+ #include "gf.h"
+ 
+ /*----- Macros ------------------------------------------------------------*/
+ 
+ #define MAX(x, y) ((x) >= (y) ? (x) : (y))
+ 
+ /*----- Main code ---------------------------------------------------------*/
+ 
+ /* --- @gf_add@ --- *
+  *
+  * Arguments:	@mp *d@ = destination
+  *		@mp *a, *b@ = sources
+  *
+  * Returns:	Result, @a@ added to @b@.
+  */
+ 
+ mp *gf_add(mp *d, mp *a, mp *b)
+ {
+   MP_DEST(d, MAX(MP_LEN(a), MP_LEN(b)), (a->f | b->f) & MP_BURN);
+   gfx_add(d->v, d->vl, a->v, a->vl, b->v, b->vl);
+   d->f = (a->f | b->f) & MP_BURN;
+   MP_SHRINK(d);
+   return (d);
+ }
+ 
+ /* --- @gf_mul@ --- *
+  *
+  * Arguments:	@mp *d@ = destination
+  *		@mp *a, *b@ = sources
+  *
+  * Returns:	Result, @a@ multiplied by @b@.
+  */
+ 
+ mp *gf_mul(mp *d, mp *a, mp *b)
+ {
+   a = MP_COPY(a);
+   b = MP_COPY(b);
+ 
+   if (MP_LEN(a) <= MPK_THRESH || MP_LEN(b) <= GFK_THRESH) {
+     MP_DEST(d, MP_LEN(a) + MP_LEN(b), a->f | b->f | MP_UNDEF);
+     gfx_mul(d->v, d->vl, a->v, a->vl, b->v, b->vl);
+   } else {
+     size_t m = MAX(MP_LEN(a), MP_LEN(b));
+     mpw *s;
+     MP_DEST(d, 2 * m, a->f | b->f | MP_UNDEF);
+     s = mpalloc(d->a, 2 * m);
+     gfx_kmul(d->v, d->vl, a->v, a->vl, b->v, b->vl, s, s + 2 * m);
+     mpfree(d->a, s);
+   }
+ 
+   d->f = (a->f | b->f) & MP_BURN;
+   MP_SHRINK(d);
+   MP_DROP(a);
+   MP_DROP(b);
+   return (d);
+ }
+ 
+ /* --- @gf_sqr@ --- *
+  *
+  * Arguments:	@mp *d@ = destination
+  *		@mp *a@ = source
+  *
+  * Returns:	Result, @a@ squared.
+  */
+ 
+ mp *gf_sqr(mp *d, mp *a)
+ {
+   MP_COPY(a);
+   MP_DEST(d, 2 * MP_LEN(a), a->f & MP_BURN);
+   gfx_sqr(d->v, d->vl, a->v, a->vl);
+   d->f = a->f & MP_BURN;
+   MP_SHRINK(d);
+   MP_DROP(a);
+   return (d);
+ }
+ 
+ /* --- @gf_div@ --- *
+  *
+  * Arguments:	@mp **qq, **rr@ = destination, quotient and remainder
+  *		@mp *a, *b@ = sources
+  *
+  * Use:		Calculates the quotient and remainder when @a@ is divided by
+  *		@b@.  The destinations @*qq@ and @*rr@ must be distinct.
+  *		Either of @qq@ or @rr@ may be null to indicate that the
+  *		result is irrelevant.  (Discarding both results is silly.)
+  *		There is a performance advantage if @a == *rr@.
+  */
+ 
+ void gf_div(mp **qq, mp **rr, mp *a, mp *b)
+  {
+   mp *r = rr ? *rr : MP_NEW;
+   mp *q = qq ? *qq : MP_NEW;
+ 
+   /* --- Set the remainder up right --- */
+ 
+   b = MP_COPY(b);
+   a = MP_COPY(a);
+   if (r)
+     MP_DROP(r);
+   r = a;
+   MP_DEST(r, MP_LEN(b) + 2, a->f | b->f);
+ 
+   /* --- Fix up the quotient too --- */
+ 
+   r = MP_COPY(r);
+   MP_DEST(q, MP_LEN(r), r->f | MP_UNDEF);
+   MP_DROP(r);
+ 
+   /* --- Perform the calculation --- */
+ 
+   gfx_div(q->v, q->vl, r->v, r->vl, b->v, b->vl);
+ 
+   /* --- Sort out the sign of the results --- *
+    *
+    * If the signs of the arguments differ, and the remainder is nonzero, I
+    * must add one to the absolute value of the quotient and subtract the
+    * remainder from @b@.
+    */
+ 
+   q->f = (r->f | b->f) & MP_BURN;
+   r->f = (r->f | b->f) & MP_BURN;
+ 
+   /* --- Store the return values --- */
+ 
+   MP_DROP(b);
+ 
+   if (!qq)
+     MP_DROP(q);
+   else {
+     MP_SHRINK(q);
+     *qq = q;
+   }
+ 
+   if (!rr)
+     MP_DROP(r);
+   else {
+     MP_SHRINK(r);
+     *rr = r;
+   }
+ }
+ 
+ /*----- Test rig ----------------------------------------------------------*/
+ 
+ #ifdef TEST_RIG
+ 
+ static int verify(const char *op, mp *expect, mp *result, mp *a, mp *b)
+ {
+   if (!MP_EQ(expect, result)) {
+     fprintf(stderr, "\n*** %s failed", op);
+     fputs("\n*** a      = ", stderr); mp_writefile(a, stderr, 16);
+     fputs("\n*** b      = ", stderr); mp_writefile(b, stderr, 16);
+     fputs("\n*** result = ", stderr); mp_writefile(result, stderr, 16);
+     fputs("\n*** expect = ", stderr); mp_writefile(expect, stderr, 16);
+     fputc('\n', stderr);
+     return (0);
+   }
+   return (1);
+ }
+ 
+ #define RIG(name, op)							\
+   static int t##name(dstr *v)						\
+   {									\
+     mp *a = *(mp **)v[0].buf;						\
+     mp *b = *(mp **)v[1].buf;						\
+     mp *r = *(mp **)v[2].buf;						\
+     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);							\
+   }
+ 
+ RIG(add, gf_add)
+ RIG(mul, gf_mul)
+ 
+ #undef RIG
+ 
+ static int tsqr(dstr *v)
+ {
+   mp *a = *(mp **)v[0].buf;
+   mp *r = *(mp **)v[1].buf;
+   mp *c = MP_NEW;
+   int ok = 1;
+   c = gf_sqr(MP_NEW, a);
+   ok &= verify("sqr", r, c, a, MP_ZERO);
+   mp_drop(a); mp_drop(r); mp_drop(c);
+   assert(mparena_count(MPARENA_GLOBAL) == 0);
+   return (ok);
+ }
+ 
+ static int tdiv(dstr *v)
+ {
+   mp *a = *(mp **)v[0].buf;
+   mp *b = *(mp **)v[1].buf;
+   mp *q = *(mp **)v[2].buf;
+   mp *r = *(mp **)v[3].buf;
+   mp *c = MP_NEW, *d = MP_NEW;
+   int ok = 1;
+   gf_div(&c, &d, a, b);
+   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);
+ }
+ 
+ static test_chunk tests[] = {
+   { "add", tadd, { &type_mp, &type_mp, &type_mp, 0 } },
+   { "mul", tmul, { &type_mp, &type_mp, &type_mp, 0 } },
+   { "sqr", tsqr, { &type_mp, &type_mp, 0 } },
+   { "div", tdiv, { &type_mp, &type_mp, &type_mp, &type_mp, 0 } },
+   { 0, 0, { 0 } },
+ };
+ 
+ int main(int argc, char *argv[])
+ {
+   sub_init();
+   test_run(argc, argv, tests, SRCDIR "/tests/gf");
+   return (0);
+ }
+ 
+ #endif
+ 
+ /*----- That's all, folks -------------------------------------------------*/
diff --cc gf-gcd.c
index 0000000,64d61d4..7c09d3a
mode 000000,100644..100644
--- a/gf-gcd.c
+++ b/gf-gcd.c
@@@ -1,0 -1,259 +1,262 @@@
+ /* -*-c-*-
+  *
 - * $Id: gf-gcd.c,v 1.1.2.1 2004/03/21 22:39:46 mdw Exp $
++ * $Id: gf-gcd.c,v 1.2 2004/03/21 22:52:06 mdw Exp $
+  *
+  * Euclidian algorithm on binary polynomials
+  *
+  * (c) 2004 Straylight/Edgeware
+  */
+ 
+ /*----- Licensing notice --------------------------------------------------* 
+  *
+  * This file is part of Catacomb.
+  *
+  * Catacomb is free software; you can redistribute it and/or modify
+  * 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: gf-gcd.c,v $
++ * Revision 1.2  2004/03/21 22:52:06  mdw
++ * Merge and close elliptic curve branch.
++ *
+  * Revision 1.1.2.1  2004/03/21 22:39:46  mdw
+  * Elliptic curves on binary fields work.
+  *
+  */
+ 
+ /*----- Header files ------------------------------------------------------*/
+ 
+ #include "gf.h"
+ 
+ /*----- Main code ---------------------------------------------------------*/
+ 
+ /* --- @gf_gcd@ --- *
+  *
+  * Arguments:	@mp **gcd, **xx, **yy@ = where to write the results
+  *		@mp *a, *b@ = sources (must be nonzero)
+  *
+  *
+  * Returns:	---
+  *
+  * 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.
+  */
+ 
+ void gf_gcd(mp **gcd, mp **xx, mp **yy, mp *a, mp *b)
+ {
+   mp *x = MP_ONE, *X = MP_ZERO;
+   mp *y = MP_ZERO, *Y = MP_ONE;
+   mp *u, *v;
+   mp *q = MP_NEW;
+   unsigned f = 0;
+ 
+ #define f_swap 1u
+ #define f_ext 2u
+ 
+   /* --- Sort out some initial flags --- */
+ 
+   if (xx || yy)
+     f |= f_ext;
+ 
+   /* --- Ensure that @a@ is larger than @b@ --- *
+    *
+    * If they're the same length we don't care which order they're in, so this
+    * unsigned comparison is fine.
+    */
+ 
+   if (MPX_UCMP(a->v, a->vl, <, b->v, b->vl)) {
+     { mp *t = a; a = b; b = t; }
+     f |= f_swap;
+   }
+ 
+   /* --- Check for zeroness --- */
+ 
+   if (MP_EQ(b, MP_ZERO)) {
+ 
+     /* --- Store %$|a|$% as the GCD --- */
+ 
+     if (gcd) {
+       if (*gcd) MP_DROP(*gcd);
+       a = MP_COPY(a);
+       *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_EQ(a, MP_ZERO))
+ 	  *xx = MP_ZERO;
+ 	else
+ 	  *xx = MP_ONE;
+       }
+       if (yy) {
+ 	if (*yy) MP_DROP(*yy);
+ 	*yy = MP_ZERO;
+       }
+     }
+     return;
+   }
+ 
+   /* --- Take a reference to the arguments --- */
+ 
+   a = MP_COPY(a);
+   b = MP_COPY(b);
+ 
+   /* --- Make sure @a@ and @b@ are not both even --- */
+ 
+   MP_SPLIT(a); a->f &= ~MP_NEG;
+   MP_SPLIT(b); b->f &= ~MP_NEG;
+ 
+   u = MP_COPY(a);
+   v = MP_COPY(b);
+ 
+   while (MP_LEN(v)) {
+     mp *t;
+     gf_div(&q, &u, u, v);
+     if (f & f_ext) {
+       t = gf_mul(MP_NEW, X, q);
+       t = gf_add(t, x, t);
+       MP_DROP(x); x = X; X = t;
+       t = gf_mul(MP_NEW, Y, q);
+       t = gf_add(t, y, t);
+       MP_DROP(y); y = Y; Y = t;
+     }
+     t = u; u = v; v = t;
+   }
+ 
+   MP_DROP(q);
+   if (!gcd)
+     MP_DROP(u);
+   else {
+     if (*gcd) MP_DROP(*gcd);
+     u->f &= ~MP_NEG;
+     *gcd = u;
+   }
+ 
+   /* --- Perform a little normalization --- */
+ 
+   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;
+     }
+ 
+     /* --- 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(v);
+   MP_DROP(X); MP_DROP(Y);
+   MP_DROP(a); MP_DROP(b);
+ }
+ 
+ /*----- Test rig ----------------------------------------------------------*/
+ 
+ #ifdef TEST_RIG
+ 
+ static int gcd(dstr *v)
+ {
+   int ok = 1;
+   mp *a = *(mp **)v[0].buf;
+   mp *b = *(mp **)v[1].buf;
+   mp *g = *(mp **)v[2].buf;
+   mp *x = *(mp **)v[3].buf;
+   mp *y = *(mp **)v[4].buf;
+ 
+   mp *gg = MP_NEW, *xx = MP_NEW, *yy = MP_NEW;
+   gf_gcd(&gg, &xx, &yy, a, b);
+   if (!MP_EQ(x, xx)) {
+     fputs("\n*** mp_gcd(x) failed", stderr);
+     fputs("\na      = ", stderr); mp_writefile(a, stderr, 16);
+     fputs("\nb      = ", stderr); mp_writefile(b, stderr, 16);
+     fputs("\nexpect = ", stderr); mp_writefile(x, stderr, 16);
+     fputs("\nresult = ", stderr); mp_writefile(xx, stderr, 16);
+     fputc('\n', stderr);
+     ok = 0;
+   }
+   if (!MP_EQ(y, yy)) {
+     fputs("\n*** mp_gcd(y) failed", stderr);
+     fputs("\na      = ", stderr); mp_writefile(a, stderr, 16);
+     fputs("\nb      = ", stderr); mp_writefile(b, stderr, 16);
+     fputs("\nexpect = ", stderr); mp_writefile(y, stderr, 16);
+     fputs("\nresult = ", stderr); mp_writefile(yy, stderr, 16);
+     fputc('\n', stderr);
+     ok = 0;
+   }
+ 
+   if (!ok) {
+     mp *ax = gf_mul(MP_NEW, a, xx);
+     mp *by = gf_mul(MP_NEW, b, yy);
+     ax = gf_add(ax, ax, by);
+     if (MP_EQ(ax, gg))
+       fputs("\n*** (Alternative result found.)\n", stderr);
+     MP_DROP(ax);
+     MP_DROP(by);
+   }
+ 
+   if (!MP_EQ(g, gg)) {
+     fputs("\n*** mp_gcd(gcd) failed", stderr);
+     fputs("\na      = ", stderr); mp_writefile(a, stderr, 16);
+     fputs("\nb      = ", stderr); mp_writefile(b, stderr, 16);
+     fputs("\nexpect = ", stderr); mp_writefile(g, stderr, 16);
+     fputs("\nresult = ", stderr); mp_writefile(gg, stderr, 16);
+     fputc('\n', stderr);
+     ok = 0;
+   }
+   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);
+ }
+ 
+ static test_chunk tests[] = {
+   { "gcd", gcd, { &type_mp, &type_mp, &type_mp, &type_mp, &type_mp, 0 } },
+   { 0, 0, { 0 } }
+ };
+ 
+ int main(int argc, char *argv[])
+ {
+   sub_init();
+   test_run(argc, argv, tests, SRCDIR "/tests/gf");
+   return (0);
+ }
+ 
+ #endif
+ 
+ /*----- That's all, folks -------------------------------------------------*/
diff --cc gf.h
index 0000000,889cd9b..ebf67f1
mode 000000,100644..100644
--- a/gf.h
+++ b/gf.h
@@@ -1,0 -1,124 +1,127 @@@
+ /* -*-c-*-
+  *
 - * $Id: gf.h,v 1.1.2.1 2004/03/21 22:39:46 mdw Exp $
++ * $Id: gf.h,v 1.2 2004/03/21 22:52:06 mdw Exp $
+  *
+  * Arithmetic on binary polynomials
+  *
+  * (c) 2004 Straylight/Edgeware
+  */
+ 
+ /*----- Licensing notice --------------------------------------------------* 
+  *
+  * This file is part of Catacomb.
+  *
+  * Catacomb is free software; you can redistribute it and/or modify
+  * 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: gf.h,v $
++ * Revision 1.2  2004/03/21 22:52:06  mdw
++ * Merge and close elliptic curve branch.
++ *
+  * Revision 1.1.2.1  2004/03/21 22:39:46  mdw
+  * Elliptic curves on binary fields work.
+  *
+  */
+ 
+ #ifndef CATACOMB_GF_H
+ #define CATACOMB_GF_H
+ 
+ #ifdef __cplusplus
+   extern "C" {
+ #endif
+ 
+ /*----- Header files ------------------------------------------------------*/
+ 
+ #ifndef CATACOMB_MP_H
+ #  include "mp.h"
+ #endif
+ 
+ #ifndef CATACOMB_GFX_H
+ #  include "gfx.h"
+ #endif
+ 
+ /*----- Functions provided ------------------------------------------------*/
+ 
+ /* --- @gf_add@ --- *
+  *
+  * Arguments:	@mp *d@ = destination
+  *		@mp *a, *b@ = sources
+  *
+  * Returns:	Result, @a@ added to @b@.
+  */
+ 
+ extern mp *gf_add(mp */*d*/, mp */*a*/, mp */*b*/);
+ #define gf_sub gf_add
+ 
+ /* --- @gf_mul@ --- *
+  *
+  * Arguments:	@mp *d@ = destination
+  *		@mp *a, *b@ = sources
+  *
+  * Returns:	Result, @a@ multiplied by @b@.
+  */
+ 
+ extern mp *gf_mul(mp */*d*/, mp */*a*/, mp */*b*/);
+ 
+ /* --- @gf_sqr@ --- *
+  *
+  * Arguments:	@mp *d@ = destination
+  *		@mp *a@ = source
+  *
+  * Returns:	Result, @a@ squared.
+  */
+ 
+ extern mp *gf_sqr(mp */*d*/, mp */*a*/);
+ 
+ /* --- @gf_div@ --- *
+  *
+  * Arguments:	@mp **qq, **rr@ = destination, quotient and remainder
+  *		@mp *a, *b@ = sources
+  *
+  * Use:		Calculates the quotient and remainder when @a@ is divided by
+  *		@b@.  The destinations @*qq@ and @*rr@ must be distinct.
+  *		Either of @qq@ or @rr@ may be null to indicate that the
+  *		result is irrelevant.  (Discarding both results is silly.)
+  *		There is a performance advantage if @a == *rr@.
+  */
+ 
+ extern void gf_div(mp **/*qq*/, mp **/*rr*/, mp */*a*/, mp */*b*/);
+ 
+ /* --- @gf_gcd@ --- *
+  *
+  * Arguments:	@mp **gcd, **xx, **yy@ = where to write the results
+  *		@mp *a, *b@ = sources (must be nonzero)
+  *
+  *
+  * Returns:	---
+  *
+  * 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.
+  */
+ 
+ extern void gf_gcd(mp **/*gcd*/, mp **/*xx*/, mp **/*yy*/,
+ 		   mp */*a*/, mp */*b*/);
+ 
+ /*----- That's all, folks -------------------------------------------------*/
+ 
+ #ifdef __cplusplus
+   }
+ #endif
+ 
+ #endif
diff --cc gfreduce-exp.h
index 0000000,0393145..f826fc7
mode 000000,100644..100644
--- a/gfreduce-exp.h
+++ b/gfreduce-exp.h
@@@ -1,0 -1,84 +1,87 @@@
+ /* -*-c-*-
+  *
 - * $Id: gfreduce-exp.h,v 1.1.2.1 2004/03/21 22:39:46 mdw Exp $
++ * $Id: gfreduce-exp.h,v 1.2 2004/03/21 22:52:06 mdw Exp $
+  *
+  * Exponentiation operations for binary field reduction
+  *
+  * (c) 2004 Straylight/Edgeware
+  */
+ 
+ /*----- Licensing notice --------------------------------------------------* 
+  *
+  * This file is part of Catacomb.
+  *
+  * Catacomb is free software; you can redistribute it and/or modify
+  * 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: gfreduce-exp.h,v $
++ * Revision 1.2  2004/03/21 22:52:06  mdw
++ * Merge and close elliptic curve branch.
++ *
+  * Revision 1.1.2.1  2004/03/21 22:39:46  mdw
+  * Elliptic curves on binary fields work.
+  *
+  */
+ 
+ #ifndef CATACOMB_GFREDUCE_EXP_H
+ #define CATACOMB_GFREDUCE_EXP_H
+ 
+ #ifdef __cplusplus
+   extern "C" {
+ #endif
+ 
+ /*----- Exponentiation definitions ----------------------------------------*/
+ 
+ #define EXP_TYPE mp *
+ 
+ #define EXP_COPY(d, x) d = MP_COPY(x)
+ #define EXP_DROP(x) MP_DROP(x)
+ 
+ #define EXP_MUL(a, x) do {						\
+   mp *t = gf_mul(spare, a, x);						\
+   spare = a;								\
+   a = gfreduce_do(gr, t, t);						\
+ } while (0)
+ 
+ #define EXP_SQR(a) do {							\
+   mp *t = gf_sqr(spare, a);						\
+   spare = a;								\
+   a = gfreduce_do(gr, t, t);						\
+ } while (0)
+ 
+ #define EXP_FIX(x)
+ 
+ #define EXP_SETMUL(d, x, y) do {					\
+   d = gf_mul(MP_NEW, x, y);						\
+   d = gfreduce_do(gr, d, d);						\
+ } while (0)
+ 
+ #define EXP_SETSQR(d, x) do {						\
+   d = gf_sqr(MP_NEW, x);						\
+   d = gfreduce_do(gr, d, d);						\
+ } while (0)
+ 
+ #include "exp.h"
+ 
+ /*----- That's all, folks -------------------------------------------------*/
+ 
+ #ifdef __cplusplus
+   }
+ #endif
+ 
+ #endif
diff --cc gfreduce.c
index 0000000,3969f11..819c276
mode 000000,100644..100644
--- a/gfreduce.c
+++ b/gfreduce.c
@@@ -1,0 -1,652 +1,655 @@@
+ /* -*-c-*-
+  *
 - * $Id: gfreduce.c,v 1.1.2.1 2004/03/21 22:39:46 mdw Exp $
++ * $Id: gfreduce.c,v 1.2 2004/03/21 22:52:06 mdw Exp $
+  *
+  * Efficient reduction modulo sparse binary polynomials
+  *
+  * (c) 2004 Straylight/Edgeware
+  */
+ 
+ /*----- Licensing notice --------------------------------------------------* 
+  *
+  * This file is part of Catacomb.
+  *
+  * Catacomb is free software; you can redistribute it and/or modify
+  * 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: gfreduce.c,v $
++ * Revision 1.2  2004/03/21 22:52:06  mdw
++ * Merge and close elliptic curve branch.
++ *
+  * Revision 1.1.2.1  2004/03/21 22:39:46  mdw
+  * Elliptic curves on binary fields work.
+  *
+  */
+ 
+ /*----- Header files ------------------------------------------------------*/
+ 
+ #include <mLib/alloc.h>
+ #include <mLib/darray.h>
+ #include <mLib/macros.h>
+ 
+ #include "gf.h"
+ #include "gfreduce.h"
+ #include "gfreduce-exp.h"
+ #include "fibrand.h"
+ #include "mprand.h"
+ 
+ /*----- Data structures ---------------------------------------------------*/
+ 
+ DA_DECL(instr_v, gfreduce_instr);
+ 
+ /*----- Main code ---------------------------------------------------------*/
+ 
+ /* --- What's going on here? --- *
+  *
+  * Let's face it, @gfx_div@ sucks.  It works (I hope), but it's not in any
+  * sense fast.  Here, we do efficient reduction modulo sparse polynomials.
+  *
+  * Suppose we have a polynomial @X@ we're trying to reduce mod @P@.  If we
+  * take the topmost nonzero word of @X@, call it @w@, then we can eliminate
+  * it by subtracting off @w P x^{k}@ for an appropriate value of @k@.  The
+  * trick is in observing that if @P@ is sparse we can do this multiplication
+  * and subtraction efficiently, just by XORing appropriate shifts of @w@ into
+  * @X@.
+  *
+  * The first tricky bit is in working out when to stop.  I'll use eight-bit
+  * words to demonstrate what I'm talking about.
+  *
+  *  xxxxxx xxxxxxxx xxxxxxxx xxxxxxxx xxxxxxxx xxxxxxxx
+  *                  001ppppp pppppppp pppppppp pppppppp
+  *                    |<rp>|
+  *                    |<------------ bp ------------->|
+  *                  |<------------ nw --------------->|
+  *
+  * The trick of taking whole words off of @X@ stops working when there are
+  * only @nw@ words left.  Then we have to mask off the bottom bits of @w@
+  * before continuing.
+  */
+ 
+ /* --- @gfreduce_create@ --- *
+  *
+  * Arguments:	@gfreduce *r@ = structure to fill in
+  *		@mp *x@ = a (hopefully sparse) polynomial
+  *
+  * Returns:	---
+  *
+  * Use:		Initializes a context structure for reduction.
+  */
+ 
+ void gfreduce_create(gfreduce *r, mp *p)
+ {
+   instr_v iv = DA_INIT;
+   unsigned long d, dw;
+   mpscan sc;
+   unsigned long i;
+   gfreduce_instr *ip;
+   unsigned f = 0;
+   size_t w, ww, wi, wl, ll;
+ 
+   /* --- Sort out the easy stuff --- */
+ 
+   d = mp_bits(p); assert(d); d--;
+   r->lim = d/MPW_BITS;
+   dw = d%MPW_BITS;
+   if (!dw)
+     r->mask = 0;
+   else {
+     r->mask = MPW(((mpw)-1) << dw);
+     r->lim++;
+   }
+   r->p = mp_copy(p);
+ 
+   /* --- Stash a new instruction --- */
+ 
+ #define INSTR(op_, arg_) do {						\
+   DA_ENSURE(&iv, 1);							\
+   DA(&iv)[DA_LEN(&iv)].op = (op_);					\
+   DA(&iv)[DA_LEN(&iv)].arg = (arg_);					\
+   DA_EXTEND(&iv, 1);							\
+ } while (0)
+ 
+ #define f_lsr 1u
+ 
+   w = (d + MPW_BITS - 1)/MPW_BITS;
+   INSTR(GFRI_LOAD, w);
+   wi = DA_LEN(&iv);
+   f = 0;
+   ll = 0;
+   for (i = 0, mp_scan(&sc, p); mp_step(&sc) && i < d; i++) {
+     if (!mp_bit(&sc))
+       continue;
+     ww = (d - i + MPW_BITS - 1)/MPW_BITS;
+     if (ww != w) {
+       wl = DA_LEN(&iv);
+       INSTR(GFRI_STORE, w);
+       if (!ll)
+ 	ll = DA_LEN(&iv);
+       if (!(f & f_lsr))
+ 	INSTR(GFRI_LOAD, ww);
+       else {
+ 	INSTR(GFRI_LOAD, w - 1);
+ 	for (; wi < wl; wi++) {
+ 	  ip = &DA(&iv)[wi];
+ 	  assert(ip->op == GFRI_LSL);
+ 	  if (ip->arg)
+ 	    INSTR(GFRI_LSR, MPW_BITS - ip->arg);
+ 	}
+ 	if (w - 1 != ww) {
+ 	  INSTR(GFRI_STORE, w - 1);
+ 	  INSTR(GFRI_LOAD, ww);
+ 	}
+ 	f &= ~f_lsr;
+       }
+       w = ww;
+       wi = DA_LEN(&iv);
+     }
+     INSTR(GFRI_LSL, (i - d)%MPW_BITS);
+     if ((i - d)%MPW_BITS)
+       f |= f_lsr;
+   }
+   wl = DA_LEN(&iv);
+   INSTR(GFRI_STORE, w);
+   if (!ll)
+     ll = DA_LEN(&iv);
+   if (f & f_lsr) {
+     INSTR(GFRI_LOAD, w - 1);
+     for (; wi < wl; wi++) {
+       ip = &DA(&iv)[wi];
+       assert(ip->op == GFRI_LSL);
+       if (ip->arg)
+ 	INSTR(GFRI_LSR, MPW_BITS - ip->arg);
+     }
+     INSTR(GFRI_STORE, w - 1);
+   }
+ 
+ #undef INSTR
+ 
+   r->in = DA_LEN(&iv);
+   r->iv = xmalloc(r->in * sizeof(gfreduce_instr));
+   r->liv = r->iv + ll;
+   memcpy(r->iv, DA(&iv), r->in * sizeof(gfreduce_instr));
+   DA_DESTROY(&iv);
+ }
+ 
+ /* --- @gfreduce_destroy@ --- *
+  *
+  * Arguments:	@gfreduce *r@ = structure to free
+  *
+  * Returns:	---
+  *
+  * Use:		Reclaims the resources from a reduction context.
+  */
+ 
+ void gfreduce_destroy(gfreduce *r)
+ {
+   mp_drop(r->p);
+   xfree(r->iv);
+ }
+ 
+ /* --- @gfreduce_dump@ --- *
+  *
+  * Arguments:	@gfreduce *r@ = structure to dump
+  *		@FILE *fp@ = file to dump on
+  *
+  * Returns:	---
+  *
+  * Use:		Dumps a reduction context.
+  */
+ 
+ void gfreduce_dump(gfreduce *r, FILE *fp)
+ {
+   size_t i;
+ 
+   fprintf(fp, "poly = "); mp_writefile(r->p, fp, 16);
+   fprintf(fp, "\n  lim = %lu; mask = %lx\n",
+ 	  (unsigned long)r->lim, (unsigned long)r->mask);
+   for (i = 0; i < r->in; i++) {
+     static const char *opname[] = { "load", "lsl", "lsr", "store" };
+     assert(r->iv[i].op < N(opname));
+     fprintf(fp, "  %s %lu\n",
+ 	    opname[r->iv[i].op],
+ 	    (unsigned long)r->iv[i].arg);
+   }
+ }
+ 
+ /* --- @gfreduce_do@ --- *
+  *
+  * Arguments:	@gfreduce *r@ = reduction context
+  *		@mp *d@ = destination
+  *		@mp *x@ = source
+  *
+  * Returns:	Destination, @x@ reduced modulo the reduction poly.
+  */
+ 
+ static void run(const gfreduce_instr *i, const gfreduce_instr *il,
+ 		mpw *v, mpw z)
+ {
+   mpw w = 0;
+ 
+   for (; i < il; i++) {
+     switch (i->op) {
+       case GFRI_LOAD: w = *(v - i->arg); break;
+       case GFRI_LSL: w ^= z << i->arg; break;
+       case GFRI_LSR: w ^= z >> i->arg; break;
+       case GFRI_STORE: *(v - i->arg) = MPW(w); break;
+       default: abort();
+     }
+   }
+ }
+ 
+ mp *gfreduce_do(gfreduce *r, mp *d, mp *x)
+ {
+   mpw *v, *vl;
+   const gfreduce_instr *il;
+   mpw z;
+ 
+   /* --- Try to reuse the source's space --- */
+ 
+   MP_COPY(x);
+   if (d) MP_DROP(d);
+   MP_DEST(x, MP_LEN(x), x->f);
+ 
+   /* --- Do the reduction --- */
+ 
+   il = r->iv + r->in;
+   if (MP_LEN(x) >= r->lim) {
+     v = x->v + r->lim;
+     vl = x->vl;
+     while (vl-- > v) {
+       while (*vl) {
+ 	z = *vl;
+ 	*vl = 0;
+ 	run(r->iv, il, vl, z);
+       }
+     }
+     if (r->mask) {
+       while (*vl & r->mask) {
+ 	z = *vl & r->mask;
+ 	*vl &= ~r->mask;
+ 	run(r->iv, il, vl, z);
+       }
+     }
+   }
+ 
+   /* --- Done --- */
+ 
+   MP_SHRINK(x);
+   return (x);
+ }
+ 
+ /* --- @gfreduce_sqrt@ --- *
+  *
+  * Arguments:	@gfreduce *r@ = pointer to reduction context
+  *		@mp *d@ = destination
+  *		@mp *x@ = some polynomial
+  *
+  * Returns:	The square root of @x@ modulo @r->p@, or null.
+  */
+ 
+ mp *gfreduce_sqrt(gfreduce *r, mp *d, mp *x)
+ {
+   mp *y = MP_COPY(x);
+   mp *z, *spare = MP_NEW;
+   unsigned long m = mp_bits(r->p) - 1;
+   unsigned long i;
+ 
+   for (i = 0; i < m - 1; i++) {
+     mp *t = gf_sqr(spare, y);
+     spare = y;
+     y = gfreduce_do(r, t, t);
+   }
+   z = gf_sqr(spare, y);
+   z = gfreduce_do(r, z, z);
+   if (!MP_EQ(x, z)) {
+     mp_drop(y);
+     y = 0;
+   }
+   mp_drop(z);
+   mp_drop(d);
+   return (y);
+ }
+ 
+ /* --- @gfreduce_trace@ --- *
+  *
+  * Arguments:	@gfreduce *r@ = pointer to reduction context
+  *		@mp *x@ = some polynomial
+  *
+  * Returns:	The trace of @x@. (%$\Tr(x)=x + x^2 + \cdots + x^{2^{m-1}}$%
+  *		if %$x \in \gf{2^m}$%).
+  */
+ 
+ int gfreduce_trace(gfreduce *r, mp *x)
+ {
+   mp *y = MP_COPY(x);
+   mp *spare = MP_NEW;
+   unsigned long m = mp_bits(r->p) - 1;
+   unsigned long i;
+   int rc;
+ 
+   for (i = 0; i < m - 1; i++) {
+     mp *t = gf_sqr(spare, y);
+     spare = y;
+     y = gfreduce_do(r, t, t);
+     y = gf_add(y, y, x);
+   }
+   rc = !MP_ISZERO(y);
+   mp_drop(spare);
+   mp_drop(y);
+   return (rc);
+ }
+ 
+ /* --- @gfreduce_halftrace@ --- *
+  *
+  * Arguments:	@gfreduce *r@ = pointer to reduction context
+  *		@mp *d@ = destination
+  *		@mp *x@ = some polynomial
+  *
+  * Returns:	The half-trace of @x@.
+  *		(%$\HfTr(x)= x + x^{2^2} + \cdots + x^{2^{m-1}}$%
+  *		if %$x \in \gf{2^m}$% with %$m$% odd).
+  */
+ 
+ mp *gfreduce_halftrace(gfreduce *r, mp *d, mp *x)
+ {
+   mp *y = MP_COPY(x);
+   mp *spare = MP_NEW;
+   unsigned long m = mp_bits(r->p) - 1;
+   unsigned long i;
+ 
+   mp_drop(d);
+   for (i = 0; i < m - 1; i += 2) {
+     mp *t = gf_sqr(spare, y);
+     spare = y;
+     y = gfreduce_do(r, t, t);
+     t = gf_sqr(spare, y);
+     spare = y;
+     y = gfreduce_do(r, t, t);
+     y = gf_add(y, y, x);
+   }
+   mp_drop(spare);
+   return (y);
+ }
+ 
+ /* --- @gfreduce_quadsolve@ --- *
+  *
+  * Arguments:	@gfreduce *r@ = pointer to reduction context
+  *		@mp *d@ = destination
+  *		@mp *x@ = some polynomial
+  *
+  * Returns:	A polynomial @y@ such that %$y^2 + y = x$%, or null.
+  */
+ 
+ mp *gfreduce_quadsolve(gfreduce *r, mp *d, mp *x)
+ {
+   unsigned long m = mp_bits(r->p) - 1;
+   mp *t;
+ 
+   MP_COPY(x);
+   if (m & 1)
+     d = gfreduce_halftrace(r, d, x);
+   else {
+     mp *z, *w, *rho = MP_NEW;
+     mp *spare = MP_NEW;
+     grand *fr = fibrand_create(0);
+     unsigned long i;
+ 
+     for (;;) {
+       rho = mprand(rho, m, fr, 0);
+       z = MP_ZERO;
+       w = MP_COPY(rho);
+       for (i = 0; i < m - 1; i++) {
+ 	t = gf_sqr(spare, z); spare = z; z = gfreduce_do(r, t, t);
+ 	t = gf_sqr(spare, w); spare = w; w = gfreduce_do(r, t, t);
+ 	t = gf_mul(spare, w, x); t = gfreduce_do(r, t, t); spare = t;
+ 	z = gf_add(z, z, t);
+ 	w = gf_add(w, w, rho);
+       }
+       if (!MP_ISZERO(w))
+ 	break;
+       MP_DROP(z);
+       MP_DROP(w);
+     }
+     if (d) MP_DROP(d);
+     MP_DROP(w);
+     MP_DROP(spare);
+     MP_DROP(rho);
+     fr->ops->destroy(fr);
+     d = z;
+   }
+ 
+   t = gf_sqr(MP_NEW, d); t = gfreduce_do(r, t, t); t = gf_add(t, t, d);
+   if (!MP_EQ(t, x)) {
+     MP_DROP(d);
+     d = 0;
+   }
+   MP_DROP(t);
+   MP_DROP(x);
+   d->v[0] &= ~(mpw)1;
+   return (d);
+ }
+ 
+ /* --- @gfreduce_exp@ --- *
+  *
+  * Arguments:	@gfreduce *gr@ = pointer to reduction context
+  *              @mp *d@ = fake destination
+  *              @mp *a@ = base
+  *              @mp *e@ = exponent
+  *
+  * Returns:     Result, %$a^e \bmod m$%.
+  */
+ 
+ mp *gfreduce_exp(gfreduce *gr, mp *d, mp *a, mp *e)
+ {
+   mp *x = MP_ONE;
+   mp *spare = (e->f & MP_BURN) ? MP_NEWSEC : MP_NEW;
+ 
+   MP_SHRINK(e);
+   if (!MP_LEN(e))
+     ;
+   else if (MP_LEN(e) < EXP_THRESH)
+     EXP_SIMPLE(x, a, e);
+   else
+     EXP_WINDOW(x, a, e);
+   mp_drop(d);
+   mp_drop(spare);
+   return (x);
+ }
+ 
+ /*----- Test rig ----------------------------------------------------------*/
+ 
+ #ifdef TEST_RIG
+ 
+ #define MP(x) mp_readstring(MP_NEW, #x, 0, 0)
+ 
+ static int vreduce(dstr *v)
+ {
+   mp *d = *(mp **)v[0].buf;
+   mp *n = *(mp **)v[1].buf;
+   mp *r = *(mp **)v[2].buf;
+   mp *c;
+   int ok = 1;
+   gfreduce rr;
+ 
+   gfreduce_create(&rr, d);
+   c = gfreduce_do(&rr, MP_NEW, n);
+   if (!MP_EQ(c, r)) {
+     fprintf(stderr, "\n*** reduction failed\n*** ");
+     gfreduce_dump(&rr, stderr);
+     fprintf(stderr, "\n*** n = "); mp_writefile(n, stderr, 16);
+     fprintf(stderr, "\n*** r = "); mp_writefile(r, stderr, 16);
+     fprintf(stderr, "\n*** c = "); mp_writefile(c, stderr, 16);
+     fprintf(stderr, "\n");
+     ok = 0;
+   }
+   gfreduce_destroy(&rr);
+   mp_drop(n); mp_drop(d); mp_drop(r); mp_drop(c);
+   assert(mparena_count(MPARENA_GLOBAL) == 0);
+   return (ok);
+ }
+ 
+ static int vmodexp(dstr *v)
+ {
+   mp *p = *(mp **)v[0].buf;
+   mp *g = *(mp **)v[1].buf;
+   mp *x = *(mp **)v[2].buf;
+   mp *r = *(mp **)v[3].buf;
+   mp *c;
+   int ok = 1;
+   gfreduce rr;
+ 
+   gfreduce_create(&rr, p);
+   c = gfreduce_exp(&rr, MP_NEW, g, x);
+   if (!MP_EQ(c, r)) {
+     fprintf(stderr, "\n*** modexp failed\n*** ");
+     fprintf(stderr, "\n*** p = "); mp_writefile(p, stderr, 16);
+     fprintf(stderr, "\n*** g = "); mp_writefile(g, stderr, 16);
+     fprintf(stderr, "\n*** x = "); mp_writefile(x, stderr, 16);
+     fprintf(stderr, "\n*** c = "); mp_writefile(c, stderr, 16);
+     fprintf(stderr, "\n*** r = "); mp_writefile(r, stderr, 16);
+     fprintf(stderr, "\n");
+     ok = 0;
+   }
+   gfreduce_destroy(&rr);
+   mp_drop(p); mp_drop(g); mp_drop(r); mp_drop(x); mp_drop(c);
+   assert(mparena_count(MPARENA_GLOBAL) == 0);
+   return (ok);
+ }
+ 
+ static int vsqrt(dstr *v)
+ {
+   mp *p = *(mp **)v[0].buf;
+   mp *x = *(mp **)v[1].buf;
+   mp *r = *(mp **)v[2].buf;
+   mp *c;
+   int ok = 1;
+   gfreduce rr;
+ 
+   gfreduce_create(&rr, p);
+   c = gfreduce_sqrt(&rr, MP_NEW, x);
+   if (!MP_EQ(c, r)) {
+     fprintf(stderr, "\n*** sqrt failed\n*** ");
+     fprintf(stderr, "\n*** p = "); mp_writefile(p, stderr, 16);
+     fprintf(stderr, "\n*** x = "); mp_writefile(x, stderr, 16);
+     fprintf(stderr, "\n*** c = "); mp_writefile(c, stderr, 16);
+     fprintf(stderr, "\n*** r = "); mp_writefile(r, stderr, 16);
+     fprintf(stderr, "\n");
+     ok = 0;
+   }
+   gfreduce_destroy(&rr);
+   mp_drop(p); mp_drop(r); mp_drop(x); mp_drop(c);
+   assert(mparena_count(MPARENA_GLOBAL) == 0);
+   return (ok);
+ }
+ 
+ static int vtr(dstr *v)
+ {
+   mp *p = *(mp **)v[0].buf;
+   mp *x = *(mp **)v[1].buf;
+   int r = *(int *)v[2].buf, c;
+   int ok = 1;
+   gfreduce rr;
+ 
+   gfreduce_create(&rr, p);
+   c = gfreduce_trace(&rr, x);
+   if (c != r) {
+     fprintf(stderr, "\n*** trace failed\n*** ");
+     fprintf(stderr, "\n*** p = "); mp_writefile(p, stderr, 16);
+     fprintf(stderr, "\n*** x = "); mp_writefile(x, stderr, 16);
+     fprintf(stderr, "\n*** c = %d", c);
+     fprintf(stderr, "\n*** r = %d", r);
+     fprintf(stderr, "\n");
+     ok = 0;
+   }
+   gfreduce_destroy(&rr);
+   mp_drop(p); mp_drop(x); 
+   assert(mparena_count(MPARENA_GLOBAL) == 0);
+   return (ok);
+ }
+ 
+ static int vhftr(dstr *v)
+ {
+   mp *p = *(mp **)v[0].buf;
+   mp *x = *(mp **)v[1].buf;
+   mp *r = *(mp **)v[2].buf;
+   mp *c;
+   int ok = 1;
+   gfreduce rr;
+ 
+   gfreduce_create(&rr, p);
+   c = gfreduce_halftrace(&rr, MP_NEW, x);
+   if (!MP_EQ(c, r)) {
+     fprintf(stderr, "\n*** halftrace failed\n*** ");
+     fprintf(stderr, "\n*** p = "); mp_writefile(p, stderr, 16);
+     fprintf(stderr, "\n*** x = "); mp_writefile(x, stderr, 16);
+     fprintf(stderr, "\n*** c = "); mp_writefile(c, stderr, 16);
+     fprintf(stderr, "\n*** r = "); mp_writefile(r, stderr, 16);
+     fprintf(stderr, "\n");
+     ok = 0;
+   }
+   gfreduce_destroy(&rr);
+   mp_drop(p); mp_drop(r); mp_drop(x); mp_drop(c);
+   assert(mparena_count(MPARENA_GLOBAL) == 0);
+   return (ok);
+ }
+ 
+ static int vquad(dstr *v)
+ {
+   mp *p = *(mp **)v[0].buf;
+   mp *x = *(mp **)v[1].buf;
+   mp *r = *(mp **)v[2].buf;
+   mp *c;
+   int ok = 1;
+   gfreduce rr;
+ 
+   gfreduce_create(&rr, p);
+   c = gfreduce_quadsolve(&rr, MP_NEW, x);
+   if (!MP_EQ(c, r)) {
+     fprintf(stderr, "\n*** quadsolve failed\n*** ");
+     fprintf(stderr, "\n*** p = "); mp_writefile(p, stderr, 16);
+     fprintf(stderr, "\n*** x = "); mp_writefile(x, stderr, 16);
+     fprintf(stderr, "\n*** c = "); mp_writefile(c, stderr, 16);
+     fprintf(stderr, "\n*** r = "); mp_writefile(r, stderr, 16);
+     fprintf(stderr, "\n");
+     ok = 0;
+   }
+   gfreduce_destroy(&rr);
+   mp_drop(p); mp_drop(r); mp_drop(x); mp_drop(c);
+   assert(mparena_count(MPARENA_GLOBAL) == 0);
+   return (ok);
+ }
+ 
+ static test_chunk defs[] = {
+   { "reduce", vreduce, { &type_mp, &type_mp, &type_mp, 0 } },
+   { "modexp", vmodexp, { &type_mp, &type_mp, &type_mp, &type_mp, 0 } },
+   { "sqrt", vsqrt, { &type_mp, &type_mp, &type_mp, 0 } },
+   { "trace", vtr, { &type_mp, &type_mp, &type_int, 0 } },
+   { "halftrace", vhftr, { &type_mp, &type_mp, &type_mp, 0 } },
+   { "quadsolve", vquad, { &type_mp, &type_mp, &type_mp, 0 } },
+   { 0, 0, { 0 } }
+ };
+ 
+ int main(int argc, char *argv[])
+ {
+   test_run(argc, argv, defs, SRCDIR"/tests/gfreduce");
+   return (0);
+ }
+ 
+ #endif
+ 
+ /*----- That's all, folks -------------------------------------------------*/
diff --cc gfreduce.h
index 0000000,2fc4c0a..9840b5e
mode 000000,100644..100644
--- a/gfreduce.h
+++ b/gfreduce.h
@@@ -1,0 -1,186 +1,189 @@@
+ /* -*-c-*-
+  *
 - * $Id: gfreduce.h,v 1.1.2.1 2004/03/21 22:39:46 mdw Exp $
++ * $Id: gfreduce.h,v 1.2 2004/03/21 22:52:06 mdw Exp $
+  *
+  * Reduction modulo sparse binary polynomials
+  *
+  * (c) 2004 Straylight/Edgeware
+  */
+ 
+ /*----- Licensing notice --------------------------------------------------* 
+  *
+  * This file is part of Catacomb.
+  *
+  * Catacomb is free software; you can redistribute it and/or modify
+  * 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: gfreduce.h,v $
++ * Revision 1.2  2004/03/21 22:52:06  mdw
++ * Merge and close elliptic curve branch.
++ *
+  * Revision 1.1.2.1  2004/03/21 22:39:46  mdw
+  * Elliptic curves on binary fields work.
+  *
+  */
+ 
+ #ifndef CATACOMB_GFREDUCE_H
+ #define CATACOMB_GFREDUCE_H
+ 
+ #ifdef __cplusplus
+   extern "C" {
+ #endif
+ 
+ /*----- Header files ------------------------------------------------------*/
+ 
+ #ifndef CATACOMB_GF_H
+ #  include "gf.h"
+ #endif
+ 
+ /*----- Data structures ---------------------------------------------------*/
+ 
+ typedef struct gfreduce_instr {
+   unsigned op;				/* Instruction opcode */
+   size_t arg;				/* Immediate argument */
+ } gfreduce_instr;
+ 
+ enum {
+   GFRI_LOAD,				/* Load @p[arg]@ */
+   GFRI_LSL,				/* XOR with @w << arg@ */
+   GFRI_LSR,				/* XOR with @w >> arg@ */
+   GFRI_STORE,				/* Store @p[arg]@ */
+   GFRI_MAX
+ };
+ 
+ typedef struct gfreduce {
+   size_t lim;				/* Word of degree bit */
+   mpw mask;				/* Mask for degree word */
+   mp *p;				/* Copy of the polynomial */
+   size_t in;				/* Number of instruction words */
+   gfreduce_instr *iv, *liv;		/* Vector of instructions */
+ } gfreduce;
+ 
+ /*----- Functions provided ------------------------------------------------*/
+ 
+ /* --- @gfreduce_create@ --- *
+  *
+  * Arguments:	@gfreduce *r@ = structure to fill in
+  *		@mp *x@ = a (hopefully sparse) polynomial
+  *
+  * Returns:	---
+  *
+  * Use:		Initializes a context structure for reduction.
+  */
+ 
+ extern void gfreduce_create(gfreduce */*r*/, mp */*p*/);
+ 
+ /* --- @gfreduce_destroy@ --- *
+  *
+  * Arguments:	@gfreduce *r@ = structure to free
+  *
+  * Returns:	---
+  *
+  * Use:		Reclaims the resources from a reduction context.
+  */
+ 
+ extern void gfreduce_destroy(gfreduce */*r*/);
+ 
+ /* --- @gfreduce_dump@ --- *
+  *
+  * Arguments:	@gfreduce *r@ = structure to dump
+  *		@FILE *fp@ = file to dump on
+  *
+  * Returns:	---
+  *
+  * Use:		Dumps a reduction context.
+  */
+ 
+ extern void gfreduce_dump(gfreduce */*r*/, FILE */*fp*/);
+ 
+ /* --- @gfreduce_do@ --- *
+  *
+  * Arguments:	@gfreduce *r@ = reduction context
+  *		@mp *d@ = destination
+  *		@mp *x@ = source
+  *
+  * Returns:	Destination, @x@ reduced modulo the reduction poly.
+  */
+ 
+ extern mp *gfreduce_do(gfreduce */*r*/, mp */*d*/, mp */*x*/);
+ 
+ /* --- @gfreduce_sqrt@ --- *
+  *
+  * Arguments:	@gfreduce *r@ = pointer to reduction context
+  *		@mp *d@ = destination
+  *		@mp *x@ = some polynomial
+  *
+  * Returns:	The square root of @x@ modulo @r->p@, or null.
+  */
+ 
+ extern mp *gfreduce_sqrt(gfreduce */*r*/, mp */*d*/, mp */*x*/);
+ 
+ /* --- @gfreduce_trace@ --- *
+  *
+  * Arguments:	@gfreduce *r@ = pointer to reduction context
+  *		@mp *x@ = some polynomial
+  *
+  * Returns:	The trace of @x@. (%$\Tr(x)=x + x^2 + \cdots + x^{2^{m-1}}$%
+  *		if %$x \in \gf{2^m}$%).
+  */
+ 
+ extern int gfreduce_trace(gfreduce */*r*/, mp */*x*/);
+ 
+ /* --- @gfreduce_halftrace@ --- *
+  *
+  * Arguments:	@gfreduce *r@ = pointer to reduction context
+  *		@mp *d@ = destination
+  *		@mp *x@ = some polynomial
+  *
+  * Returns:	The half-trace of @x@.
+  *		(%$\HfTr(x)= x + x^{2^2} + \cdots + x^{2^{m-1}}$%
+  *		if %$x \in \gf{2^m}$% with %$m$% odd).
+  */
+ 
+ extern mp *gfreduce_halftrace(gfreduce */*r*/, mp */*d*/, mp */*x*/);
+ 
+ /* --- @gfreduce_quadsolve@ --- *
+  *
+  * Arguments:	@gfreduce *r@ = pointer to reduction context
+  *		@mp *d@ = destination
+  *		@mp *x@ = some polynomial
+  *
+  * Returns:	A polynomial @y@ such that %$y^2 + y = x$%, or null.
+  */
+ 
+ extern mp *gfreduce_quadsolve(gfreduce */*r*/, mp */*d*/, mp */*x*/);
+ 
+ /* --- @gfreduce_exp@ --- *
+  *
+  * Arguments:	@gfreduce *gr@ = pointer to reduction context
+  *              @mp *d@ = fake destination
+  *              @mp *a@ = base
+  *              @mp *e@ = exponent
+  *
+  * Returns:     Result, %$a^e \bmod m$%.
+  */
+ 
+ extern mp *gfreduce_exp(gfreduce */*gr*/, mp */*d*/, mp */*a*/, mp */*e*/);
+ 
+ /*----- That's all, folks -------------------------------------------------*/
+ 
+ #ifdef __cplusplus
+   }
+ #endif
+ 
+ #endif
diff --cc gfx-sqr.c
index b253a32,778f85a..19ec574
--- a/gfx-sqr.c
+++ b/gfx-sqr.c
@@@ -1,6 -1,6 +1,6 @@@
  /* -*-c-*-
   *
-  * $Id: gfx-sqr.c,v 1.1 2000/10/08 15:49:37 mdw Exp $
 - * $Id: gfx-sqr.c,v 1.1.4.1 2004/03/21 22:39:46 mdw Exp $
++ * $Id: gfx-sqr.c,v 1.2 2004/03/21 22:52:06 mdw Exp $
   *
   * Sqaring binary polynomials
   *
@@@ -30,6 -30,9 +30,12 @@@
  /*----- Revision history --------------------------------------------------* 
   *
   * $Log: gfx-sqr.c,v $
++ * Revision 1.2  2004/03/21 22:52:06  mdw
++ * Merge and close elliptic curve branch.
++ *
+  * Revision 1.1.4.1  2004/03/21 22:39:46  mdw
+  * Elliptic curves on binary fields work.
+  *
   * Revision 1.1  2000/10/08 15:49:37  mdw
   * First glimmerings of binary polynomial arithmetic.
   *
diff --cc gfx.h
index 778baee,d525650..18ac9a5
--- a/gfx.h
+++ b/gfx.h
@@@ -1,6 -1,6 +1,6 @@@
  /* -*-c-*-
   *
-  * $Id: gfx.h,v 1.1 2000/10/08 15:49:37 mdw Exp $
 - * $Id: gfx.h,v 1.1.4.1 2004/03/21 22:39:46 mdw Exp $
++ * $Id: gfx.h,v 1.2 2004/03/21 22:52:06 mdw Exp $
   *
   * Low-level arithmetic on binary polynomials
   *
@@@ -30,6 -30,9 +30,12 @@@
  /*----- Revision history --------------------------------------------------* 
   *
   * $Log: gfx.h,v $
++ * Revision 1.2  2004/03/21 22:52:06  mdw
++ * Merge and close elliptic curve branch.
++ *
+  * Revision 1.1.4.1  2004/03/21 22:39:46  mdw
+  * Elliptic curves on binary fields work.
+  *
   * Revision 1.1  2000/10/08 15:49:37  mdw
   * First glimmerings of binary polynomial arithmetic.
   *
diff --cc mp-gcd.c
index a17a502,f55d0aa..6135e54
--- a/mp-gcd.c
+++ b/mp-gcd.c
@@@ -1,6 -1,6 +1,6 @@@
  /* -*-c-*-
   *
-  * $Id: mp-gcd.c,v 1.5 2000/10/08 12:02:41 mdw Exp $
 - * $Id: mp-gcd.c,v 1.5.4.1 2004/03/21 22:39:46 mdw Exp $
++ * $Id: mp-gcd.c,v 1.6 2004/03/21 22:52:06 mdw Exp $
   *
   * Extended GCD calculation
   *
@@@ -30,6 -30,9 +30,12 @@@
  /*----- Revision history --------------------------------------------------* 
   *
   * $Log: mp-gcd.c,v $
++ * Revision 1.6  2004/03/21 22:52:06  mdw
++ * Merge and close elliptic curve branch.
++ *
+  * Revision 1.5.4.1  2004/03/21 22:39:46  mdw
+  * Elliptic curves on binary fields work.
+  *
   * Revision 1.5  2000/10/08 12:02:41  mdw
   * Use Euclid's algorithm rather than the binary one.
   *
diff --cc mpbarrett-exp.h
index 5c34746,e34ea49..dd02637
--- a/mpbarrett-exp.h
+++ b/mpbarrett-exp.h
@@@ -1,6 -1,6 +1,6 @@@
  /* -*-c-*-
   *
-  * $Id: mpbarrett-exp.h,v 1.1 2001/06/16 12:58:12 mdw Exp $
 - * $Id: mpbarrett-exp.h,v 1.1.4.1 2004/03/20 00:20:05 mdw Exp $
++ * $Id: mpbarrett-exp.h,v 1.2 2004/03/21 22:52:06 mdw Exp $
   *
   * Exponentiation operations for Barrett reduction
   *
@@@ -30,6 -30,9 +30,12 @@@
  /*----- Revision history --------------------------------------------------* 
   *
   * $Log: mpbarrett-exp.h,v $
++ * Revision 1.2  2004/03/21 22:52:06  mdw
++ * Merge and close elliptic curve branch.
++ *
+  * Revision 1.1.4.1  2004/03/20 00:20:05  mdw
+  * Projective coordinates for prime curves
+  *
   * Revision 1.1  2001/06/16 12:58:12  mdw
   * Parameters for generic exponentiation.
   *
diff --cc mpmont-exp.h
index e732e6f,0f82b9c..5f2b31d
--- a/mpmont-exp.h
+++ b/mpmont-exp.h
@@@ -1,6 -1,6 +1,6 @@@
  /* -*-c-*-
   *
-  * $Id: mpmont-exp.h,v 1.1 2001/06/16 12:58:12 mdw Exp $
 - * $Id: mpmont-exp.h,v 1.1.4.1 2004/03/20 00:13:31 mdw Exp $
++ * $Id: mpmont-exp.h,v 1.2 2004/03/21 22:52:06 mdw Exp $
   *
   * Exponentiation operations for Montgomery reduction
   *
@@@ -30,6 -30,9 +30,12 @@@
  /*----- Revision history --------------------------------------------------* 
   *
   * $Log: mpmont-exp.h,v $
++ * Revision 1.2  2004/03/21 22:52:06  mdw
++ * Merge and close elliptic curve branch.
++ *
+  * Revision 1.1.4.1  2004/03/20 00:13:31  mdw
+  * Projective coordinates for prime curves
+  *
   * Revision 1.1  2001/06/16 12:58:12  mdw
   * Parameters for generic exponentiation.
   *
diff --cc tests/gf
index 0000000,0c3987f..bbb1514
mode 000000,100644..100644
--- a/tests/gf
+++ b/tests/gf
@@@ -1,0 -1,70 +1,70 @@@
 -# $Id: gf,v 1.1.2.1 2004/03/21 22:39:46 mdw Exp $
++# $Id: gf,v 1.2 2004/03/21 22:52:06 mdw Exp $
+ #
+ # Test cases for higher-level binary poly arithmetic.
+ 
+ add {
+   0 0 0;
+   1 1 0;
+   1 2 3;
+   4 5 1;
+   0x7fb838a8a0a95046b9d9d9fb4440f7bb
+     0xc1a7bd3b4e853fc92d4e1588719986aa
+     0xbe1f8593ee2c6f8f9497cc7335d97111;
+   0x1e2933215e1c3bba8d2b404d98f43086bfc6198a219b168f214042a5e7df6b21
+     0x1e2933215e1c3bba8d2b404d98f43086bfc6198a219b168f214042a5e7df6b22 3;  
+ }
+ 
+ mul {
+   0 0 0;
+   1 0 0;
+   0 1 0;
+   1 1 1;
+   0x7fb838a8a0a95046b9d9d9fb4440f7bb
+     0xc1a7bd3b4e853fc92d4e1588719986aa
+     0x207ccad257b4ed64447158315bfb9aca5cbc5622cfb8fcbb1380eea1bc5c624e;
+   0xc1a7bd3b4e853fc92d4e1588719986aa
+     0x283ed59f1226dcefa7ff0ef87ceff5d5
+     0x1e2933215e1c3bba8d2b404d98f43086bfc6198a219b168f214042a5e7df6b22;
+   0xbe1f8593ee2c6f8f9497cc7335d97111
+     0x35a8e33503b3695be00528f8b82db931
+     0x1e2933215e1c3bba8d2b404d98f43086bfc6198a219b168f214042a5e7df6b21;
+ }
+ 
+ sqr {
+   0 0;
+   1 1;
+   3 5;
+   0x7fb838a8a0a95046b9d9d9fb4440f7bb
+     0x1555454005404440440044411100101445415141514155451010100055154545;
+   0x01f081e69f45d3254530766ab98d55fa612c7bb27ea31bc2621d894be9c0b196b3
+     0x0155004001541441551011510504111011050015141444454140511111554414010450154545041554440501455004140401514041104554415000450141144505;
+ }
+ 
+ div {
+   0 1 0 0;
+   0x207ccad257b4ed64447158315bfb9aca5cbc5622cfb8fcbb1380eea1bc5c624e
+     0x7fb838a8a0a95046b9d9d9fb4440f7bb
+     0xc1a7bd3b4e853fc92d4e1588719986aa 0;
+   0x6e0e2a282a5411ae76767ed1103deef069ef4ed3a14ff24b
+     0x5385621c6661aaa35a24150d2c08332e
+     0x01c2334cc957151dc7
+     0x398c4111da6d06cdf3d83704ee403101;
+ }
+ 
+ gcd {
+   0xc1a7bd3b4e853fc92d4e1588719986aa
+     0xbe1f8593ee2c6f8f9497cc7335d97111
+     3
+     0x283ed59f1226dcefa7ff0ef87ceff5d5
+     0x35a8e33503b3695be00528f8b82db931;
+   0xbe1f8593ee2c6f8f9497cc7335d97111
+     0xc1a7bd3b4e853fc92d4e1588719986aa
+     3
+     0x35a8e33503b3695be00528f8b82db931
+     0x283ed59f1226dcefa7ff0ef87ceff5d5;
+   0x800000000000000000000000000000000000000c9
+     0x3f0eba16286a2d57ea0991168d4994637e8343e36
+     1
+     0xa17e704470d80cb5a78f295db0ce543dda16a169
+     0x3c8c172e24598e90b9542e6b8f6571f54be572b50;
+ }
diff --cc tests/gfreduce
index 0000000,aec5318..806ec28
mode 000000,100644..100644
--- a/tests/gfreduce
+++ b/tests/gfreduce
@@@ -1,0 -1,61 +1,61 @@@
 -# $Id: gfreduce,v 1.1.2.1 2004/03/21 22:39:46 mdw Exp $
++# $Id: gfreduce,v 1.2 2004/03/21 22:52:06 mdw Exp $
+ #
+ # Test efficient polynomial reduction
+ 
+ reduce {
+           0x10000000
+   0x4509823098098435
+            0x8098435;
+   0x100000000000000050002
+     0x4509823098098435
+     0x4509823098098435;
+   0x100000000000000050002
+     0x450982309809843545609843098560803495
+     0x144f98a2f5cbc4773cfd;
+   0xb2ca471b0867d5fae2e4f27a2d2706da
+     0xf254423fef93d5d7a76ecf22c656c1352c53257875945d33
+     0x582f783fc210f72814780e69b0bd29ff;
+ }
+ 
+ modexp {
+   0x20000000000000000000000000000000000000000000000000000000000001001
+   0x02
+   0x1ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff
+   1;
+   0x20000000000000000000000000000000000000000000000000000000000001001
+   0x435932098459080438094509845
+   0x1ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff
+   1;
+   0x10000000000000000000000000000000000000000003
+   0x02
+   0x0fffffffffffffffffffffffffffffffffffffffffff
+   1;
+   0x10000000000000000000000000000000000000000003
+   0x34235950984598345900983409845690805680985
+   0x0fffffffffffffffffffffffffffffffffffffffffff
+   1;
+ }
+ 
+ sqrt {
+   0x20000000000000000000000000000000000000000000000000000000000001001
+   0x1f081e69f45d3254530766ab98d55fa612c7bb27ea31bc2621d894be9c0b196b3
+    0x7fb838a8a0a95046b9d9d9fb4440f7bbc1a7bd3b4e853fc92d4e1588719986aa;
+   0x10000000000000000000000000000000000000000003
+    0x4594094509835690805698083560980459903450984
+    0x820291881a244a02840a2f8ece3f23f88f38bf0b3a;
+ }
+ 
+ halftrace {
+   0x20000000000000000000000000000000000000000000000000000000000001001
+   0x174e65c7d14a8ec286df8c7df17662f13f1d3563f13c8c63f23f5d0bd5d1b45cd
+    0x8d68905434b020ccb849e17a03a5c441d2a104aaf523699c1cc7a93174d21d9d;
+ }
+ 
+ quadsolve {
+   0x20000000000000000000000000000000000000000000000000000000000001001
+   0x174e65c7d14a8ec286df8c7df17662f13f1d3563f13c8c63f23f5d0bd5d1b45cd
+    0x8d68905434b020ccb849e17a03a5c441d2a104aaf523699c1cc7a93174d21d9c;
+   0x10000000000000000000000000000000000000000003
+    0x3b818b447e90713da04f13c3b07cb5e2681d08e4700
+    0x27aa17c97dfa80bbdef9f91b243c6e6ddba1a223cac;
+ }
diff --cc tests/gfx
index 17bdbf7,866bd39..1ef6081
--- a/tests/gfx
+++ b/tests/gfx
@@@ -1,6 -1,6 +1,6 @@@
  # Test vectors for low-level GF functions
  #
- # $Id: gfx,v 1.1 2000/10/08 16:01:48 mdw Exp $
 -# $Id: gfx,v 1.1.4.1 2004/03/21 22:39:46 mdw Exp $
++# $Id: gfx,v 1.2 2004/03/21 22:52:06 mdw Exp $
  
  # --- Addition (and subtraction) ---