From 8823192f6413bed15cfa884ed3a3cbbb97885657 Mon Sep 17 00:00:00 2001 From: mdw Date: Sat, 20 Mar 2004 00:13:31 +0000 Subject: [PATCH] Projective coordinates for prime curves --- calc/ecp.cal | 74 ++++++++++++++- ec-exp.h | 6 +- ec-prime.c | 300 +++++++++++++++++++++++++++++++++++++++++++++++++++++------ ec.c | 80 ++++++++++++++-- ec.h | 61 ++++++------ exp.h | 24 ++++- f-prime.c | 40 +++++++- field.h | 12 ++- mpmont-exp.h | 7 +- 9 files changed, 524 insertions(+), 80 deletions(-) diff --git a/calc/ecp.cal b/calc/ecp.cal index 82600f5..43ac1b3 100644 --- a/calc/ecp.cal +++ b/calc/ecp.cal @@ -1,6 +1,6 @@ /* -*-apcalc-*- * - * $Id: ecp.cal,v 1.1.4.1 2003/06/10 13:43:53 mdw Exp $ + * $Id: ecp.cal,v 1.1.4.2 2004/03/20 00:13:31 mdw Exp $ * * Testbed for elliptic curve arithmetic over prime fields * @@ -30,6 +30,9 @@ /*----- Revision history --------------------------------------------------* * * $Log: ecp.cal,v $ + * Revision 1.1.4.2 2004/03/20 00:13:31 mdw + * Projective coordinates for prime curves + * * Revision 1.1.4.1 2003/06/10 13:43:53 mdw * Simple (non-projective) curves over prime fields now seem to work. * @@ -42,6 +45,7 @@ obj ecp_curve { a, b, p }; obj ecp_pt { x, y, e }; +obj ecpp_pt { x, y, z, e }; /*----- Main code ---------------------------------------------------------*/ @@ -63,6 +67,72 @@ define ecp_pt(x, y, e) return (p); } +define ecpp_pt(p) +{ + local obj ecpp_pt pp; + if (istype(p, 1)) + return (0); + pp.x = p.x; + pp.y = p.y; + pp.z = 1; + pp.e = p.e; + return (pp); +} + +define ecpp_fix(pp) +{ + local obj ecp_pt p; + local e, zi, z2, z3; + if (istype(pp, 1) || pp.z == 0) + return (0); + e = pp.e; + zi = minv(pp.z, e.p); + z2 = zi * zi; + z3 = zi * z2; + p.x = pp.x * z2 % e.p; + p.y = pp.y * z3 % e.p; + p.e = e; + return (p); +} + +define ecpp_dbl(a) +{ + local m, s, t, y2; + local e; + local obj ecpp_pt d; + if (istype(a, 1) || a.y == 0) + return (0); + e = a.e; + if (e.a % e.p == e.p - 3) { + m = a.z^3 % e.p; + m = 3 * (a.x + t4) * (a.x - t4) % e.p; + } else { + m = (3 * a.x^2 - e.a * a.z^4) % e.p; + } + d.z = 2 * a.y * a.z % e.p; + y2 = a.y^2 % e.p; + s = 4 * a.x * a.y % e.p; + d.x = (m^2 - 2 * s) % e.p; + d.y = (m * (s - d.x) - y * y2^2) % e.p; + d.e = e; + return (d); +} + +define ecpp_add(a, b) +{ + if (a == 0) + d = b; + else if (b == 0) + d = a; + else if (!istype(a, b)) + quit "bad type arguments to ecp_pt_add"; + else if (a.e != b.e) + quit "points from different curves in ecp_pt_add"; + else { + e = a.e; + +} + define ecp_pt_print(a) { print "(" : a.x : ", " : a.y : ")" :; @@ -103,6 +173,8 @@ define ecp_pt_dbl(a) { local e, alpha; local obj ecp_pt d; + if (istype(a, 1)) + return (0); e = a.e; alpha = (3 * a.x^2 + e.a) * minv(2 * a.y, e.p) % e.p; d.x = (alpha^2 - 2 * a.x) % e.p; diff --git a/ec-exp.h b/ec-exp.h index 9ba7bc5..bb4e08a 100644 --- a/ec-exp.h +++ b/ec-exp.h @@ -1,6 +1,6 @@ /* -*-c-*- * - * $Id: ec-exp.h,v 1.2 2003/05/15 23:25:59 mdw Exp $ + * $Id: ec-exp.h,v 1.2.4.1 2004/03/20 00:13:31 mdw Exp $ * * Exponentiation operations for elliptic curves * @@ -30,6 +30,9 @@ /*----- Revision history --------------------------------------------------* * * $Log: ec-exp.h,v $ + * Revision 1.2.4.1 2004/03/20 00:13:31 mdw + * Projective coordinates for prime curves + * * Revision 1.2 2003/05/15 23:25:59 mdw * Make elliptic curve stuff build. * @@ -58,6 +61,7 @@ #define EXP_MUL(a, x) EC_ADD(c, &(a), &(a), &(x)) #define EXP_SQR(a) EC_DBL(c, &(a), &(a)); +#define EXP_FIX(x) EC_FIX(c, &(x), &(x)); #define EXP_SETMUL(d, x, y) do { \ EC_CREATE(&(d)); \ diff --git a/ec-prime.c b/ec-prime.c index b2bfd52..14e4c16 100644 --- a/ec-prime.c +++ b/ec-prime.c @@ -1,6 +1,6 @@ /* -*-c-*- * - * $Id: ec-prime.c,v 1.3.4.1 2003/06/10 13:43:53 mdw Exp $ + * $Id: ec-prime.c,v 1.3.4.2 2004/03/20 00:13:31 mdw Exp $ * * Elliptic curves over prime fields * @@ -30,6 +30,9 @@ /*----- Revision history --------------------------------------------------* * * $Log: ec-prime.c,v $ + * Revision 1.3.4.2 2004/03/20 00:13:31 mdw + * Projective coordinates for prime curves + * * Revision 1.3.4.1 2003/06/10 13:43:53 mdw * Simple (non-projective) curves over prime fields now seem to work. * @@ -59,7 +62,7 @@ typedef struct ecctx { /*----- Simple prime curves -----------------------------------------------*/ -static const ec_ops ec_primeops; +static const ec_ops ec_primeops, ec_primeprojops, ec_primeprojxops; static ec *ecneg(ec_curve *c, ec *d, const ec *p) { @@ -68,11 +71,33 @@ static ec *ecneg(ec_curve *c, ec *d, const ec *p) return (d); } +static ec *ecfind(ec_curve *c, ec *d, mp *x) +{ + mp *p, *q; + ecctx *cc = (ecctx *)c; + field *f = c->f; + + q = F_SQR(f, MP_NEW, x); + p = F_MUL(f, MP_NEW, x, q); + q = F_MUL(f, q, x, cc->a); + p = F_ADD(f, p, p, q); + p = F_ADD(f, p, p, cc->b); + MP_DROP(q); + p = F_SQRT(f, p, p); + if (!p) + return (0); + EC_DESTROY(d); + d->x = MP_COPY(x); + d->y = p; + d->z = MP_COPY(f->one); + return (d); +} + static ec *ecdbl(ec_curve *c, ec *d, const ec *a) { if (EC_ATINF(a)) EC_SETINF(d); - else if (!MP_LEN(a->y)) + else if (F_ZEROP(c->f, a->y)) EC_COPY(d, a); else { field *f = c->f; @@ -80,19 +105,19 @@ static ec *ecdbl(ec_curve *c, ec *d, const ec *a) mp *lambda; mp *dy, *dx; - dx = F_SQR(f, MP_NEW, a->x); - dy = F_DBL(f, MP_NEW, a->y); - dx = F_TPL(f, dx, dx); - dx = F_ADD(f, dx, dx, cc->a); - dy = F_INV(f, dy, dy); - lambda = F_MUL(f, MP_NEW, dx, dy); + dx = F_SQR(f, MP_NEW, a->x); /* %$x^2$% */ + dy = F_DBL(f, MP_NEW, a->y); /* %$2 y$% */ + dx = F_TPL(f, dx, dx); /* %$3 x^2$% */ + dx = F_ADD(f, dx, dx, cc->a); /* %$3 x^2 + A$% */ + dy = F_INV(f, dy, dy); /* %$(2 y)^{-1}$% */ + lambda = F_MUL(f, MP_NEW, dx, dy); /* %$\lambda = (3 x^2 + A)/(2 y)$% */ - dx = F_SQR(f, dx, lambda); - dy = F_DBL(f, dy, a->x); - dx = F_SUB(f, dx, dx, dy); - dy = F_SUB(f, dy, a->x, dx); - dy = F_MUL(f, dy, lambda, dy); - dy = F_SUB(f, dy, dy, a->y); + dx = F_SQR(f, dx, lambda); /* %$\lambda^2$% */ + dy = F_DBL(f, dy, a->x); /* %$2 x$% */ + dx = F_SUB(f, dx, dx, dy); /* %$x' = \lambda^2 - 2 x */ + dy = F_SUB(f, dy, a->x, dx); /* %$x - x'$% */ + dy = F_MUL(f, dy, lambda, dy); /* %$\lambda (x - x')$% */ + dy = F_SUB(f, dy, dy, a->y); /* %$y' = \lambda (x - x') - y$% */ EC_DESTROY(d); d->x = dx; @@ -103,6 +128,94 @@ static ec *ecdbl(ec_curve *c, ec *d, const ec *a) return (d); } +static ec *ecprojdbl(ec_curve *c, ec *d, const ec *a) +{ + if (EC_ATINF(a)) + EC_SETINF(d); + else if (F_ZEROP(c->f, a->y)) + EC_COPY(d, a); + else { + field *f = c->f; + ecctx *cc = (ecctx *)c; + mp *p, *q, *m, *s, *dx, *dy, *dz; + + p = F_SQR(f, MP_NEW, a->z); /* %$z^2$% */ + q = F_SQR(f, MP_NEW, p); /* %$z^4$% */ + p = F_MUL(f, p, q, cc->a); /* %$A z^4$% */ + m = F_SQR(f, MP_NEW, a->x); /* %$x^2$% */ + m = F_TPL(f, m, m); /* %$3 x^2$% */ + m = F_ADD(f, m, m, p); /* %$m = 3 x^2 + A z^4$% */ + + q = F_DBL(f, q, a->y); /* %$2 y$% */ + dz = F_MUL(f, MP_NEW, q, a->z); /* %$z' = 2 y z$% */ + + p = F_SQR(f, p, q); /* %$4 y^2$% */ + s = F_MUL(f, MP_NEW, p, a->x); /* %$s = 4 x y^2$% */ + q = F_SQR(f, q, p); /* %$16 y^4$% */ + q = F_HLV(f, q, q); /* %$t = 8 y^4$% */ + + p = F_DBL(f, p, s); /* %$2 s$% */ + dx = F_SQR(f, MP_NEW, m); /* %$m^2$% */ + dx = F_SUB(f, dx, dx, p); /* %$x' = m^2 - 2 s$% */ + + s = F_SUB(f, s, s, dx); /* %$s - x'$% */ + dy = F_MUL(f, p, m, s); /* %$m (s - x')$% */ + dy = F_SUB(f, dy, dy, q); /* %$y' = m (s - x') - t$% */ + + EC_DESTROY(d); + d->x = dx; + d->y = dy; + d->z = dz; + MP_DROP(m); + MP_DROP(q); + MP_DROP(s); + } + return (d); +} + +static ec *ecprojxdbl(ec_curve *c, ec *d, const ec *a) +{ + if (EC_ATINF(a)) + EC_SETINF(d); + else if (F_ZEROP(c->f, a->y)) + EC_COPY(d, a); + else { + field *f = c->f; + mp *p, *q, *m, *s, *dx, *dy, *dz; + + m = F_SQR(f, MP_NEW, a->z); /* %$z^2$% */ + p = F_SUB(f, MP_NEW, a->x, m); /* %$x - z^2$% */ + q = F_ADD(f, MP_NEW, a->x, m); /* %$x + z^2$% */ + m = F_MUL(f, m, p, q); /* %$x^2 - z^4$% */ + m = F_TPL(f, m, m); /* %$m = 3 x^2 - 3 z^4$% */ + + q = F_DBL(f, q, a->y); /* %$2 y$% */ + dz = F_MUL(f, MP_NEW, q, a->z); /* %$z' = 2 y z$% */ + + p = F_SQR(f, p, q); /* %$4 y^2$% */ + s = F_MUL(f, MP_NEW, p, a->x); /* %$s = 4 x y^2$% */ + q = F_SQR(f, q, p); /* %$16 y^4$% */ + q = F_HLV(f, q, q); /* %$t = 8 y^4$% */ + + p = F_DBL(f, p, s); /* %$2 s$% */ + dx = F_SQR(f, MP_NEW, m); /* %$m^2$% */ + dx = F_SUB(f, dx, dx, p); /* %$x' = m^2 - 2 s$% */ + + s = F_SUB(f, s, s, dx); /* %$s - x'$% */ + dy = F_MUL(f, p, m, s); /* %$m (s - x')$% */ + dy = F_SUB(f, dy, dy, q); /* %$y' = m (s - x') - t$% */ + + EC_DESTROY(d); + d->x = dx; + d->y = dy; + d->z = dz; + MP_DROP(m); + MP_DROP(q); + MP_DROP(s); + } + return (d); +} + static ec *ecadd(ec_curve *c, ec *d, const ec *a, const ec *b) { if (a == b) @@ -117,29 +230,32 @@ static ec *ecadd(ec_curve *c, ec *d, const ec *a, const ec *b) mp *dy, *dx; if (!MP_EQ(a->x, b->x)) { - dy = F_SUB(f, MP_NEW, a->y, b->y); - dx = F_SUB(f, MP_NEW, a->x, b->x); - dx = F_INV(f, dx, dx); + dy = F_SUB(f, MP_NEW, a->y, b->y); /* %$y_0 - y_1$% */ + dx = F_SUB(f, MP_NEW, a->x, b->x); /* %$x_0 - x_1$% */ + dx = F_INV(f, dx, dx); /* %$(x_0 - x_1)^{-1}$% */ lambda = F_MUL(f, MP_NEW, dy, dx); - } else if (!MP_LEN(a->y) || !MP_EQ(a->y, b->y)) { + /* %$\lambda = (y_0 - y1)/(x_0 - x_1)$% */ + } else if (F_ZEROP(c->f, a->y) || !MP_EQ(a->y, b->y)) { EC_SETINF(d); return (d); } else { ecctx *cc = (ecctx *)c; - dx = F_SQR(f, MP_NEW, a->x); - dx = F_TPL(f, dx, dx); - dx = F_ADD(f, dx, dx, cc->a); - dy = F_DBL(f, MP_NEW, a->y); - dy = F_INV(f, dy, dy); + dx = F_SQR(f, MP_NEW, a->x); /* %$x_0^2$% */ + dx = F_TPL(f, dx, dx); /* %$3 x_0^2$% */ + dx = F_ADD(f, dx, dx, cc->a); /* %$3 x_0^2 + A$% */ + dy = F_DBL(f, MP_NEW, a->y); /* %$2 y_0$% */ + dy = F_INV(f, dy, dy); /* %$(2 y_0)^{-1}$% */ lambda = F_MUL(f, MP_NEW, dx, dy); + /* %$\lambda = (3 x_0^2 + A)/(2 y_0)$% */ } - dx = F_SQR(f, dx, lambda); - dx = F_SUB(f, dx, dx, a->x); - dx = F_SUB(f, dx, dx, b->x); - dy = F_SUB(f, dy, b->x, dx); - dy = F_MUL(f, dy, lambda, dy); + dx = F_SQR(f, dx, lambda); /* %$\lambda^2$% */ + dx = F_SUB(f, dx, dx, a->x); /* %$\lambda^2 - x_0$% */ + dx = F_SUB(f, dx, dx, b->x); /* %$x' = \lambda^2 - x_0 - x_1$% */ + dy = F_SUB(f, dy, b->x, dx); /* %$x_1 - x'$% */ + dy = F_MUL(f, dy, lambda, dy); /* %$\lambda (x_1 - x')$% */ dy = F_SUB(f, dy, dy, b->y); + /* %$y' = \lambda (x_1 - x') - y_1$% */ EC_DESTROY(d); d->x = dx; @@ -150,6 +266,101 @@ static ec *ecadd(ec_curve *c, ec *d, const ec *a, const ec *b) return (d); } +static ec *ecprojadd(ec_curve *c, ec *d, const ec *a, const ec *b) +{ + if (a == b) + c->ops->dbl(c, d, a); + else if (EC_ATINF(a)) + EC_COPY(d, b); + else if (EC_ATINF(b)) + EC_COPY(d, a); + else { + field *f = c->f; + mp *p, *q, *r, *w, *u, *s, *dx, *dy, *dz; + + q = F_SQR(f, MP_NEW, a->z); /* %$z_0^2$% */ + u = F_MUL(f, MP_NEW, q, b->x); /* %$u = x_1 z_0^2$% */ + p = F_MUL(f, MP_NEW, q, b->y); /* %$y_1 z_0^2$% */ + s = F_MUL(f, q, p, a->z); /* %$s = y_1 z_0^3$% */ + + w = F_SUB(f, p, a->x, u); /* %$w = x_0 - u$% */ + r = F_SUB(f, MP_NEW, a->y, s); /* %$r = y_0 - s$% */ + if (F_ZEROP(f, w)) { + if (F_ZEROP(f, r)) { + MP_DROP(w); + MP_DROP(r); + MP_DROP(u); + MP_DROP(s); + return (c->ops->dbl(c, d, a)); + } else { + MP_DROP(w); + MP_DROP(r); + MP_DROP(u); + MP_DROP(s); + EC_SETINF(d); + return (d); + } + } + u = F_ADD(f, u, u, a->x); /* %$t = x_0 + u$% */ + s = F_ADD(f, s, s, a->y); /* %$m = y_0 + r$% */ + + dz = F_MUL(f, MP_NEW, a->z, w); /* %$z' = z_0 w$% */ + + p = F_SQR(f, MP_NEW, w); /* %$w^2$% */ + q = F_MUL(f, MP_NEW, p, u); /* %$t w^2$% */ + u = F_MUL(f, u, p, w); /* %$w^3$% */ + p = F_MUL(f, p, u, s); /* %$m w^3$% */ + + dx = F_SQR(f, u, r); /* %$r^2$% */ + dx = F_SUB(f, dx, dx, q); /* %$x' = r^2 - t w^2$% */ + + s = F_DBL(f, s, dx); /* %$2 x'$% */ + q = F_SUB(f, q, q, s); /* %$v = t w^2 - 2 x'$% */ + dy = F_MUL(f, s, q, r); /* %$v r$% */ + dy = F_SUB(f, dy, dy, p); /* %$v r - m w^3$% */ + dy = F_HLV(f, dy, dy); /* %$y' = (v r - m w^3)/2$% */ + + EC_DESTROY(d); + d->x = dx; + d->y = dy; + d->z = dz; + MP_DROP(p); + MP_DROP(q); + MP_DROP(r); + MP_DROP(w); + } + return (d); +} + +static int eccheck(ec_curve *c, const ec *p) +{ + ecctx *cc = (ecctx *)c; + field *f = c->f; + int rc; + mp *l = F_SQR(f, MP_NEW, p->y); + mp *x = F_SQR(f, MP_NEW, p->x); + mp *r = F_MUL(f, MP_NEW, x, p->x); + x = F_MUL(f, x, cc->a, p->x); + r = F_ADD(f, r, r, x); + r = F_ADD(f, r, r, cc->b); + rc = MP_EQ(l, r) ? 0 : -1; + mp_drop(l); + mp_drop(x); + mp_drop(r); + return (rc); +} + +static int ecprojcheck(ec_curve *c, const ec *p) +{ + ec t = EC_INIT; + int rc; + + c->ops->fix(c, &t, p); + rc = eccheck(c, &t); + EC_DESTROY(&t); + return (rc); +} + static void ecdestroy(ec_curve *c) { ecctx *cc = (ecctx *)c; @@ -180,8 +391,37 @@ extern ec_curve *ec_prime(field *f, mp *a, mp *b) return (&cc->c); } +extern ec_curve *ec_primeproj(field *f, mp *a, mp *b) +{ + ecctx *cc = CREATE(ecctx); + mp *ax; + + ax = mp_add(MP_NEW, a, MP_THREE); + ax = F_IN(f, ax, ax); + if (F_ZEROP(f, ax)) + cc->c.ops = &ec_primeprojxops; + else + cc->c.ops = &ec_primeprojops; + MP_DROP(ax); + cc->c.f = f; + cc->a = F_IN(f, MP_NEW, a); + cc->b = F_IN(f, MP_NEW, b); + return (&cc->c); +} + static const ec_ops ec_primeops = { - ecdestroy, ec_idin, ec_idout, 0, ecneg, ecadd, ec_stdsub, ecdbl + ecdestroy, ec_idin, ec_idout, ec_idfix, + 0, ecneg, ecadd, ec_stdsub, ecdbl, eccheck +}; + +static const ec_ops ec_primeprojops = { + ecdestroy, ec_projin, ec_projout, ec_projfix, + 0, ecneg, ecprojadd, ec_stdsub, ecprojdbl, ecprojcheck +}; + +static const ec_ops ec_primeprojxops = { + ecdestroy, ec_projin, ec_projout, ec_projfix, + 0, ecneg, ecprojadd, ec_stdsub, ecprojxdbl, ecprojcheck }; /*----- Test rig ----------------------------------------------------------*/ diff --git a/ec.c b/ec.c index 4111928..a2b229f 100644 --- a/ec.c +++ b/ec.c @@ -1,6 +1,6 @@ /* -*-c-*- * - * $Id: ec.c,v 1.4.4.1 2003/06/10 13:43:53 mdw Exp $ + * $Id: ec.c,v 1.4.4.2 2004/03/20 00:13:31 mdw Exp $ * * Elliptic curve definitions * @@ -30,6 +30,9 @@ /*----- Revision history --------------------------------------------------* * * $Log: ec.c,v $ + * Revision 1.4.4.2 2004/03/20 00:13:31 mdw + * Projective coordinates for prime curves + * * Revision 1.4.4.1 2003/06/10 13:43:53 mdw * Simple (non-projective) curves over prime fields now seem to work. * @@ -113,7 +116,7 @@ ec *ec_copy(ec *d, const ec *p) { EC_COPY(d, p); return (d); } /*----- Standard curve operations -----------------------------------------*/ -/* --- @ec_idin@, @ec_idout@ --- * +/* --- @ec_idin@, @ec_idout@, @ec_idfix@ --- * * * Arguments: @ec_curve *c@ = pointer to an elliptic curve * @ec *d@ = pointer to the destination @@ -152,6 +155,12 @@ ec *ec_idout(ec_curve *c, ec *d, const ec *p) return (d); } +ec *ec_idfix(ec_curve *c, ec *d, const ec *p) +{ + EC_COPY(d, p); + return (d); +} + /* --- @ec_projin@, @ec_projout@ --- * * * Arguments: @ec_curve *c@ = pointer to an elliptic curve @@ -182,12 +191,15 @@ ec *ec_projout(ec_curve *c, ec *d, const ec *p) if (EC_ATINF(p)) EC_SETINF(d); else { - mp *x, *y, *z; + mp *x, *y, *z, *zz; field *f = c->f; z = F_INV(f, MP_NEW, p->z); - x = F_MUL(f, d->x, p->x, z); + zz = F_SQR(f, MP_NEW, z); + z = F_MUL(f, z, zz, z); + x = F_MUL(f, d->x, p->x, zz); y = F_MUL(f, d->y, p->y, z); mp_drop(z); + mp_drop(zz); mp_drop(d->z); d->x = F_OUT(f, x, x); d->y = F_OUT(f, y, y); @@ -196,6 +208,28 @@ ec *ec_projout(ec_curve *c, ec *d, const ec *p) return (d); } +ec *ec_projfix(ec_curve *c, ec *d, const ec *p) +{ + if (EC_ATINF(p)) + EC_SETINF(d); + else if (d->z == c->f->one) + EC_COPY(d, p); + else { + mp *z, *zz; + field *f = c->f; + z = F_INV(f, MP_NEW, p->z); + zz = F_SQR(f, MP_NEW, z); + z = F_MUL(f, z, zz, z); + d->x = F_MUL(f, d->x, p->x, zz); + d->y = F_MUL(f, d->y, p->y, z); + mp_drop(z); + mp_drop(zz); + mp_drop(d->z); + d->z = MP_COPY(f->one); + } + return (d); +} + /* --- @ec_stdsub@ --- * * * Arguments: @ec_curve *c@ = pointer to an elliptic curve @@ -213,6 +247,7 @@ ec *ec_stdsub(ec_curve *c, ec *d, const ec *p, const ec *q) { ec t = EC_INIT; EC_NEG(c, &t, q); + EC_FIX(c, &t, &t); EC_ADD(c, d, p, &t); EC_DESTROY(&t); return (d); @@ -249,7 +284,7 @@ ec *ec_find(ec_curve *c, ec *d, mp *x) x = F_IN(c->f, MP_NEW, x); if ((d = EC_FIND(c, d, x)) != 0) EC_OUT(c, d, d); - mp_drop(x); + MP_DROP(x); return (d); } @@ -335,6 +370,29 @@ ec *ec_dbl(ec_curve *c, ec *d, const ec *p) return (EC_OUT(c, d, d)); } +/* --- @ec_check@ --- * + * + * Arguments: @ec_curve *c@ = pointer to an elliptic curve + * @const ec *p@ = pointer to the point + * + * Returns: Zero if OK, nonzero if this is an invalid point. + * + * Use: Checks that a point is actually on an elliptic curve. + */ + +int ec_check(ec_curve *c, const ec *p) +{ + ec t = EC_INIT; + int rc; + + if (EC_ATINF(p)) + return (0); + EC_IN(c, &t, p); + rc = EC_CHECK(c, &t); + EC_DESTROY(&t); + return (rc); +} + /* --- @ec_imul@, @ec_mul@ --- * * * Arguments: @ec_curve *c@ = pointer to an elliptic curve @@ -360,10 +418,14 @@ ec *ec_imul(ec_curve *c, ec *d, const ec *p, mp *n) EC_SETINF(d); if (MP_LEN(n) == 0) ; - else if (MP_LEN(n) < EXP_THRESH) - EXP_SIMPLE(*d, t, n); - else - EXP_WINDOW(*d, t, n); + else { + if (n->f & MP_NEG) + EC_NEG(c, &t, &t); + if (MP_LEN(n) < EXP_THRESH) + EXP_SIMPLE(*d, t, n); + else + EXP_WINDOW(*d, t, n); + } EC_DESTROY(&t); return (d); } diff --git a/ec.h b/ec.h index 72ee7a9..476d898 100644 --- a/ec.h +++ b/ec.h @@ -1,6 +1,6 @@ /* -*-c-*- * - * $Id: ec.h,v 1.4.4.1 2003/06/10 13:43:53 mdw Exp $ + * $Id: ec.h,v 1.4.4.2 2004/03/20 00:13:31 mdw Exp $ * * Elliptic curve definitions * @@ -30,6 +30,9 @@ /*----- Revision history --------------------------------------------------* * * $Log: ec.h,v $ + * Revision 1.4.4.2 2004/03/20 00:13:31 mdw + * Projective coordinates for prime curves + * * Revision 1.4.4.1 2003/06/10 13:43:53 mdw * Simple (non-projective) curves over prime fields now seem to work. * @@ -83,27 +86,37 @@ typedef struct ec_mulfactor { mp *exp; /* The exponent */ } ec_mulfactor; -/* --- Elliptic curve operations --- */ +/* --- Elliptic curve operations --- * + * + * All operations (apart from @destroy@ and @in@) are guaranteed to be + * performed on internal representations of points. Moreover, the second + * argument to @add@ and @mul@ is guaranteed to be the output of @in@ or + * @fix@. + */ typedef struct ec_ops { void (*destroy)(ec_curve */*c*/); ec *(*in)(ec_curve */*c*/, ec */*d*/, const ec */*p*/); ec *(*out)(ec_curve */*c*/, ec */*d*/, const ec */*p*/); + ec *(*fix)(ec_curve */*c*/, ec */*d*/, const ec */*p*/); ec *(*find)(ec_curve */*c*/, ec */*d*/, mp */*x*/); ec *(*neg)(ec_curve */*c*/, ec */*d*/, const ec */*p*/); ec *(*add)(ec_curve */*c*/, ec */*d*/, const ec */*p*/, const ec */*q*/); ec *(*sub)(ec_curve */*c*/, ec */*d*/, const ec */*p*/, const ec */*q*/); ec *(*dbl)(ec_curve */*c*/, ec */*d*/, const ec */*p*/); + int (*check)(ec_curve */*c*/, const ec */*p*/); } ec_ops; #define EC_IN(c, d, p) (c)->ops->in((c), (d), (p)) #define EC_OUT(c, d, p) (c)->ops->out((c), (d), (p)) +#define EC_FIX(c, d, p) (c)->ops->fix((c), (d), (p)) #define EC_FIND(c, d, x) (c)->ops->find((c), (d), (x)) #define EC_NEG(c, d, x) (c)->ops->neg((c), (d), (x)) #define EC_ADD(c, d, p, q) (c)->ops->add((c), (d), (p), (q)) #define EC_SUB(c, d, p, q) (c)->ops->sub((c), (d), (p), (q)) #define EC_DBL(c, d, p) (c)->ops->dbl((c), (d), (p)) +#define EC_CHECK(c, p) (c)->ops->check((c), (p)) /*----- Simple memory management things -----------------------------------*/ @@ -210,32 +223,6 @@ extern ec *ec_copy(ec */*d*/, const ec */*p*/); /*----- Interesting arithmetic --------------------------------------------*/ -/* --- @ec_in@ --- * - * - * Arguments: @ec_curve *c@ = pointer to an elliptic curve - * @ec *d@ = pointer to the destination point - * @const ec *p@ = pointer to the source point - * - * Returns: The destination point. - * - * Use: Converts a point to internal representation. - */ - -extern ec *ec_in(ec_curve */*c*/, ec */*d*/, const ec */*p*/); - -/* --- @ec_out@ --- * - * - * Arguments: @ec_curve *c@ = pointer to an elliptic curve - * @ec *d@ = pointer to the destination point - * @const ec *p@ = pointer to the source point - * - * Returns: The destination point. - * - * Use: Converts a point to external representation. - */ - -extern ec *ec_out(ec_curve */*c*/, ec */*d*/, const ec */*p*/); - /* --- @ec_find@ --- * * * Arguments: @ec_curve *c@ = pointer to an elliptic curve @@ -306,6 +293,18 @@ extern ec *ec_sub(ec_curve */*c*/, ec */*d*/, extern ec *ec_dbl(ec_curve */*c*/, ec */*d*/, const ec */*p*/); +/* --- @ec_check@ --- * + * + * Arguments: @ec_curve *c@ = pointer to an elliptic curve + * @const ec *p@ = pointer to the point + * + * Returns: Zero if OK, nonzero if this is an invalid point. + * + * Use: Checks that a point is actually on an elliptic curve. + */ + +extern int ec_check(ec_curve */*c*/, const ec */*p*/); + /* --- @ec_mul@, @ec_imul@ --- * * * Arguments: @ec_curve *c@ = pointer to an elliptic curve @@ -343,7 +342,7 @@ extern ec *ec_immul(ec_curve */*c*/, ec */*d*/, /*----- Standard curve operations -----------------------------------------*/ -/* --- @ec_idin@, @ec_idout@ --- * +/* --- @ec_idin@, @ec_idout@, @ec_idfix@ --- * * * Arguments: @ec_curve *c@ = pointer to an elliptic curve * @ec *d@ = pointer to the destination @@ -358,8 +357,9 @@ extern ec *ec_immul(ec_curve */*c*/, ec */*d*/, extern ec *ec_idin(ec_curve */*c*/, ec */*d*/, const ec */*p*/); extern ec *ec_idout(ec_curve */*c*/, ec */*d*/, const ec */*p*/); +extern ec *ec_idfix(ec_curve */*c*/, ec */*d*/, const ec */*p*/); -/* --- @ec_projin@, @ec_projout@ --- * +/* --- @ec_projin@, @ec_projout@, @ec_projfix@ --- * * * Arguments: @ec_curve *c@ = pointer to an elliptic curve * @ec *d@ = pointer to the destination @@ -373,6 +373,7 @@ extern ec *ec_idout(ec_curve */*c*/, ec */*d*/, const ec */*p*/); extern ec *ec_projin(ec_curve */*c*/, ec */*d*/, const ec */*p*/); extern ec *ec_projout(ec_curve */*c*/, ec */*d*/, const ec */*p*/); +extern ec *ec_projfix(ec_curve */*c*/, ec */*d*/, const ec */*p*/); /* --- @ec_stdsub@ --- * * diff --git a/exp.h b/exp.h index 6cfdfd8..fc9e3a9 100644 --- a/exp.h +++ b/exp.h @@ -1,6 +1,6 @@ /* -*-c-*- * - * $Id: exp.h,v 1.1 2001/06/16 13:00:59 mdw Exp $ + * $Id: exp.h,v 1.1.4.1 2004/03/20 00:13:31 mdw Exp $ * * Generalized exponentiation * @@ -30,6 +30,9 @@ /*----- Revision history --------------------------------------------------* * * $Log: exp.h,v $ + * Revision 1.1.4.1 2004/03/20 00:13:31 mdw + * Projective coordinates for prime curves + * * Revision 1.1 2001/06/16 13:00:59 mdw * New generic exponentation code. Includes sliding-window simultaneous * exponentiation. @@ -99,6 +102,10 @@ typedef struct exp_simul { * @EXP_MUL(a, x)@ Multiplies @a@ by @x@ (writing the result * back to @a@). * + * @EXP_FIX(x)@ Makes @x@ be a canonical representation of + * its value. All multiplications have the + * right argument canonical. + * * @EXP_SQR(a)@ Multiplies @a@ by itself. * * @EXP_SETMUL(d, x, y)@ Sets @d@ to be the product of @x@ and @y@. @@ -140,6 +147,7 @@ typedef struct exp_simul { \ /* --- Do the main body of the work --- */ \ \ + EXP_FIX(g); \ for (;;) { \ EXP_MUL(a, g); \ sq = 0; \ @@ -184,11 +192,15 @@ exp_simple_exit:; \ \ /* --- Do the precomputation --- */ \ \ + EXP_FIX(g); \ EXP_SETSQR(g2, g); \ + EXP_FIX(g2); \ v = xmalloc(EXP_TABSZ * sizeof(EXP_TYPE)); \ EXP_COPY(v[0], g); \ - for (i = 1; i < EXP_TABSZ; i++) \ + for (i = 1; i < EXP_TABSZ; i++) { \ EXP_SETMUL(v[i], v[i - 1], g2); \ + EXP_FIX(v[i]); \ + } \ EXP_DROP(g2); \ \ /* --- Skip top-end zero bits --- * \ @@ -286,17 +298,21 @@ exp_window_exit:; \ j = 1; \ for (i = 0; i < n; i++) { \ EXP_COPY(v[j], f[n - 1 - i].base); \ + EXP_FIX(v[j]); \ j <<= 1; \ } \ k = n * EXP_WINSZ; \ jj = 1; \ for (; i < k; i++) { \ EXP_SETSQR(v[j], v[jj]); \ + EXP_FIX(v[j]); \ j <<= 1; jj <<= 1; \ } \ for (i = 1; i < vn; i <<= 1) { \ - for (j = 1; j < i; j++) \ + for (j = 1; j < i; j++) { \ EXP_SETMUL(v[j + i], v[j], v[i]); \ + EXP_FIX(v[j + i]); \ + } \ } \ \ /* --- Set up the bitscanners --- * \ @@ -381,7 +397,7 @@ exp_window_exit:; \ \ exp_simul_done: \ while (sq--) EXP_SQR(a); \ - for (i = 1; i < vn; i++) \ + for (i = 1; i < vn; i++) \ EXP_DROP(v[i]); \ xfree(v); \ } while (0) diff --git a/f-prime.c b/f-prime.c index 0d76da0..7c1dae5 100644 --- a/f-prime.c +++ b/f-prime.c @@ -1,6 +1,6 @@ /* -*-c-*- * - * $Id: f-prime.c,v 1.3.4.1 2003/06/10 13:43:53 mdw Exp $ + * $Id: f-prime.c,v 1.3.4.2 2004/03/20 00:13:31 mdw Exp $ * * Prime fields with Montgomery arithmetic * @@ -30,6 +30,9 @@ /*----- Revision history --------------------------------------------------* * * $Log: f-prime.c,v $ + * Revision 1.3.4.2 2004/03/20 00:13:31 mdw + * Projective coordinates for prime curves + * * Revision 1.3.4.1 2003/06/10 13:43:53 mdw * Simple (non-projective) curves over prime fields now seem to work. * @@ -82,6 +85,11 @@ static mp *fout(field *ff, mp *d, mp *x) return (mpmont_reduce(&f->mm, d, x)); } +static int fzerop(field *ff, mp *x) +{ + return (!MP_LEN(x)); +} + static mp *fneg(field *ff, mp *d, mp *x) { fctx *f = (fctx *)ff; @@ -157,11 +165,37 @@ static mp *ftpl(field *ff, mp *d, mp *x) return (d); } +static mp *fqdl(field *ff, mp *d, mp *x) +{ + fctx *f = (fctx *)ff; + d = mp_lsl(d, x, 2); + while (MP_CMP(d, >, f->mm.m)) + d = mp_sub(d, d, f->mm.m); + return (d); +} + +static mp *fhlv(field *ff, mp *d, mp *x) +{ + fctx *f = (fctx *)ff; + if (!MP_LEN(x)) { + MP_COPY(x); + MP_DROP(d); + return (x); + } + if (x->v[0] & 1) { + d = mp_add(d, x, f->mm.m); + x = d; + } + return (mp_lsr(d, x, 1)); +} + static mp *fsqrt(field *ff, mp *d, mp *x) { fctx *f = (fctx *)ff; d = mpmont_reduce(&f->mm, d, x); d = mp_modsqrt(d, d, f->mm.m); + if (!d) + return (d); return (mpmont_mul(&f->mm, d, d, f->mm.r2)); } @@ -170,8 +204,8 @@ static mp *fsqrt(field *ff, mp *d, mp *x) static field_ops fops = { fdestroy, fin, fout, - fneg, fadd, fsub, fmul, fsqr, finv, freduce, - fdbl, ftpl, fsqrt + fzerop, fneg, fadd, fsub, fmul, fsqr, finv, freduce, + fdbl, ftpl, fqdl, fhlv, fsqrt }; /* --- @field_prime@ --- * diff --git a/field.h b/field.h index 5a70649..8096700 100644 --- a/field.h +++ b/field.h @@ -1,6 +1,6 @@ /* -*-c-*- * - * $Id: field.h,v 1.3 2002/01/13 13:48:44 mdw Exp $ + * $Id: field.h,v 1.3.4.1 2004/03/20 00:13:31 mdw Exp $ * * Definitions for field arithmetic * @@ -30,6 +30,9 @@ /*----- Revision history --------------------------------------------------* * * $Log: field.h,v $ + * Revision 1.3.4.1 2004/03/20 00:13:31 mdw + * Projective coordinates for prime curves + * * Revision 1.3 2002/01/13 13:48:44 mdw * Further progress. * @@ -70,6 +73,7 @@ typedef struct field_ops { mp *(*in)(field */*f*/, mp */*d*/, mp */*x*/); mp *(*out)(field */*f*/, mp */*d*/, mp */*x*/); + int (*zerop)(field */*f*/, mp */*x*/); mp *(*neg)(field */*f*/, mp */*d*/, mp */*x*/); mp *(*add)(field */*f*/, mp */*d*/, mp */*x*/, mp */*y*/); mp *(*sub)(field */*f*/, mp */*d*/, mp */*x*/, mp */*y*/); @@ -82,6 +86,8 @@ typedef struct field_ops { mp *(*dbl)(field */*f*/, mp */*d*/, mp */*x*/); mp *(*tpl)(field */*f*/, mp */*d*/, mp */*x*/); + mp *(*qdl)(field */*f*/, mp */*d*/, mp */*x*/); + mp *(*hlv)(field */*f*/, mp */*d*/, mp */*x*/); mp *(*sqrt)(field */*f*/, mp */*d*/, mp */*x*/); } field_ops; @@ -90,6 +96,8 @@ typedef struct field_ops { #define F_IN(f, d, x) (f)->ops->in((f), (d), (x)) #define F_OUT(f, d, x) (f)->ops->out((f), (d), (x)) + +#define F_ZEROP(f, x) (f)->ops->zerop((f), (x)) #define F_NEG(f, d, x) (f)->ops->neg((f), (d), (x)) #define F_ADD(f, d, x, y) (f)->ops->add((f), (d), (x), (y)) #define F_SUB(f, d, x, y) (f)->ops->sub((f), (d), (x), (y)) @@ -100,6 +108,8 @@ typedef struct field_ops { #define F_DBL(f, d, x) (f)->ops->dbl((f), (d), (x)) #define F_TPL(f, d, x) (f)->ops->tpl((f), (d), (x)) +#define F_QDL(f, d, x) (f)->ops->qdl((f), (d), (x)) +#define F_HLV(f, d, x) (f)->ops->hlv((f), (d), (x)) #define F_SQRT(f, d, x) (f)->ops->sqrt((f), (d), (x)) /*----- Helpful field operations ------------------------------------------*/ diff --git a/mpmont-exp.h b/mpmont-exp.h index e732e6f..0f82b9c 100644 --- a/mpmont-exp.h +++ b/mpmont-exp.h @@ -1,6 +1,6 @@ /* -*-c-*- * - * $Id: mpmont-exp.h,v 1.1 2001/06/16 12:58:12 mdw Exp $ + * $Id: mpmont-exp.h,v 1.1.4.1 2004/03/20 00:13:31 mdw Exp $ * * Exponentiation operations for Montgomery reduction * @@ -30,6 +30,9 @@ /*----- Revision history --------------------------------------------------* * * $Log: mpmont-exp.h,v $ + * Revision 1.1.4.1 2004/03/20 00:13:31 mdw + * Projective coordinates for prime curves + * * Revision 1.1 2001/06/16 12:58:12 mdw * Parameters for generic exponentiation. * @@ -61,6 +64,8 @@ a = mpmont_reduce(mm, t, t); \ } while (0) +#define EXP_FIX(x) + #define EXP_SETMUL(d, x, y) d = mpmont_mul(mm, MP_NEW, x, y) #define EXP_SETSQR(d, x) do { \ -- 2.11.0