[u/mdw/catacomb] / ec-bin.c

/* -*-c-*-
 *
 * $Id: ec-bin.c,v 1.8 2004/04/03 03:32:05 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.8  2004/04/03 03:32:05  mdw
 * General robustification.
 *
 * Revision 1.7  2004/04/01 21:28:41  mdw
 * Normal basis support (translates to poly basis internally).  Rewrite
 * EC and prime group table generators in awk, so that they can reuse data
 * for repeated constants.
 *
 * Revision 1.6  2004/04/01 12:50:09  mdw
 * Add cyclic group abstraction, with test code.  Separate off exponentation
 * functions for better static linking.  Fix a buttload of bugs on the way.
 * Generally ensure that negative exponents do inversion correctly.  Add
 * table of standard prime-field subgroups.  (Binary field subgroups are
 * currently unimplemented but easy to add if anyone ever finds a good one.)
 *
 * Revision 1.5  2004/03/27 17:54:11  mdw
 * Standard curves and curve checking.
 *
 * Revision 1.4  2004/03/23 15:19:32  mdw
 * Test elliptic curves more thoroughly.
 *
 * Revision 1.3  2004/03/22 02:19:09  mdw
 * Rationalise the sliding-window threshold.  Drop guarantee that right
 * arguments to EC @add@ are canonical, and fix up projective implementations
 * to cope.
 *
 * 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 *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)
{
  field *f = c->f;
  mp *y, *u, *v;
  
  if (F_ZEROP(f, x))
    y = F_SQRT(f, MP_NEW, c->b);
  else {
    u = F_SQR(f, MP_NEW, x);		/* %$x^2$% */
    y = F_MUL(f, MP_NEW, u, c->a);	/* %$a x^2$% */
    y = F_ADD(f, y, y, c->b);		/* %$a x^2 + b$% */
    v = F_MUL(f, MP_NEW, u, x);		/* %$x^3$% */
    y = F_ADD(f, y, y, v);		/* %$A = x^3 + a x^2 + b$% */
    if (!F_ZEROP(f, y)) {
      u = F_INV(f, u, u);		/* %$x^{-2}$% */
      v = F_MUL(f, v, u, y);	    /* %$B = A x^{-2} = x + a + b x^{-2}$% */
      y = F_QUADSOLVE(f, y, v);		/* %$z^2 + z = B$% */
      if (y) y = F_MUL(f, y, y, x);	/* %$y = z x$% */
    }
    MP_DROP(u);
    MP_DROP(v);
  }
  if (!y) return (0);
  EC_DESTROY(d);
  d->x = MP_COPY(x);
  d->y = y;
  d->z = MP_COPY(f->one);
  return (d);
}

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;
    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, c->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);
  }
  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;
    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, c->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, c->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;
    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, c->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)
{
  field *f = c->f;
  int rc;
  mp *u, *v;

  if (EC_ATINF(p)) return (0);
  v = F_SQR(f, MP_NEW, p->x);
  u = F_MUL(f, MP_NEW, v, p->x);
  v = F_MUL(f, v, v, c->a);
  u = F_ADD(f, u, u, v);
  u = F_ADD(f, u, u, c->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) ? 0 : -1;
  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->c.a);
  MP_DROP(cc->c.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, or null.
 *
 * 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->c.a = F_IN(f, MP_NEW, a);
  cc->c.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->c.a = F_IN(f, MP_NEW, a);
  cc->c.b = F_IN(f, MP_NEW, b);
  cc->bb = F_SQRT(f, MP_NEW, cc->c.b);
  if (cc->bb)
    cc->bb = F_SQRT(f, cc->bb, cc->bb);
  if (!cc->bb) {
    MP_DROP(cc->c.a);
    MP_DROP(cc->c.b);
    DESTROY(cc);
    return (0);
  }
  return (&cc->c);
}

static const ec_ops ec_binops = {
  ecdestroy, ec_stdsamep, ec_idin, ec_idout, ec_idfix,
  ecfind, ecneg, ecadd, ec_stdsub, ecdbl, eccheck
};

static const ec_ops ec_binprojops = {
  ecdestroy, ec_stdsamep, ec_projin, ec_projout, ec_projfix,
  ecfind, 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, *beta;
  int i, n = argc == 1 ? 1 : atoi(argv[1]);

  printf("ec-bin: ");
  fflush(stdout);
  a = MP(0x7ffffffffffffffffffffffffffffffffffffffff);
  b = MP(0x6645f3cacf1638e139c6cd13ef61734fbc9e3d9fb);
  p = MP(0x800000000000000000000000000000000000000c9);
  beta = MP(0x715169c109c612e390d347c748342bcd3b02a0bef);
  r = MP(0x040000000000000000000292fe77e70c12a4234c32);

  f = field_binnorm(p, beta);
  c = ec_binproj(f, a, b);
  g.x = MP(0x0311103c17167564ace77ccb09c681f886ba54ee8);
  g.y = MP(0x333ac13c6447f2e67613bf7009daf98c87bb50c7f);

  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);
      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); MP_DROP(beta);
  assert(!mparena_count(&mparena_global));
  printf("ok\n");
  return (0);
}

#endif

/*----- That's all, folks -------------------------------------------------*/
Commit	Line	Data
ceb3f0c0	1	/* --c--
ceb3f0c0	2	*
02d7884d	3	* $Id: ec-bin.c,v 1.8 2004/04/03 03:32:05 mdw Exp $
ceb3f0c0	4	*
	5	* Arithmetic for elliptic curves over binary fields
	6	*
	7	* (c) 2004 Straylight/Edgeware
	8	*/
	9
	10	/----- Licensing notice --------------------------------------------------
	11	*
	12	* This file is part of Catacomb.
	13	*
	14	* Catacomb is free software; you can redistribute it and/or modify
	15	* it under the terms of the GNU Library General Public License as
	16	* published by the Free Software Foundation; either version 2 of the
	17	* License, or (at your option) any later version.
	18	*
	19	* Catacomb is distributed in the hope that it will be useful,
	20	* but WITHOUT ANY WARRANTY; without even the implied warranty of
	21	* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
	22	* GNU Library General Public License for more details.
	23	*
	24	* You should have received a copy of the GNU Library General Public
	25	* License along with Catacomb; if not, write to the Free
	26	* Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
	27	* MA 02111-1307, USA.
	28	*/
	29
	30	/----- Revision history --------------------------------------------------
	31	*
	32	* $Log: ec-bin.c,v $
02d7884d	33	* Revision 1.8 2004/04/03 03:32:05 mdw
	34	* General robustification.
	35	*
4edc47b8	36	* Revision 1.7 2004/04/01 21:28:41 mdw
	37	* Normal basis support (translates to poly basis internally). Rewrite
	38	* EC and prime group table generators in awk, so that they can reuse data
	39	* for repeated constants.
	40	*
34e4f738	41	* Revision 1.6 2004/04/01 12:50:09 mdw
	42	* Add cyclic group abstraction, with test code. Separate off exponentation
	43	* functions for better static linking. Fix a buttload of bugs on the way.
	44	* Generally ensure that negative exponents do inversion correctly. Add
	45	* table of standard prime-field subgroups. (Binary field subgroups are
	46	* currently unimplemented but easy to add if anyone ever finds a good one.)
	47	*
432c4e18	48	* Revision 1.5 2004/03/27 17:54:11 mdw
	49	* Standard curves and curve checking.
	50	*
bc985cef	51	* Revision 1.4 2004/03/23 15:19:32 mdw
	52	* Test elliptic curves more thoroughly.
	53	*
391faf42	54	* Revision 1.3 2004/03/22 02:19:09 mdw
	55	* Rationalise the sliding-window threshold. Drop guarantee that right
	56	* arguments to EC @add@ are canonical, and fix up projective implementations
	57	* to cope.
	58	*
c3caa2fa	59	* Revision 1.2 2004/03/21 22:52:06 mdw
	60	* Merge and close elliptic curve branch.
	61	*
ceb3f0c0	62	* Revision 1.1.2.1 2004/03/21 22:39:46 mdw
	63	* Elliptic curves on binary fields work.
	64	*
	65	*/
	66
	67	/----- Header files ------------------------------------------------------/
	68
	69	#include <mLib/sub.h>
	70
	71	#include "ec.h"
	72
	73	/----- Data structures ---------------------------------------------------/
	74
	75	typedef struct ecctx {
	76	ec_curve c;
ceb3f0c0	77	mp *bb;
	78	} ecctx;
	79
	80	/----- Main code ---------------------------------------------------------/
	81
	82	static const ec_ops ec_binops, ec_binprojops;
	83
	84	static ec ecneg(ec_curve c, ec d, const ec p)
	85	{
	86	EC_COPY(d, p);
	87	if (d->x)
	88	d->y = F_ADD(c->f, d->y, d->y, d->x);
	89	return (d);
	90	}
	91
	92	static ec ecprojneg(ec_curve c, ec d, const ec p)
	93	{
	94	EC_COPY(d, p);
	95	if (d->x) {
	96	mp *t = F_MUL(c->f, MP_NEW, d->x, d->z);
	97	d->y = F_ADD(c->f, d->y, d->y, t);
	98	MP_DROP(t);
	99	}
	100	return (d);
	101	}
	102
	103	static ec ecfind(ec_curve c, ec d, mp x)
	104	{
bc985cef	105	field *f = c->f;
bc985cef	106	mp y, u, *v;
	107
	108	if (F_ZEROP(f, x))
432c4e18	109	y = F_SQRT(f, MP_NEW, c->b);
bc985cef	110	else {
bc985cef	111	u = F_SQR(f, MP_NEW, x); /* %$x^2$% */
432c4e18	112	y = F_MUL(f, MP_NEW, u, c->a); /* %$a x^2$% */
432c4e18	113	y = F_ADD(f, y, y, c->b); /* %$a x^2 + b$% */
bc985cef	114	v = F_MUL(f, MP_NEW, u, x); /* %$x^3$% */
	115	y = F_ADD(f, y, y, v); /* %$A = x^3 + a x^2 + b$% */
	116	if (!F_ZEROP(f, y)) {
	117	u = F_INV(f, u, u); /* %$x^{-2}$% */
	118	v = F_MUL(f, v, u, y); /* %$B = A x^{-2} = x + a + b x^{-2}$% */
	119	y = F_QUADSOLVE(f, y, v); /* %$z^2 + z = B$% */
	120	if (y) y = F_MUL(f, y, y, x); /* %$y = z x$% */
	121	}
	122	MP_DROP(u);
	123	MP_DROP(v);
	124	}
	125	if (!y) return (0);
	126	EC_DESTROY(d);
	127	d->x = MP_COPY(x);
	128	d->y = y;
	129	d->z = MP_COPY(f->one);
	130	return (d);
ceb3f0c0	131	}
	132
	133	static ec ecdbl(ec_curve c, ec d, const ec a)
	134	{
	135	if (EC_ATINF(a) \|\| F_ZEROP(c->f, a->x))
	136	EC_SETINF(d);
	137	else {
	138	field *f = c->f;
ceb3f0c0	139	mp *lambda;
	140	mp dx, dy;
	141
	142	dx = F_INV(f, MP_NEW, a->x); /* %$x^{-1}$% */
	143	dy = F_MUL(f, MP_NEW, dx, a->y); /* %$y/x$% */
	144	lambda = F_ADD(f, dy, dy, a->x); /* %$\lambda = x + y/x$% */
	145
	146	dx = F_SQR(f, dx, lambda); /* %$\lambda^2$% */
	147	dx = F_ADD(f, dx, dx, lambda); /* %$\lambda^2 + \lambda$% */
432c4e18	148	dx = F_ADD(f, dx, dx, c->a); /* %$x' = a + \lambda^2 + \lambda$% */
ceb3f0c0	149
	150	dy = F_ADD(f, MP_NEW, a->x, dx); /* %$ x + x' $% */
	151	dy = F_MUL(f, dy, dy, lambda); /* %$ (x + x') \lambda$% */
	152	dy = F_ADD(f, dy, dy, a->y); /* %$ (x + x') \lambda + y$% */
	153	dy = F_ADD(f, dy, dy, dx); /* %$ y' = (x + x') \lambda + y + x'$% */
	154
	155	EC_DESTROY(d);
	156	d->x = dx;
	157	d->y = dy;
	158	d->z = 0;
	159	MP_DROP(lambda);
	160	}
	161	return (d);
	162	}
	163
	164	static ec ecprojdbl(ec_curve c, ec d, const ec a)
	165	{
	166	if (EC_ATINF(a) \|\| F_ZEROP(c->f, a->x))
	167	EC_SETINF(d);
	168	else {
	169	field *f = c->f;
	170	ecctx cc = (ecctx )c;
	171	mp dx, dy, dz, u, *v;
	172
	173	dy = F_SQR(f, MP_NEW, a->z); /* %$z^2$% */
	174	dx = F_MUL(f, MP_NEW, dy, cc->bb); /* %$c z^2$% */
	175	dx = F_ADD(f, dx, dx, a->x); /* %$x + c z^2$% */
	176	dz = F_SQR(f, MP_NEW, dx); /* %$(x + c z^2)^2$% */
	177	dx = F_SQR(f, dx, dz); /* %$x' = (x + c z^2)^4$% */
	178
	179	dz = F_MUL(f, dz, dy, a->x); /* %$z' = x z^2$% */
	180
	181	dy = F_SQR(f, dy, a->x); /* %$x^2$% */
	182	u = F_MUL(f, MP_NEW, a->y, a->z); /* %$y z$% */
	183	u = F_ADD(f, u, u, dz); /* %$z' + y z$% */
	184	u = F_ADD(f, u, u, dy); /* %$u = z' + x^2 + y z$% */
	185
	186	v = F_SQR(f, MP_NEW, dy); /* %$x^4$% */
	187	dy = F_MUL(f, dy, v, dz); /* %$x^4 z'$% */
	188	v = F_MUL(f, v, u, dx); /* %$u x'$% */
	189	dy = F_ADD(f, dy, dy, v); /* %$y' = x^4 z' + u x'$% */
	190
	191	EC_DESTROY(d);
	192	d->x = dx;
	193	d->y = dy;
	194	d->z = dz;
	195	MP_DROP(u);
	196	MP_DROP(v);
ceb3f0c0	197	}
	198	return (d);
	199	}
	200
	201	static ec ecadd(ec_curve c, ec d, const ec a, const ec *b)
	202	{
	203	if (a == b)
	204	ecdbl(c, d, a);
	205	else if (EC_ATINF(a))
	206	EC_COPY(d, b);
	207	else if (EC_ATINF(b))
	208	EC_COPY(d, a);
	209	else {
	210	field *f = c->f;
ceb3f0c0	211	mp *lambda;
	212	mp dx, dy;
	213
	214	if (!MP_EQ(a->x, b->x)) {
	215	dx = F_ADD(f, MP_NEW, a->x, b->x); /* %$x_0 + x_1$% */
	216	dy = F_INV(f, MP_NEW, dx); /* %$(x_0 + x_1)^{-1}$% */
	217	dx = F_ADD(f, dx, a->y, b->y); /* %$y_0 + y_1$% */
	218	lambda = F_MUL(f, MP_NEW, dy, dx);
	219	/* %$\lambda = (y_0 + y_1)/(x_0 + x_1)$% */
	220
	221	dx = F_SQR(f, dx, lambda); /* %$\lambda^2$% */
	222	dx = F_ADD(f, dx, dx, lambda); /* %$\lambda^2 + \lambda$% */
432c4e18	223	dx = F_ADD(f, dx, dx, c->a); /* %$a + \lambda^2 + \lambda$% */
ceb3f0c0	224	dx = F_ADD(f, dx, dx, a->x); /* %$a + \lambda^2 + \lambda + x_0$% */
	225	dx = F_ADD(f, dx, dx, b->x);
	226	/* %$x' = a + \lambda^2 + \lambda + x_0 + x_1$% */
	227	} else if (!MP_EQ(a->y, b->y) \|\| F_ZEROP(f, a->x)) {
	228	EC_SETINF(d);
	229	return (d);
	230	} else {
	231	dx = F_INV(f, MP_NEW, a->x); /* %$x^{-1}$% */
	232	dy = F_MUL(f, MP_NEW, dx, a->y); /* %$y/x$% */
	233	lambda = F_ADD(f, dy, dy, a->x); /* %$\lambda = x + y/x$% */
	234
	235	dx = F_SQR(f, dx, lambda); /* %$\lambda^2$% */
	236	dx = F_ADD(f, dx, dx, lambda); /* %$\lambda^2 + \lambda$% */
432c4e18	237	dx = F_ADD(f, dx, dx, c->a); /* %$x' = a + \lambda^2 + \lambda$% */
ceb3f0c0	238	dy = MP_NEW;
	239	}
	240
	241	dy = F_ADD(f, dy, a->x, dx); /* %$ x + x' $% */
	242	dy = F_MUL(f, dy, dy, lambda); /* %$ (x + x') \lambda$% */
	243	dy = F_ADD(f, dy, dy, a->y); /* %$ (x + x') \lambda + y$% */
	244	dy = F_ADD(f, dy, dy, dx); /* %$ y' = (x + x') \lambda + y + x'$% */
	245
	246	EC_DESTROY(d);
	247	d->x = dx;
	248	d->y = dy;
	249	d->z = 0;
	250	MP_DROP(lambda);
	251	}
	252	return (d);
	253	}
	254
	255	static ec ecprojadd(ec_curve c, ec d, const ec a, const ec *b)
	256	{
	257	if (a == b)
	258	c->ops->dbl(c, d, a);
	259	else if (EC_ATINF(a))
	260	EC_COPY(d, b);
	261	else if (EC_ATINF(b))
	262	EC_COPY(d, a);
	263	else {
	264	field *f = c->f;
ceb3f0c0	265	mp dx, dy, dz, u, uu, v, t, s, ss, r, w, l;
	266
	267	dz = F_SQR(f, MP_NEW, b->z); /* %$z_1^2$% */
	268	u = F_MUL(f, MP_NEW, dz, a->x); /* %$u_0 = x_0 z_1^2$% */
	269	t = F_MUL(f, MP_NEW, dz, b->z); /* %$z_1^3$% */
	270	s = F_MUL(f, MP_NEW, t, a->y); /* %$s_0 = y_0 z_1^3$% */
	271
	272	dz = F_SQR(f, dz, a->z); /* %$z_0^2$% */
	273	uu = F_MUL(f, MP_NEW, dz, b->x); /* %$u_1 = x_1 z_0^2$% */
	274	t = F_MUL(f, t, dz, a->z); /* %$z_0^3$% */
	275	ss = F_MUL(f, MP_NEW, t, b->y); /* %$s_1 = y_1 z_0^3$% */
	276
	277	w = F_ADD(f, u, u, uu); /* %$r = u_0 + u_1$% */
	278	r = F_ADD(f, s, s, ss); /* %$w = s_0 + s_1$% */
	279	if (F_ZEROP(f, w)) {
	280	MP_DROP(w);
	281	MP_DROP(uu);
	282	MP_DROP(ss);
	283	MP_DROP(t);
	284	MP_DROP(dz);
	285	if (F_ZEROP(f, r)) {
	286	MP_DROP(r);
	287	return (c->ops->dbl(c, d, a));
	288	} else {
	289	MP_DROP(r);
	290	EC_SETINF(d);
	291	return (d);
	292	}
	293	}
	294
	295	l = F_MUL(f, t, a->z, w); /* %$l = z_0 w$% */
	296
	297	dz = F_MUL(f, dz, l, b->z); /* %$z' = l z_1$% */
	298
	299	ss = F_MUL(f, ss, r, b->x); /* %$r x_1$% */
	300	t = F_MUL(f, uu, l, b->y); /* %$l y_1$% */
	301	v = F_ADD(f, ss, ss, t); /* %$v = r x_1 + l y_1$% */
	302
	303	t = F_ADD(f, t, r, dz); /* %$t = r + z'$% */
	304
	305	uu = F_SQR(f, MP_NEW, dz); /* %$z'^2$% */
432c4e18	306	dx = F_MUL(f, MP_NEW, uu, c->a); /* %$a z'^2$% */
ceb3f0c0	307	uu = F_MUL(f, uu, t, r); /* %$t r$% */
	308	dx = F_ADD(f, dx, dx, uu); /* %$a z'^2 + t r$% */
	309	r = F_SQR(f, r, w); /* %$w^2$% */
	310	uu = F_MUL(f, uu, r, w); /* %$w^3$% */
	311	dx = F_ADD(f, dx, dx, uu); /* %$x' = a z'^2 + t r + w^3$% */
	312
	313	r = F_SQR(f, r, l); /* %$l^2$% */
	314	dy = F_MUL(f, uu, v, r); /* %$v l^2$% */
	315	l = F_MUL(f, l, t, dx); /* %$t x'$% */
	316	dy = F_ADD(f, dy, dy, l); /* %$y' = t x' + v l^2$% */
	317
	318	EC_DESTROY(d);
	319	d->x = dx;
	320	d->y = dy;
	321	d->z = dz;
	322	MP_DROP(l);
	323	MP_DROP(r);
	324	MP_DROP(w);
	325	MP_DROP(t);
	326	MP_DROP(v);
	327	}
	328	return (d);
	329	}
	330
	331	static int eccheck(ec_curve c, const ec p)
	332	{
ceb3f0c0	333	field *f = c->f;
	334	int rc;
	335	mp u, v;
	336
34e4f738	337	if (EC_ATINF(p)) return (0);
ceb3f0c0	338	v = F_SQR(f, MP_NEW, p->x);
ceb3f0c0	339	u = F_MUL(f, MP_NEW, v, p->x);
432c4e18	340	v = F_MUL(f, v, v, c->a);
ceb3f0c0	341	u = F_ADD(f, u, u, v);
432c4e18	342	u = F_ADD(f, u, u, c->b);
ceb3f0c0	343	v = F_MUL(f, v, p->x, p->y);
	344	u = F_ADD(f, u, u, v);
	345	v = F_SQR(f, v, p->y);
	346	u = F_ADD(f, u, u, v);
bc985cef	347	rc = F_ZEROP(f, u) ? 0 : -1;
ceb3f0c0	348	mp_drop(u);
	349	mp_drop(v);
	350	return (rc);
	351	}
	352
	353	static int ecprojcheck(ec_curve c, const ec p)
	354	{
	355	ec t = EC_INIT;
	356	int rc;
	357
	358	c->ops->fix(c, &t, p);
	359	rc = eccheck(c, &t);
	360	EC_DESTROY(&t);
	361	return (rc);
	362	}
	363
	364	static void ecdestroy(ec_curve *c)
	365	{
	366	ecctx cc = (ecctx )c;
432c4e18	367	MP_DROP(cc->c.a);
432c4e18	368	MP_DROP(cc->c.b);
ceb3f0c0	369	if (cc->bb) MP_DROP(cc->bb);
	370	DESTROY(cc);
	371	}
	372
	373	/* --- @ec_bin@, @ec_binproj@ --- *
	374	*
	375	* Arguments: @field *f@ = the underlying field for this elliptic curve
	376	* @mp a, b@ = the coefficients for this curve
	377	*
02d7884d	378	* Returns: A pointer to the curve, or null.
ceb3f0c0	379	*
	380	* Use: Creates a curve structure for an elliptic curve defined over
	381	* a binary field. The @binproj@ variant uses projective
	382	* coordinates, which can be a win.
	383	*/
	384
	385	ec_curve ec_bin(field f, mp a, mp b)
	386	{
	387	ecctx *cc = CREATE(ecctx);
	388	cc->c.ops = &ec_binops;
	389	cc->c.f = f;
432c4e18	390	cc->c.a = F_IN(f, MP_NEW, a);
432c4e18	391	cc->c.b = F_IN(f, MP_NEW, b);
ceb3f0c0	392	cc->bb = 0;
	393	return (&cc->c);
	394	}
	395
	396	ec_curve ec_binproj(field f, mp a, mp b)
	397	{
	398	ecctx *cc = CREATE(ecctx);
	399	cc->c.ops = &ec_binprojops;
	400	cc->c.f = f;
432c4e18	401	cc->c.a = F_IN(f, MP_NEW, a);
432c4e18	402	cc->c.b = F_IN(f, MP_NEW, b);
4edc47b8	403	cc->bb = F_SQRT(f, MP_NEW, cc->c.b);
02d7884d	404	if (cc->bb)
	405	cc->bb = F_SQRT(f, cc->bb, cc->bb);
	406	if (!cc->bb) {
	407	MP_DROP(cc->c.a);
	408	MP_DROP(cc->c.b);
	409	DESTROY(cc);
	410	return (0);
	411	}
ceb3f0c0	412	return (&cc->c);
	413	}
	414
	415	static const ec_ops ec_binops = {
34e4f738	416	ecdestroy, ec_stdsamep, ec_idin, ec_idout, ec_idfix,
bc985cef	417	ecfind, ecneg, ecadd, ec_stdsub, ecdbl, eccheck
ceb3f0c0	418	};
	419
	420	static const ec_ops ec_binprojops = {
34e4f738	421	ecdestroy, ec_stdsamep, ec_projin, ec_projout, ec_projfix,
bc985cef	422	ecfind, ecprojneg, ecprojadd, ec_stdsub, ecprojdbl, ecprojcheck
ceb3f0c0	423	};
	424
	425	/----- Test rig ----------------------------------------------------------/
	426
	427	#ifdef TEST_RIG
	428
	429	#define MP(x) mp_readstring(MP_NEW, #x, 0, 0)
	430
	431	int main(int argc, char *argv[])
	432	{
	433	field *f;
	434	ec_curve *c;
	435	ec g = EC_INIT, d = EC_INIT;
4edc47b8	436	mp p, a, b, r, *beta;
ceb3f0c0	437	int i, n = argc == 1 ? 1 : atoi(argv[1]);
	438
	439	printf("ec-bin: ");
	440	fflush(stdout);
4edc47b8	441	a = MP(0x7ffffffffffffffffffffffffffffffffffffffff);
	442	b = MP(0x6645f3cacf1638e139c6cd13ef61734fbc9e3d9fb);
	443	p = MP(0x800000000000000000000000000000000000000c9);
	444	beta = MP(0x715169c109c612e390d347c748342bcd3b02a0bef);
	445	r = MP(0x040000000000000000000292fe77e70c12a4234c32);
ceb3f0c0	446
4edc47b8	447	f = field_binnorm(p, beta);
ceb3f0c0	448	c = ec_binproj(f, a, b);
4edc47b8	449	g.x = MP(0x0311103c17167564ace77ccb09c681f886ba54ee8);
4edc47b8	450	g.y = MP(0x333ac13c6447f2e67613bf7009daf98c87bb50c7f);
ceb3f0c0	451
	452	for (i = 0; i < n; i++) {
	453	ec_mul(c, &d, &g, r);
	454	if (EC_ATINF(&d)) {
	455	fprintf(stderr, "zero too early\n");
	456	return (1);
	457	}
	458	ec_add(c, &d, &d, &g);
	459	if (!EC_ATINF(&d)) {
	460	fprintf(stderr, "didn't reach zero\n");
	461	MP_EPRINTX("d.x", d.x);
	462	MP_EPRINTX("d.y", d.y);
ceb3f0c0	463	return (1);
	464	}
	465	ec_destroy(&d);
	466	}
	467
	468	ec_destroy(&g);
	469	ec_destroycurve(c);
	470	F_DESTROY(f);
4edc47b8	471	MP_DROP(p); MP_DROP(a); MP_DROP(b); MP_DROP(r); MP_DROP(beta);
ceb3f0c0	472	assert(!mparena_count(&mparena_global));
	473	printf("ok\n");
	474	return (0);
	475	}
	476
	477	#endif
	478
	479	/----- That's all, folks -------------------------------------------------/