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

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