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

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