git.distorted.org.uk Git - u/mdw/catacomb/blob - ec-prime.c

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