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

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