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

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