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

/* -*-c-*-
 *
 * $Id$
 *
 * 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.
 */

/*----- Header files ------------------------------------------------------*/

#include <mLib/sub.h>

#include "ec.h"
#include "ec-guts.h"

/*----- 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_bin *cc = (ecctx_bin *)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_bin *cc = (ecctx_bin *)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_bin *cc = CREATE(ecctx_bin);
  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_bin *cc = CREATE(ecctx_bin);
  int i;
  mp *c, *d;

  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);

  c = MP_COPY(cc->c.b);
  for (i = 0; i < f->nbits - 2; i++)
    c = F_SQR(f, c, c);
  d = F_SQR(f, MP_NEW, c); d = F_SQR(f, d, d);
  if (!MP_EQ(d, cc->c.b)) {
    MP_DROP(c);
    MP_DROP(d);
    MP_DROP(cc->c.a);
    MP_DROP(cc->c.b);
    DESTROY(cc);
    return (0);
  }
  cc->bb = c;
  MP_DROP(d);
  return (&cc->c);
}

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

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