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

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