chiark - git - mdw - catacomb/blob - ec-prime.c

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