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

/* -*-c-*-
 *
 * $Id: ec-bin.c,v 1.4 2004/03/23 15:19:32 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.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 *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)
{
  field *f = c->f;
  ecctx *cc = (ecctx *)c;
  mp *y, *u, *v;
  
  if (F_ZEROP(f, x))
    y = F_SQRT(f, MP_NEW, cc->b);
  else {
    u = F_SQR(f, MP_NEW, x);		/* %$x^2$% */
    y = F_MUL(f, MP_NEW, u, cc->a);	/* %$a x^2$% */
    y = F_ADD(f, y, y, cc->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;
    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) ? 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->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,
  ecfind, ecneg, ecadd, ec_stdsub, ecdbl, eccheck
};

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