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

Merge and close elliptic curve branch.
---

c3caa2face1cda7002eb58245ad75865bf437455
diff --combined BRANCHES
index 2f32e82d,2f32e82d..53108159
--- 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 --combined Makefile.m4
index 814fd340,4a950442..9be48d51
--- 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'.
@@@ -272,7 -260,7 +281,7 @@@ define(`cipher_modes', `_(ecb) _(cbc) _
  
  define(`hashes', `dnl
  _(md5) _(md4) _(md2) _(tiger) dnl
 -_(sha) _(sha256) _(sha384) _(sha512) dnl
 +_(sha) _(sha224) _(sha256) _(sha384) _(sha512) dnl
  _(rmd128) _(rmd160) _(rmd256) _(rmd320)')
  define(`hash_modes', `_(mgf) _(hmac)')
  
@@@ -297,8 -285,7 +306,8 @@@ _(gfshare) _(gfx-sqr)'
  autoheaders: addsuffix(`gen_tables', `-tab.h') primetab.h mptypes.h
  define(`emit', `
  _item`'-tab.h: _item`'-mktab
 -	./_item`'-mktab >_item`'-tab.h')dnl
 +	./_item`'-mktab >_item`'-tab.h.new
 +	mv _item`'-tab.h.new _item`'-tab.h')dnl
  gen_tables
  
  primetab.h: primetab.c
@@@ -307,8 -294,7 +316,8 @@@ primetab.c: genprime
  		-t "unsigned short" -i primetab
  archinclude_HEADERS = mptypes.h
  mptypes.h: mptypes
 -	./mptypes >mptypes.h
 +	./mptypes >mptypes.h.new
 +	mv mptypes.h.new mptypes.h
  
  BUILT_SOURCES = \
  	getdate.c modes-stamp \
@@@ -321,24 -307,26 +330,26 @@@
  
  lib_LTLIBRARIES = libcatacomb.la
  
- libcatacomb_la_LDFLAGS = -version-info 2:1:0
+ libcatacomb_la_LDFLAGS = -version-info 3:0:1
  ## Middle number is the patchlevel.  Final number is the minor version.  The
  ## difference between the first and last numbers is major version.
  
  pkginclude_HEADERS = \
 -	arena.h paranoia.h \
 +	arena.h paranoia.h buf.h \
  	blkc.h hash.h gcipher.h ghash.h gmac.h grand.h ghash-def.h \
  	lcrand.h fibrand.h rc4.h seal.h rand.h noise.h fipstest.h maurer.h \
  	key.h key-data.h passphrase.h pixie.h lmem.h \
  	mpx.h bitops.h mpw.h mpscan.h mparena.h mp.h mptext.h mpint.h \
- 	exp.h mpbarrett.h mpmont.h mpcrt.h mprand.h mpmul.h \
- 	gfx.h \
+ 	exp.h mpbarrett.h mpbarrett-exp.h mpmont.h mpmont-exp.h \
+ 	mpcrt.h mprand.h mpmul.h \
+ 	gfx.h gf.h \
  	primetab.h pfilt.h rabin.h \
  	pgen.h prim.h strongprime.h limlee.h keycheck.h \
  	bbs.h rsa.h dh.h dsarand.h dsa.h \
  	oaep.h pkcs1.h pss.h tlsprf.h sslprf.h \
  	gfshare.h share.h \
  	rho.h \
+ 	field.h ec.h ec-exp.h \
  	allwithsuffix(`ciphers', `cipher_modes', `.h') \
  	allwithsuffix(`hashes', `hash_modes', `.h') \
  	addsuffix(`cipher_modes', `-def.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 \
@@@ -468,10 -459,7 +482,10 @@@ man_MANS = key.1 hashsum.1 keyring.5 pi
  
  EXTRA_DIST = \
  	Makefile.m4 genmodes $(man_MANS) xpixie \
 -	README.cipher README.hash README.random README.mp
 +	README.cipher README.hash README.random README.mp \
 +	debian/rules debian/copyright debian/control debian/changelog \
 +	debian/catacomb-bin.postinst debian/catacomb-bin.config \
 +	debian/catacomb-bin.prerm debian/catacomb-bin.templates
  
  dist-hook:
  	@ln getdate.c $(distdir) || ln $(srcdir)/getdate.c $(distdir) || true
@@@ -526,7 -514,13 +540,13 @@@ CTESTRIG(mpmont-mexp
  CTESTRIG(mpcrt)
  CTESTRIG(mpmul)
  CTESTRIG(gfx)
+ CTESTRIG(gfx-sqr)
  CTESTRIG(gfx-kmul)
+ CTESTRIG(gf-arith)
+ CTESTRIG(gf-gcd)
+ CTESTRIG(gfreduce)
+ CTESTRIG(ec-prime)
+ CTESTRIG(ec-bin)
  CTESTRIG(pgen)
  CTESTRIG(dsa-gen)
  CTESTRIG(dsa-sign)
@@@ -547,8 -541,7 +567,8 @@@ CLEANFILES = 
  ## --- Makefile building (haha!) ---
  
  $(srcdir)/Makefile.am: $(srcdir)/Makefile.m4
 -	m4 $(srcdir)/Makefile.m4 >$(srcdir)/Makefile.am
 +	m4 $(srcdir)/Makefile.m4 >$(srcdir)/Makefile.am.new
 +	mv $(srcdir)/Makefile.am.new $(srcdir)/Makefile.am
  
  DISTCLEANFILES = libtool
  
diff --combined calc/ec2.cal
index 00000000,0d7ceb91..3e890347
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 --combined calc/ecp.cal
index 04971aac,43ac1b3e..7c560c53
--- 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.
   *
@@@ -39,6 -45,7 +48,7 @@@
  
  obj ecp_curve { a, b, p };
  obj ecp_pt { x, y, e };
+ obj ecpp_pt { x, y, z, e };
  
  /*----- Main code ---------------------------------------------------------*/
  
@@@ -60,6 -67,72 +70,72 @@@ define ecp_pt(x, y, e
    return (p);
  }
  
+ define ecpp_pt(p)
+ {
+   local obj ecpp_pt pp;
+   if (istype(p, 1))
+     return (0);
+   pp.x = p.x;
+   pp.y = p.y;
+   pp.z = 1;
+   pp.e = p.e;
+   return (pp);
+ }
+ 
+ define ecpp_fix(pp)
+ {
+   local obj ecp_pt p;
+   local e, zi, z2, z3;
+   if (istype(pp, 1) || pp.z == 0)
+     return (0);
+   e = pp.e;
+   zi = minv(pp.z, e.p);
+   z2 = zi * zi;
+   z3 = zi * z2;
+   p.x = pp.x * z2 % e.p;
+   p.y = pp.y * z3 % e.p;
+   p.e = e;
+   return (p);
+ }
+ 
+ define ecpp_dbl(a)
+ {
+   local m, s, t, y2;
+   local e;
+   local obj ecpp_pt d;
+   if (istype(a, 1) || a.y == 0)
+     return (0);
+   e = a.e;
+   if (e.a % e.p == e.p - 3) {
+     m = a.z^3 % e.p;
+     m = 3 * (a.x + t4) * (a.x - t4) % e.p;
+   } else {
+     m = (3 * a.x^2 - e.a * a.z^4) % e.p;
+   }
+   d.z = 2 * a.y * a.z % e.p;
+   y2 = a.y^2 % e.p;
+   s = 4 * a.x * a.y % e.p;
+   d.x = (m^2 - 2 * s) % e.p;
+   d.y = (m * (s - d.x) - y * y2^2) % e.p;
+   d.e = e;
+   return (d);
+ }
+ 
+ define ecpp_add(a, b)
+ {
+   if (a == 0)
+     d = b;
+   else if (b == 0)
+     d = a;
+   else if (!istype(a, b))
+     quit "bad type arguments to ecp_pt_add";
+   else if (a.e != b.e)
+     quit "points from different curves in ecp_pt_add";
+   else {
+     e = a.e;
+     
+ }
+ 
  define ecp_pt_print(a)
  {
    print "(" : a.x : ", " : a.y : ")" :;
@@@ -96,6 -169,20 +172,20 @@@ define ecp_pt_add(a, b
    return (d);
  }
  
+ define ecp_pt_dbl(a)
+ {
+   local e, alpha;
+   local obj ecp_pt d;
+   if (istype(a, 1))
+     return (0);
+   e = a.e;
+   alpha = (3 * a.x^2 + e.a) * minv(2 * a.y, e.p) % e.p;
+   d.x = (alpha^2 - 2 * a.x) % e.p;
+   d.y = (-a.y + alpha * (a.x - d.x)) % e.p;
+   d.e = e;
+   return (d);
+ }
+ 
  define ecp_pt_neg(a)
  {
    local obj ecp_pt d;
@@@ -105,6 -192,15 +195,15 @@@
    return (d);
  }
  
+ define ecp_pt_check(a)
+ {
+   local e;
+ 
+   e = a.e;
+   if (a.y^2 % e.p != (a.x^3 + e.a * a.x + e.b) % e.p)
+     quit "bad curve point";
+ }
+ 
  define ecp_pt_mul(a, b)
  {
    local p, n;
@@@ -124,10 -220,18 +223,18 @@@
      if (n & 1)
        d += p;
      n >>= 1;
-     p += p;
+     p = ecp_pt_dbl(p);
    }
    return (d);
  }
  
+ /*----- FIPS186-2 standard curves -----------------------------------------*/
+ 
+ p192 = ecp_curve(-3, 0x64210519e59c80e70fa7e9ab72243049feb8deecc146b9b1,
+ 		 6277101735386680763835789423207666416083908700390324961279);
+ p192_r = 6277101735386680763835789423176059013767194773182842284081;
+ p192_g = ecp_pt(0x188da80eb03090f67cbf20eb43a18800f4ff0afd82ff1012,
+ 		0x07192b95ffc8da78631011ed6b24cdd573f977a11e794811, p192);
+ 
  /*----- That's all, folks -------------------------------------------------*/
  
diff --combined calc/gfx.cal
index 8d8fd00e,5b19cb32..446061ed
--- 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.
   *
@@@ -79,7 -82,9 +85,9 @@@ define gf_mul(x, y
  define gfx_div(rx, dx)
  {
    local r = gfint(rx), d = gfint(dx), i;
-   local q = 0, dbits = highbit(d), rbits = highbit(r);
+   local q = 0, dbits, rbits;
+   dbits = highbit(d);
+   rbits = highbit(r);
    for (i = rbits - dbits; i >= 0; i--) {
      if (bit(r, i + dbits)) {
        r = xor(r, d << i);
@@@ -91,14 -96,36 +99,36 @@@
  
  define gf_div(x, y)
  {
-   local l = gfx_div(x, y);
+   local l;
+   l = gfx_div(x, y);
    return gf(l[[0]]);
  }
  
  define gf_mod(x, y)
  {
-   local l = gfx_div(x, y);
+   local l;
+   l = gfx_div(x, y);
    return gf(l[[1]]);
  }
  
+ define gf_inv(a, b)
+ {
+   local g, x, y, X, Y, u, v, t, q, r;
+   x = gf(1); X = gf(0);
+   y = gf(0); Y = gf(1);
+ 
+   if (b == gf(0)) { g = a; } else if (a == gf(0)) { g = b; }
+   else {
+     while (b != gf(0)) {
+       q = gf_div(b, a); r = gf_mod(b, a);
+       t = X * q + x; x = X; X = t;
+       t = Y * q + y; y = Y; Y = t;
+       b = a; a = r;
+     }
+     g = a;
+   }
+   if (g != gf(1)) quit "not coprime in gf_inv";
+   return Y;
+ }
+ 
  /*----- That's all, folks -------------------------------------------------*/
diff --combined configure.in
index 555e0efb,599cc282..7e2245d9
--- 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 --combined debian/changelog
index 5717a50f,00000000..ce65b69c
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 --combined ec-bin.c
index 00000000,4f79c3dc..3e85e65b
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 --combined ec-exp.h
index 9ba7bc5e,bb4e08a2..0daf7171
--- 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.
   *
@@@ -58,6 -61,7 +64,7 @@@
  
  #define EXP_MUL(a, x) EC_ADD(c, &(a), &(a), &(x))
  #define EXP_SQR(a) EC_DBL(c, &(a), &(a));
+ #define EXP_FIX(x) EC_FIX(c, &(x), &(x));
  
  #define EXP_SETMUL(d, x, y) do {					\
    EC_CREATE(&(d));							\
diff --combined ec-prime.c
index 46118552,40f487e4..bdc63683
--- 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.
   *
@@@ -54,14 -63,37 +66,37 @@@ typedef struct ecctx 
    mp *a, *b;
  } ecctx;
  
- /*----- Main code ---------------------------------------------------------*/
+ /*----- Simple prime curves -----------------------------------------------*/
  
- static const ec_ops ec_primeops;
+ static const ec_ops ec_primeops, ec_primeprojops, ec_primeprojxops;
  
  static ec *ecneg(ec_curve *c, ec *d, const ec *p)
  {
    EC_COPY(d, p);
-   d->y = F_NEG(c->f, d->y, d->y);
+   if (d->y)
+     d->y = F_NEG(c->f, d->y, d->y);
+   return (d);
+ }
+ 
+ static ec *ecfind(ec_curve *c, ec *d, mp *x)
+ {
+   mp *p, *q;
+   ecctx *cc = (ecctx *)c;
+   field *f = c->f;
+ 
+   q = F_SQR(f, MP_NEW, x);
+   p = F_MUL(f, MP_NEW, x, q);
+   q = F_MUL(f, q, x, cc->a);
+   p = F_ADD(f, p, p, q);
+   p = F_ADD(f, p, p, cc->b);
+   MP_DROP(q);
+   p = F_SQRT(f, p, p);
+   if (!p)
+     return (0);
+   EC_DESTROY(d);
+   d->x = MP_COPY(x);
+   d->y = p;
+   d->z = MP_COPY(f->one);
    return (d);
  }
  
@@@ -69,7 -101,7 +104,7 @@@ static ec *ecdbl(ec_curve *c, ec *d, co
  {
    if (EC_ATINF(a))
      EC_SETINF(d);
-   else if (!MP_LEN(a->y))
+   else if (F_ZEROP(c->f, a->y))
      EC_COPY(d, a);
    else {
      field *f = c->f;
@@@ -77,19 -109,19 +112,19 @@@
      mp *lambda;
      mp *dy, *dx;
  
-     dx = F_SQR(f, MP_NEW, a->x);
-     dy = F_DBL(f, MP_NEW, a->y);
-     dx = F_TPL(f, dx, dx);
-     dx = F_ADD(f, dx, dx, cc->a);
-     dy = F_INV(f, dy, dy);
-     lambda = F_MUL(f, MP_NEW, dx, dy);
+     dx = F_SQR(f, MP_NEW, a->x);	/* %$x^2$% */
+     dy = F_DBL(f, MP_NEW, a->y);	/* %$2 y$% */
+     dx = F_TPL(f, dx, dx);		/* %$3 x^2$% */
+     dx = F_ADD(f, dx, dx, cc->a);	/* %$3 x^2 + A$% */
+     dy = F_INV(f, dy, dy);		/* %$(2 y)^{-1}$% */
+     lambda = F_MUL(f, MP_NEW, dx, dy);	/* %$\lambda = (3 x^2 + A)/(2 y)$% */
  
-     dx = F_SQR(f, dx, lambda);
-     dy = F_DBL(f, dy, a->x);
-     dx = F_SUB(f, dx, dx, dy);
-     dy = F_SUB(f, dy, a->x, dx);
-     dy = F_MUL(f, dy, lambda, dy);
-     dy = F_SUB(f, dy, dy, a->y);
+     dx = F_SQR(f, dx, lambda);		/* %$\lambda^2$% */
+     dy = F_DBL(f, dy, a->x);		/* %$2 x$% */
+     dx = F_SUB(f, dx, dx, dy);		/* %$x' = \lambda^2 - 2 x */
+     dy = F_SUB(f, dy, a->x, dx);	/* %$x - x'$% */
+     dy = F_MUL(f, dy, lambda, dy);	/* %$\lambda (x - x')$% */
+     dy = F_SUB(f, dy, dy, a->y);	/* %$y' = \lambda (x - x') - y$% */
  
      EC_DESTROY(d);
      d->x = dx;
@@@ -100,6 -132,94 +135,94 @@@
    return (d);
  }
  
+ static ec *ecprojdbl(ec_curve *c, ec *d, const ec *a)
+ {
+   if (EC_ATINF(a))
+     EC_SETINF(d);
+   else if (F_ZEROP(c->f, a->y))
+     EC_COPY(d, a);
+   else {
+     field *f = c->f;
+     ecctx *cc = (ecctx *)c;
+     mp *p, *q, *m, *s, *dx, *dy, *dz;
+ 
+     p = F_SQR(f, MP_NEW, a->z);		/* %$z^2$% */
+     q = F_SQR(f, MP_NEW, p);		/* %$z^4$% */
+     p = F_MUL(f, p, q, cc->a);		/* %$A z^4$% */
+     m = F_SQR(f, MP_NEW, a->x);		/* %$x^2$% */
+     m = F_TPL(f, m, m);			/* %$3 x^2$% */
+     m = F_ADD(f, m, m, p);		/* %$m = 3 x^2 + A z^4$% */
+ 
+     q = F_DBL(f, q, a->y);		/* %$2 y$% */
+     dz = F_MUL(f, MP_NEW, q, a->z);	/* %$z' = 2 y z$% */
+ 
+     p = F_SQR(f, p, q);			/* %$4 y^2$% */
+     s = F_MUL(f, MP_NEW, p, a->x);	/* %$s = 4 x y^2$% */
+     q = F_SQR(f, q, p);			/* %$16 y^4$% */
+     q = F_HLV(f, q, q);			/* %$t = 8 y^4$% */
+ 
+     p = F_DBL(f, p, s);			/* %$2 s$% */
+     dx = F_SQR(f, MP_NEW, m);		/* %$m^2$% */
+     dx = F_SUB(f, dx, dx, p);		/* %$x' = m^2 - 2 s$% */
+ 
+     s = F_SUB(f, s, s, dx);		/* %$s - x'$% */
+     dy = F_MUL(f, p, m, s);		/* %$m (s - x')$% */
+     dy = F_SUB(f, dy, dy, q);		/* %$y' = m (s - x') - t$% */
+ 
+     EC_DESTROY(d);
+     d->x = dx;
+     d->y = dy;
+     d->z = dz;
+     MP_DROP(m);
+     MP_DROP(q);
+     MP_DROP(s);
+   }
+   return (d);
+ }
+ 
+ static ec *ecprojxdbl(ec_curve *c, ec *d, const ec *a)
+ {
+   if (EC_ATINF(a))
+     EC_SETINF(d);
+   else if (F_ZEROP(c->f, a->y))
+     EC_COPY(d, a);
+   else {
+     field *f = c->f;
+     mp *p, *q, *m, *s, *dx, *dy, *dz;
+ 
+     m = F_SQR(f, MP_NEW, a->z);		/* %$z^2$% */
+     p = F_SUB(f, MP_NEW, a->x, m);	/* %$x - z^2$% */
+     q = F_ADD(f, MP_NEW, a->x, m);	/* %$x + z^2$% */
+     m = F_MUL(f, m, p, q);		/* %$x^2 - z^4$% */
+     m = F_TPL(f, m, m);			/* %$m = 3 x^2 - 3 z^4$% */
+ 
+     q = F_DBL(f, q, a->y);		/* %$2 y$% */
+     dz = F_MUL(f, MP_NEW, q, a->z);	/* %$z' = 2 y z$% */
+ 
+     p = F_SQR(f, p, q);			/* %$4 y^2$% */
+     s = F_MUL(f, MP_NEW, p, a->x);	/* %$s = 4 x y^2$% */
+     q = F_SQR(f, q, p);			/* %$16 y^4$% */
+     q = F_HLV(f, q, q);			/* %$t = 8 y^4$% */
+ 
+     p = F_DBL(f, p, s);			/* %$2 s$% */
+     dx = F_SQR(f, MP_NEW, m);		/* %$m^2$% */
+     dx = F_SUB(f, dx, dx, p);		/* %$x' = m^2 - 2 s$% */
+ 
+     s = F_SUB(f, s, s, dx);		/* %$s - x'$% */
+     dy = F_MUL(f, p, m, s);		/* %$m (s - x')$% */
+     dy = F_SUB(f, dy, dy, q);		/* %$y' = m (s - x') - t$% */
+ 
+     EC_DESTROY(d);
+     d->x = dx;
+     d->y = dy;
+     d->z = dz;
+     MP_DROP(m);
+     MP_DROP(q);
+     MP_DROP(s);
+   }
+   return (d);
+ }
+ 
  static ec *ecadd(ec_curve *c, ec *d, const ec *a, const ec *b)
  {
    if (a == b)
@@@ -114,29 -234,31 +237,31 @@@
      mp *dy, *dx;
  
      if (!MP_EQ(a->x, b->x)) {
-       dy = F_SUB(f, MP_NEW, a->y, b->y);
-       dx = F_SUB(f, MP_NEW, a->x, b->x);
-       dx = F_INV(f, dx, dx);
+       dy = F_SUB(f, MP_NEW, a->y, b->y); /* %$y_0 - y_1$% */
+       dx = F_SUB(f, MP_NEW, a->x, b->x); /* %$x_0 - x_1$% */
+       dx = F_INV(f, dx, dx);		/* %$(x_0 - x_1)^{-1}$% */
        lambda = F_MUL(f, MP_NEW, dy, dx);
-     } else if (!MP_LEN(a->y) || !MP_EQ(a->y, b->y)) {
+ 				   /* %$\lambda = (y_0 - y1)/(x_0 - x_1)$% */
+     } else if (F_ZEROP(c->f, a->y) || !MP_EQ(a->y, b->y)) {
        EC_SETINF(d);
        return (d);
      } else {
        ecctx *cc = (ecctx *)c;
-       dx = F_SQR(f, MP_NEW, a->x);
-       dx = F_TPL(f, dx, dx);
-       dx = F_ADD(f, dx, dx, cc->a);
-       dy = F_DBL(f, MP_NEW, a->y);
-       dy = F_INV(f, dy, dy);
+       dx = F_SQR(f, MP_NEW, a->x);	/* %$x_0^2$% */
+       dx = F_TPL(f, dx, dx);		/* %$3 x_0^2$% */
+       dx = F_ADD(f, dx, dx, cc->a);	/* %$3 x_0^2 + A$% */
+       dy = F_DBL(f, MP_NEW, a->y);	/* %$2 y_0$% */
+       dy = F_INV(f, dy, dy);		/* %$(2 y_0)^{-1}$% */
        lambda = F_MUL(f, MP_NEW, dx, dy);
+ 				    /* %$\lambda = (3 x_0^2 + A)/(2 y_0)$% */
      }
  
-     dx = F_SQR(f, dx, lambda);
-     dx = F_SUB(f, dx, dx, a->x);
-     dx = F_SUB(f, dx, dx, b->x);
-     dy = F_SUB(f, dy, b->x, dx);
-     dy = F_MUL(f, dy, lambda, dy);
-     dy = F_SUB(f, dy, dy, b->y);
+     dx = F_SQR(f, dx, lambda);		/* %$\lambda^2$% */
+     dx = F_SUB(f, dx, dx, a->x);	/* %$\lambda^2 - x_0$% */
+     dx = F_SUB(f, dx, dx, b->x);	/* %$x' = \lambda^2 - x_0 - x_1$% */
+     dy = F_SUB(f, dy, b->x, dx);	/* %$x_1 - x'$% */
+     dy = F_MUL(f, dy, lambda, dy);	/* %$\lambda (x_1 - x')$% */
+     dy = F_SUB(f, dy, dy, b->y);      /* %$y' = \lambda (x_1 - x') - y_1$% */
  
      EC_DESTROY(d);
      d->x = dx;
@@@ -147,6 -269,98 +272,98 @@@
    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;
+     mp *p, *q, *r, *w, *u, *s, *dx, *dy, *dz;
+ 
+     q = F_SQR(f, MP_NEW, a->z);		/* %$z_0^2$% */
+     u = F_MUL(f, MP_NEW, q, b->x);	/* %$u = x_1 z_0^2$% */
+     p = F_MUL(f, MP_NEW, q, b->y);	/* %$y_1 z_0^2$% */
+     s = F_MUL(f, q, p, a->z);		/* %$s = y_1 z_0^3$% */
+ 
+     w = F_SUB(f, p, a->x, u);		/* %$w = x_0 - u$% */
+     r = F_SUB(f, MP_NEW, a->y, s);	/* %$r = y_0 - s$% */
+     if (F_ZEROP(f, w)) {
+       MP_DROP(w);
+       MP_DROP(u);
+       MP_DROP(s);
+       if (F_ZEROP(f, r)) {
+ 	MP_DROP(r);
+ 	return (c->ops->dbl(c, d, a));
+       } else {
+ 	MP_DROP(r);
+ 	EC_SETINF(d);
+ 	return (d);
+       }
+     }
+     u = F_ADD(f, u, u, a->x);		/* %$t = x_0 + u$% */
+     s = F_ADD(f, s, s, a->y);		/* %$m = y_0 + r$% */
+ 
+     dz = F_MUL(f, MP_NEW, a->z, w);	/* %$z' = z_0 w$% */
+ 
+     p = F_SQR(f, MP_NEW, w);		/* %$w^2$% */
+     q = F_MUL(f, MP_NEW, p, u);		/* %$t w^2$% */
+     u = F_MUL(f, u, p, w);		/* %$w^3$% */
+     p = F_MUL(f, p, u, s);		/* %$m w^3$% */
+     
+     dx = F_SQR(f, u, r);		/* %$r^2$% */
+     dx = F_SUB(f, dx, dx, q);		/* %$x' = r^2 - t w^2$% */
+ 
+     s = F_DBL(f, s, dx);		/* %$2 x'$% */
+     q = F_SUB(f, q, q, s);		/* %$v = t w^2 - 2 x'$% */
+     dy = F_MUL(f, s, q, r);		/* %$v r$% */
+     dy = F_SUB(f, dy, dy, p);		/* %$v r - m w^3$% */
+     dy = F_HLV(f, dy, dy);		/* %$y' = (v r - m w^3)/2$% */
+ 
+     EC_DESTROY(d);
+     d->x = dx;
+     d->y = dy;
+     d->z = dz;
+     MP_DROP(p);
+     MP_DROP(q);
+     MP_DROP(r);
+     MP_DROP(w);
+   }
+   return (d);
+ }
+ 
+ static int eccheck(ec_curve *c, const ec *p)
+ {
+   ecctx *cc = (ecctx *)c;
+   field *f = c->f;
+   int rc;
+   mp *l = F_SQR(f, MP_NEW, p->y);
+   mp *x = F_SQR(f, MP_NEW, p->x);
+   mp *r = F_MUL(f, MP_NEW, x, p->x);
+   x = F_MUL(f, x, cc->a, p->x);
+   r = F_ADD(f, r, r, x);
+   r = F_ADD(f, r, r, cc->b);
+   rc = MP_EQ(l, r) ? 0 : -1;
+   mp_drop(l);
+   mp_drop(x);
+   mp_drop(r);
+   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;
@@@ -157,7 -371,7 +374,7 @@@
  
  /* --- @ec_prime@, @ec_primeproj@ --- *
   *
-  * Arguments:	@field *f@ = the underyling field for this elliptic curve
+  * Arguments:	@field *f@ = the underlying field for this elliptic curve
   *		@mp *a, *b@ = the coefficients for this curve
   *
   * Returns:	A pointer to the curve.
@@@ -172,13 -386,42 +389,42 @@@ extern ec_curve *ec_prime(field *f, mp 
    ecctx *cc = CREATE(ecctx);
    cc->c.ops = &ec_primeops;
    cc->c.f = f;
-   cc->a = MP_COPY(a);
-   cc->b = MP_COPY(b);
+   cc->a = F_IN(f, MP_NEW, a);
+   cc->b = F_IN(f, MP_NEW, b);
+   return (&cc->c);
+ }
+ 
+ extern ec_curve *ec_primeproj(field *f, mp *a, mp *b)
+ {
+   ecctx *cc = CREATE(ecctx);
+   mp *ax;
+ 
+   ax = mp_add(MP_NEW, a, MP_THREE);
+   ax = F_IN(f, ax, ax);
+   if (F_ZEROP(f, ax))
+     cc->c.ops = &ec_primeprojxops;
+   else
+     cc->c.ops = &ec_primeprojops;
+   MP_DROP(ax);
+   cc->c.f = f;
+   cc->a = F_IN(f, MP_NEW, a);
+   cc->b = F_IN(f, MP_NEW, b);
    return (&cc->c);
  }
  
  static const ec_ops ec_primeops = {
-   ecdestroy, ec_idin, ec_idout, 0, ecneg, ecadd, ec_stdsub, ecdbl
+   ecdestroy, ec_idin, ec_idout, ec_idfix,
+   0, ecneg, ecadd, ec_stdsub, ecdbl, eccheck
+ };
+ 
+ static const ec_ops ec_primeprojops = {
+   ecdestroy, ec_projin, ec_projout, ec_projfix,
+   0, ecneg, ecprojadd, ec_stdsub, ecprojdbl, ecprojcheck
+ };
+ 
+ static const ec_ops ec_primeprojxops = {
+   ecdestroy, ec_projin, ec_projout, ec_projfix,
+   0, ecneg, ecprojadd, ec_stdsub, ecprojxdbl, ecprojcheck
  };
  
  /*----- Test rig ----------------------------------------------------------*/
@@@ -187,33 -430,48 +433,48 @@@
  
  #define MP(x) mp_readstring(MP_NEW, #x, 0, 0)
  
- int main(void)
+ 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-prime: ");
+   fflush(stdout);
    a = MP(-3);
    b = MP(0x64210519e59c80e70fa7e9ab72243049feb8deecc146b9b1);
    p = MP(6277101735386680763835789423207666416083908700390324961279);
-   r = MP(6277101735386680763835789423176059013767194773182842284081);
+   r = MP(6277101735386680763835789423176059013767194773182842284080);
  
    f = field_prime(p);
-   c = ec_prime(f, a, b);
+   c = ec_primeproj(f, a, b);
    
    g.x = MP(0x188da80eb03090f67cbf20eb43a18800f4ff0afd82ff1012);
    g.y = MP(0x07192b95ffc8da78631011ed6b24cdd573f977a11e794811);
  
-   ec_mul(c, &d, &g, r);
-   MP_PRINT("d.x", d.x);
-   MP_PRINT("d.y", d.y);
- 
-   ec_destroy(&d);
+   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_EPRINT("d.x", d.x);
+       MP_EPRINT("d.y", 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);
  }
  
diff --combined ec.c
index 47c01a9f,a2b229f7..c95333f6
--- 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.
   *
@@@ -110,7 -116,7 +119,7 @@@ ec *ec_copy(ec *d, const ec *p) { EC_CO
  
  /*----- Standard curve operations -----------------------------------------*/
  
- /* --- @ec_idin@, @ec_idout@ --- *
+ /* --- @ec_idin@, @ec_idout@, @ec_idfix@ --- *
   *
   * Arguments:	@ec_curve *c@ = pointer to an elliptic curve
   *		@ec *d@ = pointer to the destination
@@@ -149,6 -155,12 +158,12 @@@ ec *ec_idout(ec_curve *c, ec *d, const 
    return (d);
  }
  
+ ec *ec_idfix(ec_curve *c, ec *d, const ec *p)
+ {
+   EC_COPY(d, p);
+   return (d);
+ }
+ 
  /* --- @ec_projin@, @ec_projout@ --- *
   *
   * Arguments:	@ec_curve *c@ = pointer to an elliptic curve
@@@ -179,12 -191,15 +194,15 @@@ ec *ec_projout(ec_curve *c, ec *d, cons
    if (EC_ATINF(p))
      EC_SETINF(d);
    else {
-     mp *x, *y, *z;
+     mp *x, *y, *z, *zz;
      field *f = c->f;
      z = F_INV(f, MP_NEW, p->z);
-     x = F_MUL(f, d->x, p->x, z);
+     zz = F_SQR(f, MP_NEW, z);
+     z = F_MUL(f, z, zz, z);
+     x = F_MUL(f, d->x, p->x, zz);
      y = F_MUL(f, d->y, p->y, z);
      mp_drop(z);
+     mp_drop(zz);
      mp_drop(d->z);
      d->x = F_OUT(f, x, x);
      d->y = F_OUT(f, y, y);
@@@ -193,6 -208,28 +211,28 @@@
    return (d);
  }
  
+ ec *ec_projfix(ec_curve *c, ec *d, const ec *p)
+ {
+   if (EC_ATINF(p))
+     EC_SETINF(d);
+   else if (d->z == c->f->one)
+     EC_COPY(d, p);
+   else {
+     mp *z, *zz;
+     field *f = c->f;
+     z = F_INV(f, MP_NEW, p->z);
+     zz = F_SQR(f, MP_NEW, z);
+     z = F_MUL(f, z, zz, z);
+     d->x = F_MUL(f, d->x, p->x, zz);
+     d->y = F_MUL(f, d->y, p->y, z);
+     mp_drop(z);
+     mp_drop(zz);
+     mp_drop(d->z);
+     d->z = MP_COPY(f->one);
+   }
+   return (d);  
+ }
+ 
  /* --- @ec_stdsub@ --- *
   *
   * Arguments:	@ec_curve *c@ = pointer to an elliptic curve
@@@ -210,6 -247,7 +250,7 @@@ ec *ec_stdsub(ec_curve *c, ec *d, cons
  {
    ec t = EC_INIT;
    EC_NEG(c, &t, q);
+   EC_FIX(c, &t, &t);
    EC_ADD(c, d, p, &t);
    EC_DESTROY(&t);
    return (d);
@@@ -246,10 -284,28 +287,28 @@@ ec *ec_find(ec_curve *c, ec *d, mp *x
    x = F_IN(c->f, MP_NEW, x);
    if ((d = EC_FIND(c, d, x)) != 0)
      EC_OUT(c, d, d);
-   mp_drop(x);
+   MP_DROP(x);
    return (d);
  }
  
+ /* --- @ec_neg@ --- *
+  *
+  * Arguments:	@ec_curve *c@ = pointer to an elliptic curve
+  *		@ec *d@ = pointer to the destination point
+  *		@const ec *p@ = pointer to the operand point
+  *
+  * Returns:	The destination point.
+  *
+  * Use:		Computes the negation of the given point.
+  */
+ 
+ ec *ec_neg(ec_curve *c, ec *d, const ec *p)
+ {
+   EC_IN(c, d, p);
+   EC_NEG(c, d, d);
+   return (EC_OUT(c, d, d));
+ }
+ 
  /* --- @ec_add@ --- *
   *
   * Arguments:	@ec_curve *c@ = pointer to an elliptic curve
@@@ -273,6 -329,29 +332,29 @@@ ec *ec_add(ec_curve *c, ec *d, const e
    return (d);
  }
  
+ /* --- @ec_sub@ --- *
+  *
+  * Arguments:	@ec_curve *c@ = pointer to an elliptic curve
+  *		@ec *d@ = pointer to the destination point
+  *		@const ec *p, *q@ = pointers to the operand points
+  *
+  * Returns:	The destination @d@.
+  *
+  * Use:		Subtracts one point from another on an elliptic curve.
+  */
+ 
+ ec *ec_sub(ec_curve *c, ec *d, const ec *p, const ec *q)
+ {
+   ec pp, qq;
+   EC_IN(c, &pp, p);
+   EC_IN(c, &qq, q);
+   EC_SUB(c, d, &qq, &qq);
+   EC_OUT(c, d, d);
+   EC_DESTROY(&pp);
+   EC_DESTROY(&qq);
+   return (d);
+ }
+ 
  /* --- @ec_dbl@ --- *
   *
   * Arguments:	@ec_curve *c@ = pointer to an elliptic curve
@@@ -291,6 -370,29 +373,29 @@@ ec *ec_dbl(ec_curve *c, ec *d, const e
    return (EC_OUT(c, d, d));
  }
  
+ /* --- @ec_check@ --- *
+  *
+  * Arguments:	@ec_curve *c@ = pointer to an elliptic curve
+  *		@const ec *p@ = pointer to the point
+  *
+  * Returns:	Zero if OK, nonzero if this is an invalid point.
+  *
+  * Use:		Checks that a point is actually on an elliptic curve.
+  */
+ 
+ int ec_check(ec_curve *c, const ec *p)
+ {
+   ec t = EC_INIT;
+   int rc;
+ 
+   if (EC_ATINF(p))
+     return (0);
+   EC_IN(c, &t, p);
+   rc = EC_CHECK(c, &t);
+   EC_DESTROY(&t);
+   return (rc);
+ }
+ 
  /* --- @ec_imul@, @ec_mul@ --- *
   *
   * Arguments:	@ec_curve *c@ = pointer to an elliptic curve
@@@ -316,10 -418,15 +421,15 @@@ ec *ec_imul(ec_curve *c, ec *d, const e
    EC_SETINF(d);
    if (MP_LEN(n) == 0)
      ;
-   else if (MP_LEN(n) < EXP_THRESH)
-     EXP_SIMPLE(*d, t, n);
-   else
-     EXP_WINDOW(*d, t, n);
+   else {
+     if (n->f & MP_NEG)
+       EC_NEG(c, &t, &t);
+     if (MP_LEN(n) < EXP_THRESH)
+       EXP_SIMPLE(*d, t, n);
+     else
+       EXP_WINDOW(*d, t, n);
+   }
+   EC_DESTROY(&t);
    return (d);
  }
  
diff --combined ec.h
index 105838cc,bb5e2017..07f1468e
--- 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.
   *
@@@ -80,27 -89,37 +92,37 @@@ typedef struct ec_mulfactor 
    mp *exp;				/* The exponent */
  } ec_mulfactor;
  
- /* --- Elliptic curve operations --- */
+ /* --- Elliptic curve operations --- *
+  *
+  * All operations (apart from @destroy@ and @in@) are guaranteed to be
+  * performed on internal representations of points.  Moreover, the second
+  * argument to @add@ and @mul@ is guaranteed to be the output of @in@ or
+  * @fix@.
+  */
  
  typedef struct ec_ops {
    void (*destroy)(ec_curve */*c*/);
    ec *(*in)(ec_curve */*c*/, ec */*d*/, const ec */*p*/);
    ec *(*out)(ec_curve */*c*/, ec */*d*/, const ec */*p*/);
+   ec *(*fix)(ec_curve */*c*/, ec */*d*/, const ec */*p*/);
    ec *(*find)(ec_curve */*c*/, ec */*d*/, mp */*x*/);
    ec *(*neg)(ec_curve */*c*/, ec */*d*/, const ec */*p*/);
    ec *(*add)(ec_curve */*c*/, ec */*d*/, const ec */*p*/, const ec */*q*/);
    ec *(*sub)(ec_curve */*c*/, ec */*d*/, const ec */*p*/, const ec */*q*/);
    ec *(*dbl)(ec_curve */*c*/, ec */*d*/, const ec */*p*/);
+   int (*check)(ec_curve */*c*/, const ec */*p*/);
  } ec_ops;
  
  #define EC_IN(c, d, p)		(c)->ops->in((c), (d), (p))
- #define EC_OUT(c, d, p)		(c)->ops->in((c), (d), (p))
+ #define EC_OUT(c, d, p)		(c)->ops->out((c), (d), (p))
+ #define EC_FIX(c, d, p)		(c)->ops->fix((c), (d), (p))
  
  #define EC_FIND(c, d, x)	(c)->ops->find((c), (d), (x))
  #define EC_NEG(c, d, x)		(c)->ops->neg((c), (d), (x))
  #define EC_ADD(c, d, p, q)	(c)->ops->add((c), (d), (p), (q))
  #define EC_SUB(c, d, p, q)	(c)->ops->sub((c), (d), (p), (q))
  #define EC_DBL(c, d, p)		(c)->ops->dbl((c), (d), (p))
+ #define EC_CHECK(c, p)		(c)->ops->check((c), (p))
  
  /*----- Simple memory management things -----------------------------------*/
  
@@@ -207,32 -226,6 +229,6 @@@ extern ec *ec_copy(ec */*d*/, const ec 
  
  /*----- Interesting arithmetic --------------------------------------------*/
  
- /* --- @ec_in@ --- *
-  *
-  * Arguments:	@ec_curve *c@ = pointer to an elliptic curve
-  *		@ec *d@ = pointer to the destination point
-  *		@const ec *p@ = pointer to the source point
-  *
-  * Returns:	The destination point.
-  *
-  * Use:		Converts a point to internal representation.
-  */
- 
- extern ec *ec_in(ec_curve */*c*/, ec */*d*/, const ec */*p*/);
- 
- /* --- @ec_out@ --- *
-  *
-  * Arguments:	@ec_curve *c@ = pointer to an elliptic curve
-  *		@ec *d@ = pointer to the destination point
-  *		@const ec *p@ = pointer to the source point
-  *
-  * Returns:	The destination point.
-  *
-  * Use:		Converts a point to external representation.
-  */
- 
- extern ec *ec_out(ec_curve */*c*/, ec */*d*/, const ec */*p*/);
- 
  /* --- @ec_find@ --- *
   *
   * Arguments:	@ec_curve *c@ = pointer to an elliptic curve
@@@ -303,6 -296,18 +299,18 @@@ extern ec *ec_sub(ec_curve */*c*/, ec *
  
  extern ec *ec_dbl(ec_curve */*c*/, ec */*d*/, const ec */*p*/);
  
+ /* --- @ec_check@ --- *
+  *
+  * Arguments:	@ec_curve *c@ = pointer to an elliptic curve
+  *		@const ec *p@ = pointer to the point
+  *
+  * Returns:	Zero if OK, nonzero if this is an invalid point.
+  *
+  * Use:		Checks that a point is actually on an elliptic curve.
+  */
+ 
+ extern int ec_check(ec_curve */*c*/, const ec */*p*/);
+ 
  /* --- @ec_mul@, @ec_imul@ --- *
   *
   * Arguments:	@ec_curve *c@ = pointer to an elliptic curve
@@@ -340,7 -345,7 +348,7 @@@ extern ec *ec_immul(ec_curve */*c*/, e
  
  /*----- Standard curve operations -----------------------------------------*/
  
- /* --- @ec_idin@, @ec_idout@ --- *
+ /* --- @ec_idin@, @ec_idout@, @ec_idfix@ --- *
   *
   * Arguments:	@ec_curve *c@ = pointer to an elliptic curve
   *		@ec *d@ = pointer to the destination
@@@ -355,8 -360,9 +363,9 @@@
  
  extern ec *ec_idin(ec_curve */*c*/, ec */*d*/, const ec */*p*/);
  extern ec *ec_idout(ec_curve */*c*/, ec */*d*/, const ec */*p*/);
+ extern ec *ec_idfix(ec_curve */*c*/, ec */*d*/, const ec */*p*/);
  
- /* --- @ec_projin@, @ec_projout@ --- *
+ /* --- @ec_projin@, @ec_projout@, @ec_projfix@ --- *
   *
   * Arguments:	@ec_curve *c@ = pointer to an elliptic curve
   *		@ec *d@ = pointer to the destination
@@@ -370,6 -376,7 +379,7 @@@
  
  extern ec *ec_projin(ec_curve */*c*/, ec */*d*/, const ec */*p*/);
  extern ec *ec_projout(ec_curve */*c*/, ec */*d*/, const ec */*p*/);
+ extern ec *ec_projfix(ec_curve */*c*/, ec */*d*/, const ec */*p*/);
  
  /* --- @ec_stdsub@ --- *
   *
@@@ -402,7 -409,7 +412,7 @@@ extern void ec_destroycurve(ec_curve */
  
  /* --- @ec_prime@, @ec_primeproj@ --- *
   *
-  * Arguments:	@field *f@ = the underyling field for this elliptic curve
+  * Arguments:	@field *f@ = the underlying field for this elliptic curve
   *		@mp *a, *b@ = the coefficients for this curve
   *
   * Returns:	A pointer to the curve.
@@@ -415,18 -422,20 +425,20 @@@
  extern ec_curve *ec_prime(field */*f*/, mp */*a*/, mp */*b*/);
  extern ec_curve *ec_primeproj(field */*f*/, mp */*a*/, mp */*b*/);
  
- /* --- @ec_bin@ --- *
+ /* --- @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 a non-supersingular elliptic
-  *		curve defined over a binary field.
+  * 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.
   */
  
  extern ec_curve *ec_bin(field */*f*/, mp */*a*/, mp */*b*/);
+ extern ec_curve *ec_binproj(field */*f*/, mp */*a*/, mp */*b*/);
  
  /*----- That's all, folks -------------------------------------------------*/
  
diff --combined exp.h
index 6cfdfd82,fc9e3a92..6bfd686b
--- 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.
@@@ -99,6 -102,10 +105,10 @@@ typedef struct exp_simul 
   * @EXP_MUL(a, x)@		Multiplies @a@ by @x@ (writing the result
   *				back to @a@).
   *
+  * @EXP_FIX(x)@			Makes @x@ be a canonical representation of
+  *				its value.  All multiplications have the
+  *				right argument canonical.
+  *
   * @EXP_SQR(a)@			Multiplies @a@ by itself.
   *
   * @EXP_SETMUL(d, x, y)@	Sets @d@ to be the product of @x@ and @y@.
@@@ -140,6 -147,7 +150,7 @@@
  									\
    /* --- Do the main body of the work --- */				\
  									\
+   EXP_FIX(g);								\
    for (;;) {								\
      EXP_MUL(a, g);							\
      sq = 0;								\
@@@ -184,11 -192,15 +195,15 @@@ exp_simple_exit:;							
  									\
    /* --- Do the precomputation --- */					\
  									\
+   EXP_FIX(g);								\
    EXP_SETSQR(g2, g);							\
+   EXP_FIX(g2);								\
    v = xmalloc(EXP_TABSZ * sizeof(EXP_TYPE));				\
    EXP_COPY(v[0], g);							\
-   for (i = 1; i < EXP_TABSZ; i++)					\
+   for (i = 1; i < EXP_TABSZ; i++) {					\
      EXP_SETMUL(v[i], v[i - 1], g2);					\
+     EXP_FIX(v[i]);							\
+   }									\
    EXP_DROP(g2);								\
  									\
    /* --- Skip top-end zero bits --- *					\
@@@ -286,17 -298,21 +301,21 @@@ exp_window_exit:;							
    j = 1;								\
    for (i = 0; i < n; i++) {						\
      EXP_COPY(v[j], f[n - 1 - i].base);					\
+     EXP_FIX(v[j]);							\
      j <<= 1;								\
    }									\
    k = n * EXP_WINSZ;							\
    jj = 1;								\
    for (; i < k; i++) {							\
      EXP_SETSQR(v[j], v[jj]);						\
+     EXP_FIX(v[j]);							\
      j <<= 1; jj <<= 1;							\
    }									\
    for (i = 1; i < vn; i <<= 1) {					\
-     for (j = 1; j < i; j++)						\
+     for (j = 1; j < i; j++) {						\
        EXP_SETMUL(v[j + i], v[j], v[i]);					\
+       EXP_FIX(v[j + i]);						\
+     }									\
    }									\
  									\
    /* --- Set up the bitscanners --- *					\
@@@ -381,7 -397,7 +400,7 @@@
  									\
  exp_simul_done:								\
    while (sq--) EXP_SQR(a);						\
-   for (i = 1; i < vn; i++)					\
+   for (i = 1; i < vn; i++)						\
      EXP_DROP(v[i]);							\
    xfree(v);								\
  } while (0)
diff --combined f-binpoly.c
index 00000000,509efc4f..02e683dd
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 --combined f-prime.c
index f4bf3a74,7215ec85..84549020
--- 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.
   *
@@@ -69,7 -78,8 +81,8 @@@ static void fdestroy(field *ff
  static mp *fin(field *ff, mp *d, mp *x)
  {
    fctx *f = (fctx *)ff;
-   return (mpmont_mul(&f->mm, d, x, f->mm.r2));
+   mp_div(0, &d, x, f->mm.m);
+   return (mpmont_mul(&f->mm, d, d, f->mm.r2));
  }
  
  static mp *fout(field *ff, mp *d, mp *x)
@@@ -78,6 -88,11 +91,11 @@@
    return (mpmont_reduce(&f->mm, d, x));
  }
  
+ static int fzerop(field *ff, mp *x)
+ {
+   return (!MP_LEN(x));
+ }
+ 
  static mp *fneg(field *ff, mp *d, mp *x)
  {
    fctx *f = (fctx *)ff;
@@@ -86,12 -101,24 +104,24 @@@
  
  static mp *fadd(field *ff, mp *d, mp *x, mp *y)
  {
-   return (mp_add(d, x, y));
+   fctx *f = (fctx *)ff;
+   d = mp_add(d, x, y);
+   if (d->f & MP_NEG)
+     d = mp_add(d, d, f->mm.m);
+   else if (MP_CMP(d, >, f->mm.m))
+     d = mp_sub(d, d, f->mm.m);
+   return (d);
  }
  
  static mp *fsub(field *ff, mp *d, mp *x, mp *y)
  {
-   return (mp_sub(d, x, y));
+   fctx *f = (fctx *)ff;
+   d = mp_sub(d, x, y);
+   if (d->f & MP_NEG)
+     d = mp_add(d, d, f->mm.m);
+   else if (MP_CMP(d, >, f->mm.m))
+     d = mp_sub(d, d, f->mm.m);
+   return (d);
  }
  
  static mp *fmul(field *ff, mp *d, mp *x, mp *y)
@@@ -122,26 -149,57 +152,57 @@@ static mp *freduce(field *ff, mp *d, m
    return (d);
  }
  
+ static mp *fsqrt(field *ff, mp *d, mp *x)
+ {
+   fctx *f = (fctx *)ff;
+   d = mpmont_reduce(&f->mm, d, x);
+   d = mp_modsqrt(d, d, f->mm.m);
+   if (!d)
+     return (d);
+   return (mpmont_mul(&f->mm, d, d, f->mm.r2));
+ }
+ 
  static mp *fdbl(field *ff, mp *d, mp *x)
  {
- /*   fctx *f = (fctx *)ff; */
-   return (mp_lsl(d, x, 1));
+   fctx *f = (fctx *)ff;
+   d = mp_lsl(d, x, 1);
+   if (MP_CMP(d, >, f->mm.m))
+     d = mp_sub(d, d, f->mm.m);
+   return (d);
  }
  
  static mp *ftpl(field *ff, mp *d, mp *x)
  {
- /*   fctx *f = (fctx *)ff; */
+   fctx *f = (fctx *)ff;
    MP_DEST(d, MP_LEN(x) + 1, x->f);
    MPX_UMULN(d->v, d->vl, x->v, x->vl, 3);
+   while (MP_CMP(d, >, f->mm.m))
+     d = mp_sub(d, d, f->mm.m);
    return (d);
  }
  
- static mp *fsqrt(field *ff, mp *d, mp *x)
+ static mp *fqdl(field *ff, mp *d, mp *x)
  {
    fctx *f = (fctx *)ff;
-   d = mpmont_reduce(&f->mm, d, x);
-   d = mp_modsqrt(d, d, f->mm.m);
-   return (mpmont_mul(&f->mm, d, d, f->mm.r2));
+   d = mp_lsl(d, x, 2);
+   while (MP_CMP(d, >, f->mm.m))
+     d = mp_sub(d, d, f->mm.m);
+   return (d);
+ }
+ 
+ static mp *fhlv(field *ff, mp *d, mp *x)
+ {
+   fctx *f = (fctx *)ff;
+   if (!MP_LEN(x)) {
+     MP_COPY(x);
+     MP_DROP(d);
+     return (x);
+   }
+   if (x->v[0] & 1) {
+     d = mp_add(d, x, f->mm.m);
+     x = d;
+   }
+   return (mp_lsr(d, x, 1));
  }
  
  /* --- Field operations table --- */
@@@ -149,8 -207,9 +210,9 @@@
  static field_ops fops = {
    fdestroy,
    fin, fout,
-   fneg, fadd, fsub, fmul, fsqr, finv, freduce,
-   fdbl, ftpl, fsqrt
+   fzerop, fneg, fadd, fsub, fmul, fsqr, finv, freduce, fsqrt,
+   0,
+   fdbl, ftpl, fqdl, fhlv
  };
  
  /* --- @field_prime@ --- *
diff --combined field.c
index c7e7559a,0f02fa7b..c310b938
--- 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 --combined field.h
index 5a706498,ea019c5d..dd674c92
--- 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.
   *
@@@ -70,6 -76,7 +79,7 @@@ typedef struct field_ops 
    mp *(*in)(field */*f*/, mp */*d*/, mp */*x*/);
    mp *(*out)(field */*f*/, mp */*d*/, mp */*x*/);
  
+   int (*zerop)(field */*f*/, mp */*x*/);
    mp *(*neg)(field */*f*/, mp */*d*/, mp */*x*/);
    mp *(*add)(field */*f*/, mp */*d*/, mp */*x*/, mp */*y*/);
    mp *(*sub)(field */*f*/, mp */*d*/, mp */*x*/, mp */*y*/);
@@@ -77,12 -84,18 +87,18 @@@
    mp *(*sqr)(field */*f*/, mp */*d*/, mp */*x*/);
    mp *(*inv)(field */*f*/, mp */*d*/, mp */*x*/);
    mp *(*reduce)(field */*f*/, mp */*d*/, mp */*x*/);
+   mp *(*sqrt)(field */*f*/, mp */*d*/, mp */*x*/);
+ 
+   /* --- Operations for binary fields only --- */
+ 
+   mp *(*quadsolve)(field */*f*/, mp */*d*/, mp */*x*/);
  
    /* --- Operations for prime fields only --- */
  
    mp *(*dbl)(field */*f*/, mp */*d*/, mp */*x*/);
    mp *(*tpl)(field */*f*/, mp */*d*/, mp */*x*/);
-   mp *(*sqrt)(field */*f*/, mp */*d*/, mp */*x*/);
+   mp *(*qdl)(field */*f*/, mp */*d*/, mp */*x*/);
+   mp *(*hlv)(field */*f*/, mp */*d*/, mp */*x*/);
  
  } field_ops;
  
@@@ -90,6 -103,8 +106,8 @@@
  
  #define F_IN(f, d, x)		(f)->ops->in((f), (d), (x))
  #define F_OUT(f, d, x)		(f)->ops->out((f), (d), (x))
+ 
+ #define F_ZEROP(f, x)		(f)->ops->zerop((f), (x))
  #define F_NEG(f, d, x)		(f)->ops->neg((f), (d), (x))
  #define F_ADD(f, d, x, y)	(f)->ops->add((f), (d), (x), (y))
  #define F_SUB(f, d, x, y)	(f)->ops->sub((f), (d), (x), (y))
@@@ -97,10 -112,14 +115,14 @@@
  #define F_SQR(f, d, x)		(f)->ops->sqr((f), (d), (x))
  #define F_INV(f, d, x)		(f)->ops->inv((f), (d), (x))
  #define F_REDUCE(f, d, x)	(f)->ops->reduce((f), (d), (x))
+ #define F_SQRT(f, d, x)		(f)->ops->sqrt((f), (d), (x))
+ 
+ #define F_QUADSOLVE(f, d, x)	(f)->ops->quadsolve((f), (d), (x))
  
  #define F_DBL(f, d, x)		(f)->ops->dbl((f), (d), (x))
  #define F_TPL(f, d, x)		(f)->ops->tpl((f), (d), (x))
- #define F_SQRT(f, d, x)		(f)->ops->sqrt((f), (d), (x))
+ #define F_QDL(f, d, x)		(f)->ops->qdl((f), (d), (x))
+ #define F_HLV(f, d, x)		(f)->ops->hlv((f), (d), (x))
  
  /*----- Helpful field operations ------------------------------------------*/
  
@@@ -144,20 -163,6 +166,6 @@@ extern field *field_prime(mp */*p*/)
  
  extern field *field_binpoly(mp */*p*/);
  
- /* --- @field_binsparsepoly@ --- *
-  *
-  * Arguments:	@unsigned deg@ = the degree of the field polynomial
-  *		@const unsigned *c@ = degrees of nonzero coefficients
-  *
-  * Returns:	A pointer to the field.
-  *
-  * Use:		Creates a field structure for a binary field taking advantage
-  *		of the sparseness of the given polynomial to improve
-  *		reduction performance.
-  */
- 
- extern field *field_binsparsepoly(unsigned /*deg*/, const unsigned */*c*/);
- 
  /*----- That's all, folks -------------------------------------------------*/
  
  #ifdef __cplusplus
diff --combined gf-arith.c
index 00000000,6838e44f..debfba14
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 --combined gf-gcd.c
index 00000000,64d61d4a..7c09d3ab
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 --combined gf.h
index 00000000,889cd9b1..ebf67f1b
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 --combined gfreduce-exp.h
index 00000000,0393145e..f826fc7a
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 --combined gfreduce.c
index 00000000,3969f117..819c2762
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 --combined gfreduce.h
index 00000000,2fc4c0ad..9840b5e1
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 --combined gfx-sqr.c
index b253a320,778f85a5..19ec5745
--- 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.
   *
@@@ -38,7 -41,7 +44,7 @@@
  /*----- Header files ------------------------------------------------------*/
  
  #include "mpx.h"
- /* #include "gfx.h" */
+ #include "gfx.h"
  #include "gfx-sqr-tab.h"
  
  /*----- Static variables --------------------------------------------------*/
@@@ -99,7 -102,7 +105,7 @@@ void gfx_sqr(mpw *dv, mpw *dvl, const m
  
      /* --- Output buffering --- */
  
-     if (bb > MPW_BITS) {
+     if (bb >= MPW_BITS) {
        *dv++ = MPW(aa);
        if (dv >= dvl)
  	return;
diff --combined gfx.h
index 778baeed,d525650d..18ac9a55
--- 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.
   *
@@@ -112,6 -115,19 +118,19 @@@ extern void gfx_mul(mpw */*dv*/, mpw */
  		    const mpw */*av*/, const mpw */*avl*/,
  		    const mpw */*bv*/, const mpw */*bvl*/);
  
+ /* --- @gfx_sqr@ --- *
+  *
+  * Arguments:	@mpw *dv, *dvl@ = destination vector base and limit
+  *		@const mpw *av, *avl@ = argument vector base and limit
+  *
+  * Returns:	---
+  *
+  * Use:		Performs squaring of binary polynomials.
+  */
+ 
+ extern void gfx_sqr(mpw */*dv*/, mpw */*dvl*/,
+ 		    const mpw */*av*/, const mpw */*avl*/);
+ 
  /* --- @gfx_div@ --- *
   *
   * Arguments:	@mpw *qv, *qvl@ = quotient vector base and limit
diff --combined mp-gcd.c
index a17a502f,f55d0aab..6135e545
--- 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.
   *
@@@ -74,12 -77,10 +80,10 @@@ void mp_gcd(mp **gcd, mp **xx, mp **yy
    mp *q = MP_NEW;
    unsigned f = 0;
  
-   enum {
-     f_swap = 1u,
-     f_aneg = 2u,
-     f_bneg = 4u,
-     f_ext = 8u
-   };
+ #define f_swap 1u
+ #define f_aneg 2u
+ #define f_bneg 4u
+ #define f_ext 8u
  
    /* --- Sort out some initial flags --- */
  
diff --combined mpbarrett-exp.h
index 5c34746a,e34ea498..dd02637a
--- 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.
   *
@@@ -61,6 -64,8 +67,8 @@@
    a = mpbarrett_reduce(mb, t, t);					\
  } while (0)
  
+ #define EXP_FIX(x)
+ 
  #define EXP_SETMUL(d, x, y) do {					\
    d = mp_mul(MP_NEW, x, y);						\
    d = mpbarrett_reduce(mb, d, d);					\
diff --combined mpmont-exp.h
index e732e6fa,0f82b9c4..5f2b31d2
--- 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.
   *
@@@ -61,6 -64,8 +67,8 @@@
    a = mpmont_reduce(mm, t, t);						\
  } while (0)
  
+ #define EXP_FIX(x)
+ 
  #define EXP_SETMUL(d, x, y) d = mpmont_mul(mm, MP_NEW, x, y)
  
  #define EXP_SETSQR(d, x) do {						\
diff --combined tests/gf
index 00000000,0c3987fa..bbb1514e
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 --combined tests/gfreduce
index 00000000,aec5318c..806ec284
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 --combined tests/gfx
index 17bdbf7e,866bd39b..1ef60816
--- 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) ---
  
@@@ -356,6 -356,11 +356,11 @@@ sqr 
      1500551145554555115501404405515454551004050045401511450544004455;
    93c9240145720297b96a404941792a21
      4105504104100001101115040004411545411444100010411001154104440401;
+ 
+   # --- Regression ---
+ 
+   01f081e69f45d3254530766ab98d55fa612c7bb27ea31bc2621d894be9c0b196b3
+     0155004001541441551011510504111011050015141444454140511111554414010450154545041554440501455004140401514041104554415000450141144505;
  }
  
  div {