mdw@git.distorted.org.uk Git - catacomb/blob - math/ec-prime.c

   1 /* -*-c-*-
   2  *
   3  * Elliptic curves over prime fields
   4  *
   5  * (c) 2001 Straylight/Edgeware
   6  */
   7
   8 /*----- Licensing notice --------------------------------------------------*
   9  *
  10  * This file is part of Catacomb.
  11  *
  12  * Catacomb is free software; you can redistribute it and/or modify
  13  * it under the terms of the GNU Library General Public License as
  14  * published by the Free Software Foundation; either version 2 of the
  15  * License, or (at your option) any later version.
  16  *
  17  * Catacomb is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  20  * GNU Library General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Library General Public
  23  * License along with Catacomb; if not, write to the Free
  24  * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
  25  * MA 02111-1307, USA.
  26  */
  27
  28 /*----- Header files ------------------------------------------------------*/
  29
  30 #include <mLib/sub.h>
  31
  32 #include "ec.h"
  33
  34 /*----- Simple prime curves -----------------------------------------------*/
  35
  36 static const ec_ops ec_primeops, ec_primeprojops, ec_primeprojxops;
  37
  38 static ec *ecneg(ec_curve *c, ec *d, const ec *p)
  39 {
  40   EC_COPY(d, p);
  41   if (d->y)
  42     d->y = F_NEG(c->f, d->y, d->y);
  43   return (d);
  44 }
  45
  46 static ec *ecfind(ec_curve *c, ec *d, mp *x)
  47 {
  48   mp *p, *q;
  49   field *f = c->f;
  50
  51   q = F_SQR(f, MP_NEW, x);
  52   p = F_MUL(f, MP_NEW, x, q);
  53   q = F_MUL(f, q, x, c->a);
  54   p = F_ADD(f, p, p, q);
  55   p = F_ADD(f, p, p, c->b);
  56   MP_DROP(q);
  57   p = F_SQRT(f, p, p);
  58   if (!p)
  59     return (0);
  60   EC_DESTROY(d);
  61   d->x = MP_COPY(x);
  62   d->y = p;
  63   d->z = MP_COPY(f->one);
  64   return (d);
  65 }
  66
  67 static ec *ecdbl(ec_curve *c, ec *d, const ec *a)
  68 {
  69   if (EC_ATINF(a) || F_ZEROP(c->f, a->y))
  70     EC_SETINF(d);
  71   else {
  72     field *f = c->f;
  73     mp *lambda;
  74     mp *dy, *dx;
  75
  76     dx = F_SQR(f, MP_NEW, a->x);        /* %$x^2$% */
  77     dy = F_DBL(f, MP_NEW, a->y);        /* %$2 y$% */
  78     dx = F_TPL(f, dx, dx);              /* %$3 x^2$% */
  79     dx = F_ADD(f, dx, dx, c->a);        /* %$3 x^2 + A$% */
  80     dy = F_INV(f, dy, dy);              /* %$(2 y)^{-1}$% */
  81     lambda = F_MUL(f, MP_NEW, dx, dy);  /* %$\lambda = (3 x^2 + A)/(2 y)$% */
  82
  83     dx = F_SQR(f, dx, lambda);          /* %$\lambda^2$% */
  84     dy = F_DBL(f, dy, a->x);            /* %$2 x$% */
  85     dx = F_SUB(f, dx, dx, dy);          /* %$x' = \lambda^2 - 2 x */
  86     dy = F_SUB(f, dy, a->x, dx);        /* %$x - x'$% */
  87     dy = F_MUL(f, dy, lambda, dy);      /* %$\lambda (x - x')$% */
  88     dy = F_SUB(f, dy, dy, a->y);        /* %$y' = \lambda (x - x') - y$% */
  89
  90     EC_DESTROY(d);
  91     d->x = dx;
  92     d->y = dy;
  93     d->z = 0;
  94     MP_DROP(lambda);
  95   }
  96   return (d);
  97 }
  98
  99 static ec *ecprojdbl(ec_curve *c, ec *d, const ec *a)
 100 {
 101   if (EC_ATINF(a) || F_ZEROP(c->f, a->y))
 102     EC_SETINF(d);
 103   else {
 104     field *f = c->f;
 105     mp *p, *q, *m, *s, *dx, *dy, *dz;
 106
 107     q = F_ADD(f, MP_NEW, a->y, a->z);   /* %$y + z$% */
 108     dz = F_SQR(f, MP_NEW, q);         /* %$(y + z)^2 = y^2 + 2 y z + z^2$% */
 109     p = F_SQR(f, MP_NEW, a->z);         /* %$z^2$% */
 110     dz = F_SUB(f, dz, dz, p);           /* %$y^2 + 2 y z$% */
 111     q = F_SQR(f, q, p);                 /* %$z^4$% */
 112     p = F_MUL(f, p, q, c->a);           /* %$A z^4$% */
 113     m = F_SQR(f, MP_NEW, a->x);         /* %$x^2$% */
 114     m = F_TPL(f, m, m);                 /* %$3 x^2$% */
 115     m = F_ADD(f, m, m, p);              /* %$m = 3 x^2 + A z^4$% */
 116
 117     p = F_SQR(f, p, a->y);              /* %$y^2$% */
 118     dz = F_SUB(f, dz, dz, p);           /* %$z' = 2 y z$% */
 119
 120     p = F_QDL(f, p, p);                 /* %$4 y^2$% */
 121     s = F_MUL(f, MP_NEW, p, a->x);      /* %$s = 4 x y^2$% */
 122     q = F_SQR(f, q, p);                 /* %$16 y^4$% */
 123     q = F_HLV(f, q, q);                 /* %$t = 8 y^4$% */
 124
 125     p = F_DBL(f, p, s);                 /* %$2 s$% */
 126     dx = F_SQR(f, MP_NEW, m);           /* %$m^2$% */
 127     dx = F_SUB(f, dx, dx, p);           /* %$x' = m^2 - 2 s$% */
 128
 129     s = F_SUB(f, s, s, dx);             /* %$s - x'$% */
 130     dy = F_MUL(f, p, m, s);             /* %$m (s - x')$% */
 131     dy = F_SUB(f, dy, dy, q);           /* %$y' = m (s - x') - t$% */
 132
 133     EC_DESTROY(d);
 134     d->x = dx;
 135     d->y = dy;
 136     d->z = dz;
 137     MP_DROP(m);
 138     MP_DROP(q);
 139     MP_DROP(s);
 140   }
 141   return (d);
 142 }
 143
 144 /* #define USE_SQUARE */
 145
 146 static ec *ecprojxdbl(ec_curve *c, ec *d, const ec *a)
 147 {
 148   if (EC_ATINF(a) || F_ZEROP(c->f, a->y))
 149     EC_SETINF(d);
 150   else {
 151     field *f = c->f;
 152     mp *p, *q, *m, *s, *dx, *dy, *dz;
 153
 154 #ifdef USE_SQUARE
 155     q = F_ADD(f, MP_NEW, a->y, a->z);   /* %$y + z$% */
 156     dz = F_SQR(f, MP_NEW, q);         /* %$(y + z)^2 = y^2 + 2 y z + z^2$% */
 157 #else
 158     q = MP_NEW;
 159 #endif
 160     m = F_SQR(f, MP_NEW, a->z);         /* %$z^2$% */
 161 #ifdef USE_SQUARE
 162     dz = F_SUB(f, dz, dz, m);           /* %$y^2 + 2 y z$% */
 163 #endif
 164     p = F_SUB(f, MP_NEW, a->x, m);      /* %$x - z^2$% */
 165     q = F_ADD(f, q, a->x, m);           /* %$x + z^2$% */
 166     m = F_MUL(f, m, p, q);              /* %$x^2 - z^4$% */
 167     m = F_TPL(f, m, m);                 /* %$m = 3 x^2 - 3 z^4$% */
 168
 169 #ifdef USE_SQUARE
 170     p = F_SQR(f, p, a->y);              /* %$y^2$% */
 171     dz = F_SUB(f, dz, dz, p);           /* %$z' = 2 y z$% */
 172     p = F_QDL(f, p, p);                 /* %$4 y^2$% */
 173 #else
 174     q = F_DBL(f, q, a->y);              /* %$2 y% */
 175     dz = F_MUL(f, MP_NEW, a->z, q);     /* %$z' = 2 y z$% */
 176     p = F_SQR(f, p, q);                 /* %$4 y^2% */
 177 #endif
 178
 179     s = F_MUL(f, MP_NEW, p, a->x);      /* %$s = 4 x y^2$% */
 180     q = F_SQR(f, q, p);                 /* %$16 y^4$% */
 181     q = F_HLV(f, q, q);                 /* %$t = 8 y^4$% */
 182
 183     p = F_DBL(f, p, s);                 /* %$2 s$% */
 184     dx = F_SQR(f, MP_NEW, m);           /* %$m^2$% */
 185     dx = F_SUB(f, dx, dx, p);           /* %$x' = m^2 - 2 s$% */
 186
 187     s = F_SUB(f, s, s, dx);             /* %$s - x'$% */
 188     dy = F_MUL(f, p, m, s);             /* %$m (s - x')$% */
 189     dy = F_SUB(f, dy, dy, q);           /* %$y' = m (s - x') - t$% */
 190
 191     EC_DESTROY(d);
 192     d->x = dx;
 193     d->y = dy;
 194     d->z = dz;
 195     MP_DROP(m);
 196     MP_DROP(q);
 197     MP_DROP(s);
 198   }
 199   return (d);
 200 }
 201
 202 static ec *ecadd(ec_curve *c, ec *d, const ec *a, const ec *b)
 203 {
 204   if (a == b)
 205     ecdbl(c, d, a);
 206   else if (EC_ATINF(a))
 207     EC_COPY(d, b);
 208   else if (EC_ATINF(b))
 209     EC_COPY(d, a);
 210   else {
 211     field *f = c->f;
 212     mp *lambda;
 213     mp *dy, *dx;
 214
 215     if (!MP_EQ(a->x, b->x)) {
 216       dy = F_SUB(f, MP_NEW, a->y, b->y); /* %$y_0 - y_1$% */
 217       dx = F_SUB(f, MP_NEW, a->x, b->x); /* %$x_0 - x_1$% */
 218       dx = F_INV(f, dx, dx);            /* %$(x_0 - x_1)^{-1}$% */
 219       lambda = F_MUL(f, MP_NEW, dy, dx);
 220                                    /* %$\lambda = (y_0 - y1)/(x_0 - x_1)$% */
 221     } else if (F_ZEROP(c->f, a->y) || !MP_EQ(a->y, b->y)) {
 222       EC_SETINF(d);
 223       return (d);
 224     } else {
 225       dx = F_SQR(f, MP_NEW, a->x);      /* %$x_0^2$% */
 226       dx = F_TPL(f, dx, dx);            /* %$3 x_0^2$% */
 227       dx = F_ADD(f, dx, dx, c->a);      /* %$3 x_0^2 + A$% */
 228       dy = F_DBL(f, MP_NEW, a->y);      /* %$2 y_0$% */
 229       dy = F_INV(f, dy, dy);            /* %$(2 y_0)^{-1}$% */
 230       lambda = F_MUL(f, MP_NEW, dx, dy);
 231                                     /* %$\lambda = (3 x_0^2 + A)/(2 y_0)$% */
 232     }
 233
 234     dx = F_SQR(f, dx, lambda);          /* %$\lambda^2$% */
 235     dx = F_SUB(f, dx, dx, a->x);        /* %$\lambda^2 - x_0$% */
 236     dx = F_SUB(f, dx, dx, b->x);        /* %$x' = \lambda^2 - x_0 - x_1$% */
 237     dy = F_SUB(f, dy, b->x, dx);        /* %$x_1 - x'$% */
 238     dy = F_MUL(f, dy, lambda, dy);      /* %$\lambda (x_1 - x')$% */
 239     dy = F_SUB(f, dy, dy, b->y);      /* %$y' = \lambda (x_1 - x') - y_1$% */
 240
 241     EC_DESTROY(d);
 242     d->x = dx;
 243     d->y = dy;
 244     d->z = 0;
 245     MP_DROP(lambda);
 246   }
 247   return (d);
 248 }
 249
 250 static ec *ecprojadd(ec_curve *c, ec *d, const ec *a, const ec *b)
 251 {
 252   if (a == b)
 253     c->ops->dbl(c, d, a);
 254   else if (EC_ATINF(a))
 255     EC_COPY(d, b);
 256   else if (EC_ATINF(b))
 257     EC_COPY(d, a);
 258   else {
 259     field *f = c->f;
 260     mp *p, *q, *r, *w, *u, *uu, *s, *ss, *dx, *dy, *dz;
 261
 262     q = F_SQR(f, MP_NEW, a->z);         /* %$z_0^2$% */
 263     u = F_MUL(f, MP_NEW, q, b->x);      /* %$u = x_1 z_0^2$% */
 264     p = F_MUL(f, MP_NEW, q, b->y);      /* %$y_1 z_0^2$% */
 265     s = F_MUL(f, q, p, a->z);           /* %$s = y_1 z_0^3$% */
 266
 267     q = F_SQR(f, MP_NEW, b->z);         /* %$z_1^2$% */
 268     uu = F_MUL(f, MP_NEW, q, a->x);     /* %$uu = x_0 z_1^2$%*/
 269     p = F_MUL(f, p, q, a->y);           /* %$y_0 z_1^2$% */
 270     ss = F_MUL(f, q, p, b->z);          /* %$ss = y_0 z_1^3$% */
 271
 272     w = F_SUB(f, p, uu, u);             /* %$w = uu - u$% */
 273     r = F_SUB(f, MP_NEW, ss, s);        /* %$r = ss - s$% */
 274     if (F_ZEROP(f, w)) {
 275       MP_DROP(w);
 276       MP_DROP(u);
 277       MP_DROP(s);
 278       MP_DROP(uu);
 279       MP_DROP(ss);
 280       if (F_ZEROP(f, r)) {
 281         MP_DROP(r);
 282         return (c->ops->dbl(c, d, a));
 283       } else {
 284         MP_DROP(r);
 285         EC_SETINF(d);
 286         return (d);
 287       }
 288     }
 289     u = F_ADD(f, u, u, uu);             /* %$t = uu + u$% */
 290     s = F_ADD(f, s, s, ss);             /* %$m = ss + s$% */
 291
 292     uu = F_MUL(f, uu, a->z, w);         /* %$z_0 w$% */
 293     dz = F_MUL(f, ss, uu, b->z);        /* %$z' = z_0 z_1 w$% */
 294
 295     p = F_SQR(f, uu, w);                /* %$w^2$% */
 296     q = F_MUL(f, MP_NEW, p, u);         /* %$t w^2$% */
 297     u = F_MUL(f, u, p, w);              /* %$w^3$% */
 298     p = F_MUL(f, p, u, s);              /* %$m w^3$% */
 299
 300     dx = F_SQR(f, u, r);                /* %$r^2$% */
 301     dx = F_SUB(f, dx, dx, q);           /* %$x' = r^2 - t w^2$% */
 302
 303     s = F_DBL(f, s, dx);                /* %$2 x'$% */
 304     q = F_SUB(f, q, q, s);              /* %$v = t w^2 - 2 x'$% */
 305     dy = F_MUL(f, s, q, r);             /* %$v r$% */
 306     dy = F_SUB(f, dy, dy, p);           /* %$v r - m w^3$% */
 307     dy = F_HLV(f, dy, dy);              /* %$y' = (v r - m w^3)/2$% */
 308
 309     EC_DESTROY(d);
 310     d->x = dx;
 311     d->y = dy;
 312     d->z = dz;
 313     MP_DROP(p);
 314     MP_DROP(q);
 315     MP_DROP(r);
 316     MP_DROP(w);
 317   }
 318   return (d);
 319 }
 320
 321 static int eccheck(ec_curve *c, const ec *p)
 322 {
 323   field *f = c->f;
 324   mp *l, *x, *r;
 325   int rc;
 326   if (EC_ATINF(p)) return (0);
 327   l = F_SQR(f, MP_NEW, p->y);
 328   x = F_SQR(f, MP_NEW, p->x);
 329   r = F_MUL(f, MP_NEW, x, p->x);
 330   x = F_MUL(f, x, c->a, p->x);
 331   r = F_ADD(f, r, r, x);
 332   r = F_ADD(f, r, r, c->b);
 333   rc = MP_EQ(l, r) ? 0 : -1;
 334   mp_drop(l);
 335   mp_drop(x);
 336   mp_drop(r);
 337   return (rc);
 338 }
 339
 340 static int ecprojcheck(ec_curve *c, const ec *p)
 341 {
 342   ec t = EC_INIT;
 343   int rc;
 344
 345   c->ops->fix(c, &t, p);
 346   rc = eccheck(c, &t);
 347   EC_DESTROY(&t);
 348   return (rc);
 349 }
 350
 351 static void ecdestroy(ec_curve *c)
 352 {
 353   MP_DROP(c->a);
 354   MP_DROP(c->b);
 355   DESTROY(c);
 356 }
 357
 358 /* --- @ec_prime@, @ec_primeproj@ --- *
 359  *
 360  * Arguments:   @field *f@ = the underlying field for this elliptic curve
 361  *              @mp *a, *b@ = the coefficients for this curve
 362  *
 363  * Returns:     A pointer to the curve, or null.
 364  *
 365  * Use:         Creates a curve structure for an elliptic curve defined over
 366  *              a prime field.  The @primeproj@ variant uses projective
 367  *              coordinates, which can be a win.
 368  */
 369
 370 extern ec_curve *ec_prime(field *f, mp *a, mp *b)
 371 {
 372   ec_curve *c = CREATE(ec_curve);
 373   c->ops = &ec_primeops;
 374   c->f = f;
 375   c->a = F_IN(f, MP_NEW, a);
 376   c->b = F_IN(f, MP_NEW, b);
 377   return (c);
 378 }
 379
 380 extern ec_curve *ec_primeproj(field *f, mp *a, mp *b)
 381 {
 382   ec_curve *c = CREATE(ec_curve);
 383   mp *ax;
 384
 385   ax = mp_add(MP_NEW, a, MP_THREE);
 386   ax = F_IN(f, ax, ax);
 387   if (F_ZEROP(f, ax))
 388     c->ops = &ec_primeprojxops;
 389   else
 390     c->ops = &ec_primeprojops;
 391   MP_DROP(ax);
 392   c->f = f;
 393   c->a = F_IN(f, MP_NEW, a);
 394   c->b = F_IN(f, MP_NEW, b);
 395   return (c);
 396 }
 397
 398 static const ec_ops ec_primeops = {
 399   "prime",
 400   ecdestroy, ec_stdsamep, ec_idin, ec_idout, ec_idfix,
 401   ecfind, ecneg, ecadd, ec_stdsub, ecdbl, eccheck
 402 };
 403
 404 static const ec_ops ec_primeprojops = {
 405   "primeproj",
 406   ecdestroy, ec_stdsamep, ec_projin, ec_projout, ec_projfix,
 407   ecfind, ecneg, ecprojadd, ec_stdsub, ecprojdbl, ecprojcheck
 408 };
 409
 410 static const ec_ops ec_primeprojxops = {
 411   "primeproj",
 412   ecdestroy, ec_stdsamep, ec_projin, ec_projout, ec_projfix,
 413   ecfind, ecneg, ecprojadd, ec_stdsub, ecprojxdbl, ecprojcheck
 414 };
 415
 416 /*----- Test rig ----------------------------------------------------------*/
 417
 418 #ifdef TEST_RIG
 419
 420 #define MP(x) mp_readstring(MP_NEW, #x, 0, 0)
 421
 422 int main(int argc, char *argv[])
 423 {
 424   field *f;
 425   ec_curve *c;
 426   ec g = EC_INIT, d = EC_INIT;
 427   mp *p, *a, *b, *r;
 428   int i, n = argc == 1 ? 1 : atoi(argv[1]);
 429
 430   printf("ec-prime: ");
 431   fflush(stdout);
 432   a = MP(-3);
 433   b = MP(0xb3312fa7e23ee7e4988e056be3f82d19181d9c6efe8141120314088f5013875ac656398d8a2ed19d2a85c8edd3ec2aef);
 434   p = MP(39402006196394479212279040100143613805079739270465446667948293404245721771496870329047266088258938001861606973112319);
 435   r = MP(39402006196394479212279040100143613805079739270465446667946905279627659399113263569398956308152294913554433653942642);
 436
 437   f = field_niceprime(p);
 438   c = ec_primeproj(f, a, b);
 439
 440   g.x = MP(0xaa87ca22be8b05378eb1c71ef320ad746e1d3b628ba79b9859f741e082542a385502f25dbf55296c3a545e3872760ab7);
 441   g.y = MP(0x3617de4a96262c6f5d9e98bf9292dc29f8f41dbd289a147ce9da3113b5f0b8c00a60b1ce1d7e819d7a431d7c90ea0e5f);
 442
 443   for (i = 0; i < n; i++) {
 444     ec_mul(c, &d, &g, r);
 445     if (EC_ATINF(&d)) {
 446       fprintf(stderr, "zero too early\n");
 447       return (1);
 448     }
 449     ec_add(c, &d, &d, &g);
 450     if (!EC_ATINF(&d)) {
 451       fprintf(stderr, "didn't reach zero\n");
 452       MP_EPRINT("d.x", d.x);
 453       MP_EPRINT("d.y", d.y);
 454       return (1);
 455     }
 456     ec_destroy(&d);
 457   }
 458   ec_destroy(&g);
 459   ec_destroycurve(c);
 460   F_DESTROY(f);
 461   MP_DROP(p); MP_DROP(a); MP_DROP(b); MP_DROP(r);
 462   assert(!mparena_count(&mparena_global));
 463   printf("ok\n");
 464   return (0);
 465 }
 466
 467 #endif
 468
 469 /*----- That's all, folks -------------------------------------------------*/