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

/* -*-c-*-
 *
 * $Id: ec-bin.c,v 1.6 2004/04/01 12:50:09 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.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);
    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;
    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.
 *
 * 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, b);
  cc->bb = F_SQRT(f, cc->bb, cc->bb);
  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;
  int i, n = argc == 1 ? 1 : atoi(argv[1]);

  printf("ec-bin: ");
  fflush(stdout);
  a = MP(1);
  b = MP(0x021a5c2c8ee9feb5c4b9a753b7b476b7fd6422ef1f3dd674761fa99d6ac27c8a9a197b272822f6cd57a55aa4f50ae317b13545f);
  p = MP(0x2000000000000000000000000000000000000000000000000000000000000000000000000000000008000000000000000000001);
  r =
  MP(661055968790248598951915308032771039828404682964281219284648798304157774827374805208143723762179110965979867288366567526770);

  f = field_binpoly(p);
  c = ec_binproj(f, a, b);
  
  g.x = MP(0x15d4860d088ddb3496b0c6064756260441cde4af1771d4db01ffe5b34e59703dc255a868a1180515603aeab60794e54bb7996a7);
  g.y = MP(0x061b1cfab6be5f32bbfa78324ed106a7636b9c5a7bd198d0158aa4f5488d08f38514f1fdf4b4f40d2181b3681c364ba0273c706);

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

#endif

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