libsecurity_cryptkit/lib/ellipticProj.c

   1 /* Copyright (c) 1998 Apple Computer, Inc.  All rights reserved.
   2  *
   3  * NOTICE: USE OF THE MATERIALS ACCOMPANYING THIS NOTICE IS SUBJECT
   4  * TO THE TERMS OF THE SIGNED "FAST ELLIPTIC ENCRYPTION (FEE) REFERENCE
   5  * SOURCE CODE EVALUATION AGREEMENT" BETWEEN APPLE COMPUTER, INC. AND THE
   6  * ORIGINAL LICENSEE THAT OBTAINED THESE MATERIALS FROM APPLE COMPUTER,
   7  * INC.  ANY USE OF THESE MATERIALS NOT PERMITTED BY SUCH AGREEMENT WILL
   8  * EXPOSE YOU TO LIABILITY.
   9  ***************************************************************************
  10  *
  11  * ellipticProj.c - elliptic projective algebra routines.
  12  *
  13  * Revision History
  14  * ----------------
  15  * 1 Sep 1998   Doug Mitchell and Richard Crandall at Apple
  16  *      Created.
  17  *
  18  **************************************************************
  19
  20    PROJECTIVE FORMAT
  21
  22    Functions are supplied herein for projective format
  23    of points.  Alternative formats differ in their
  24    range of applicability, efficiency, and so on.
  25    Primary advantages of the projective format herein are:
  26     -- No explicit inversions (until perhaps one such at the end of
  27        an elliptic multiply operation)
  28     -- Fairly low operation count (~11 muls for point doubling,
  29        ~16 muls for point addition)
  30
  31    The elliptic curve is over F_p, with p > 3 prime, and Weierstrass
  32    parameterization:
  33
  34       y^2 = x^3 + a x + b
  35
  36    The projective-format coordinates are actually stored in
  37    the form {X, Y, Z}, with true x,y
  38    coordinates on the curve given by {x,y} = {X/Z^2, Y/Z^3}.
  39    The function normalizeProj() performs the
  40    transformation from projective->true.
  41    (The other direction is trivial, i.e. {x,y} -> {x,y,1} will do.)
  42    The basic point multiplication function is
  43
  44       ellMulProj()
  45
  46    which obtains the result k * P for given point P and integer
  47    multiplier k.  If true {x,y} are required for a multiple, one
  48    passes a point P = {X, Y, 1} to ellMulProj(), then afterwards
  49    calls normalizeProj(),
  50
  51    Projective format is an answer to the typical sluggishness of
  52    standard elliptic arithmetic, whose explicit inversion in the
  53    field is, depending of course on the machinery and programmer,
  54    costly.  Projective format is thus especially interesting for
  55    cryptography.
  56
  57    REFERENCES
  58
  59    Crandall R and Pomerance C 1998, "Prime numbers: a computational
  60                 perspective," Springer-Verlag, manuscript
  61
  62    Solinas J 1998, IEEE P1363 Annex A (draft standard)
  63
  64 ***********************************************************/
  65
  66 #include "ckconfig.h"
  67 #if CRYPTKIT_ELL_PROJ_ENABLE
  68
  69 #include "ellipticProj.h"
  70 #include "falloc.h"
  71 #include "curveParams.h"
  72 #include "elliptic.h"
  73 #include "feeDebug.h"
  74
  75 /*
  76  * convert REC-style smulg to generic imulg
  77  */
  78 #define smulg(s, g) imulg((unsigned)s, g)
  79
  80 pointProj newPointProj(unsigned numDigits)
  81 {
  82         pointProj pt;
  83
  84         pt = (pointProj) fmalloc(sizeof(pointProjStruct));
  85         pt->x = newGiant(numDigits);
  86         pt->y = newGiant(numDigits);
  87         pt->z = newGiant(numDigits);
  88         return(pt);
  89 }
  90
  91 void freePointProj(pointProj pt)
  92 {
  93         clearGiant(pt->x); freeGiant(pt->x);
  94         clearGiant(pt->y); freeGiant(pt->y);
  95         clearGiant(pt->z); freeGiant(pt->z);
  96         ffree(pt);
  97 }
  98
  99 void ptopProj(pointProj pt1, pointProj pt2)
 100 {
 101         gtog(pt1->x, pt2->x);
 102         gtog(pt1->y, pt2->y);
 103         gtog(pt1->z, pt2->z);
 104 }
 105
 106 /**************************************************************
 107  *
 108  *      Elliptic curve operations
 109  *
 110  **************************************************************/
 111
 112 /* Begin projective-format functions for
 113
 114    y^2 = x^3 + a x + b.
 115
 116    These are useful in elliptic curve cryptography (ECC).
 117    A point is kept as a triple {X,Y,Z}, with the true (x,y)
 118    coordinates given by
 119
 120         {x,y} = {X/Z^2, Y/Z^3}
 121
 122    The function normalizeProj() performs the inverse conversion to get
 123    the true (x,y) pair.
 124  */
 125
 126 void ellDoubleProj(pointProj pt, curveParams *cp)
 127 /* pt := 2 pt on the curve. */
 128 {
 129         giant x = pt->x, y = pt->y, z = pt->z;
 130         giant t1;
 131         giant t2;
 132         giant t3;
 133
 134         if(isZero(y) || isZero(z)) {
 135                 int_to_giant(1,x); int_to_giant(1,y); int_to_giant(0,z);
 136                 return;
 137         }
 138         t1 = borrowGiant(cp->maxDigits);
 139         t2 = borrowGiant(cp->maxDigits);
 140         t3 = borrowGiant(cp->maxDigits);
 141
 142         if((cp->a->sign >= 0) || (cp->a->n[0] != 3)) { /* Path prior to Apr2001. */
 143                 gtog(z,t1); gsquare(t1); feemod(cp, t1);
 144                 gsquare(t1); feemod(cp, t1);
 145                 mulg(cp->a, t1); feemod(cp, t1);                        /* t1 := a z^4. */
 146                 gtog(x, t2); gsquare(t2); feemod(cp, t2);
 147                         smulg(3, t2);                           /* t2 := 3x^2. */
 148                 addg(t2, t1); feemod(cp, t1);                           /* t1 := slope m. */
 149         } else { /* New optimization for a = -3 (post Apr 2001). */
 150                 gtog(z, t1); gsquare(t1); feemod(cp, t1);   /* t1 := z^2. */
 151                 gtog(x, t2); subg(t1, t2);                  /* t2 := x-z^2. */
 152                 addg(x, t1); smulg(3, t1);                  /* t1 := 3(x+z^2). */
 153                 mulg(t2, t1); feemod(cp, t1);               /* t1 := slope m. */
 154         }
 155         mulg(y, z); addg(z,z); feemod(cp, z);           /* z := 2 y z. */
 156         gtog(y, t2); gsquare(t2); feemod(cp, t2);       /* t2 := y^2. */
 157         gtog(t2, t3); gsquare(t3); feemod(cp, t3);      /* t3 := y^4. */
 158         gshiftleft(3, t3);                              /* t3 := 8 y^4. */
 159         mulg(x, t2); gshiftleft(2, t2); feemod(cp, t2); /* t2 := 4xy^2. */
 160         gtog(t1, x); gsquare(x); feemod(cp, x);
 161         subg(t2, x); subg(t2, x); feemod(cp, x);        /* x done. */
 162         gtog(t1, y); subg(x, t2); mulg(t2, y); subg(t3, y);
 163         feemod(cp, y);
 164         returnGiant(t1);
 165         returnGiant(t2);
 166         returnGiant(t3);
 167 }
 168
 169 void ellAddProj(pointProj pt0, pointProj pt1, curveParams *cp)
 170 /* pt0 := pt0 + pt1 on the curve. */
 171 {
 172         giant x0 = pt0->x, y0 = pt0->y, z0 = pt0->z,
 173                   x1 = pt1->x, y1 = pt1->y, z1 = pt1->z;
 174         giant t1;
 175         giant t2;
 176         giant t3;
 177         giant t4;
 178         giant t5;
 179         giant t6;
 180         giant t7;
 181
 182         if(isZero(z0)) {
 183                 gtog(x1,x0); gtog(y1,y0); gtog(z1,z0);
 184                 return;
 185         }
 186         if(isZero(z1)) return;
 187
 188         t1 = borrowGiant(cp->maxDigits);
 189         t2 = borrowGiant(cp->maxDigits);
 190         t3 = borrowGiant(cp->maxDigits);
 191         t4 = borrowGiant(cp->maxDigits);
 192         t5 = borrowGiant(cp->maxDigits);
 193         t6 = borrowGiant(cp->maxDigits);
 194         t7 = borrowGiant(cp->maxDigits);
 195
 196         gtog(x0, t1); gtog(y0,t2); gtog(z0, t3);
 197         gtog(x1, t4); gtog(y1, t5);
 198         if(!isone(z1)) {
 199                 gtog(z1, t6);
 200                 gtog(t6, t7); gsquare(t7); feemod(cp, t7);
 201                 mulg(t7, t1); feemod(cp, t1);
 202                 mulg(t6, t7); feemod(cp, t7);
 203                 mulg(t7, t2); feemod(cp, t2);
 204         }
 205         gtog(t3, t7); gsquare(t7); feemod(cp, t7);
 206         mulg(t7, t4); feemod(cp, t4);
 207         mulg(t3, t7); feemod(cp, t7);
 208         mulg(t7, t5); feemod(cp, t5);
 209         negg(t4); addg(t1, t4); feemod(cp, t4);
 210         negg(t5); addg(t2, t5); feemod(cp, t5);
 211         if(isZero(t4)) {
 212                 if(isZero(t5)) {
 213                         ellDoubleProj(pt0, cp);
 214                 } else {
 215                         int_to_giant(1, x0); int_to_giant(1, y0);
 216                         int_to_giant(0, z0);
 217                 }
 218                 goto out;
 219         }
 220         addg(t1, t1); subg(t4, t1); feemod(cp, t1);
 221         addg(t2, t2); subg(t5, t2); feemod(cp, t2);
 222         if(!isone(z1)) {
 223                 mulg(t6, t3); feemod(cp, t3);
 224         }
 225         mulg(t4, t3); feemod(cp, t3);
 226         gtog(t4, t7); gsquare(t7); feemod(cp, t7);
 227         mulg(t7, t4); feemod(cp, t4);
 228         mulg(t1, t7); feemod(cp, t7);
 229         gtog(t5, t1); gsquare(t1); feemod(cp, t1);
 230         subg(t7, t1); feemod(cp, t1);
 231         subg(t1, t7); subg(t1, t7); feemod(cp, t7);
 232         mulg(t7, t5); feemod(cp, t5);
 233         mulg(t2, t4); feemod(cp, t4);
 234         gtog(t5, t2); subg(t4,t2); feemod(cp, t2);
 235         if(t2->n[0] & 1) { /* Test if t2 is odd. */
 236                 addg(cp->basePrime, t2);
 237         }
 238         gshiftright(1, t2);
 239         gtog(t1, x0); gtog(t2, y0); gtog(t3, z0);
 240 out:
 241         returnGiant(t1);
 242         returnGiant(t2);
 243         returnGiant(t3);
 244         returnGiant(t4);
 245         returnGiant(t5);
 246         returnGiant(t6);
 247         returnGiant(t7);
 248 }
 249
 250
 251 void ellNegProj(pointProj pt, curveParams *cp)
 252 /* pt := -pt on the curve. */
 253 {
 254         negg(pt->y); feemod(cp, pt->y);
 255 }
 256
 257 void ellSubProj(pointProj pt0, pointProj pt1, curveParams *cp)
 258 /* pt0 := pt0 - pt1 on the curve. */
 259 {
 260         ellNegProj(pt1, cp);
 261         ellAddProj(pt0, pt1,cp);
 262         ellNegProj(pt1, cp);
 263 }
 264
 265 /*
 266  * Simple projective multiply.
 267  *
 268  * pt := pt * k, result normalized.
 269  */
 270 void ellMulProjSimple(pointProj pt0, giant k, curveParams *cp)
 271 {
 272         pointProjStruct pt1;    // local, giants borrowed
 273
 274         CKASSERT(isone(pt0->z));
 275         CKASSERT(cp->curveType == FCT_Weierstrass);
 276
 277         /* ellMulProj assumes constant pt0, can't pass as src and dst */
 278         pt1.x = borrowGiant(cp->maxDigits);
 279         pt1.y = borrowGiant(cp->maxDigits);
 280         pt1.z = borrowGiant(cp->maxDigits);
 281         ellMulProj(pt0, &pt1, k, cp);
 282         normalizeProj(&pt1, cp);
 283         CKASSERT(isone(pt1.z));
 284
 285         ptopProj(&pt1, pt0);
 286         returnGiant(pt1.x);
 287         returnGiant(pt1.y);
 288         returnGiant(pt1.z);
 289 }
 290
 291 void ellMulProj(pointProj pt0, pointProj pt1, giant k, curveParams *cp)
 292 /* General elliptic multiplication;
 293    pt1 := k*pt0 on the curve,
 294    with k an arbitrary integer.
 295  */
 296 {
 297         giant x = pt0->x, y = pt0->y, z = pt0->z,
 298                   xx = pt1->x, yy = pt1->y, zz = pt1->z;
 299         int ksign, hlen, klen, b, hb, kb;
 300         giant t0;
 301
 302         CKASSERT(cp->curveType == FCT_Weierstrass);
 303         if(isZero(k)) {
 304                 int_to_giant(1, xx);
 305                 int_to_giant(1, yy);
 306                 int_to_giant(0, zz);
 307                 return;
 308         }
 309         t0 = borrowGiant(cp->maxDigits);
 310         ksign = k->sign;
 311         if(ksign < 0) negg(k);
 312         gtog(x,xx); gtog(y,yy); gtog(z,zz);
 313         gtog(k, t0); addg(t0, t0); addg(k, t0); /* t0 := 3k. */
 314         hlen = bitlen(t0);
 315         klen = bitlen(k);
 316         for(b = hlen-2; b > 0; b--) {
 317                 ellDoubleProj(pt1,cp);
 318                 hb = bitval(t0, b);
 319                 if(b < klen) kb = bitval(k, b); else kb = 0;
 320                 if((hb != 0) && (kb == 0))
 321                         ellAddProj(pt1, pt0, cp);
 322                 else if((hb == 0) && (kb !=0))
 323                         ellSubProj(pt1, pt0, cp);
 324         }
 325         if(ksign < 0) {
 326                 ellNegProj(pt1, cp);
 327                 k->sign = -k->sign;
 328         }
 329         returnGiant(t0);
 330 }
 331
 332 void normalizeProj(pointProj pt, curveParams *cp)
 333 /* Obtain actual x,y coords via normalization:
 334    {x,y,z} := {x/z^2, y/z^3, 1}.
 335  */
 336
 337 {       giant x = pt->x, y = pt->y, z = pt->z;
 338         giant t1;
 339
 340         CKASSERT(cp->curveType == FCT_Weierstrass);
 341         if(isZero(z)) {
 342                 int_to_giant(1,x); int_to_giant(1,y);
 343                 return;
 344         }
 345         t1 = borrowGiant(cp->maxDigits);
 346         binvg_cp(cp, z);                // was binvaux(p, z);
 347                 gtog(z, t1);
 348         gsquare(z); feemod(cp, z);
 349         mulg(z, x); feemod(cp, x);
 350         mulg(t1, z); mulg(z, y); feemod(cp, y);
 351         int_to_giant(1, z);
 352         returnGiant(t1);
 353 }
 354
 355 static int
 356 jacobi_symbol(giant a, curveParams *cp)
 357 /* Standard Jacobi symbol (a/cp->basePrime).
 358   basePrime must be odd, positive. */
 359 {
 360         int t = 1, u;
 361         giant t5 = borrowGiant(cp->maxDigits);
 362         giant t6 = borrowGiant(cp->maxDigits);
 363         giant t7 = borrowGiant(cp->maxDigits);
 364         int rtn;
 365
 366         gtog(a, t5); feemod(cp, t5);
 367         gtog(cp->basePrime, t6);
 368         while(!isZero(t5)) {
 369             u = (t6->n[0]) & 7;
 370                 while((t5->n[0] & 1) == 0) {
 371                         gshiftright(1, t5);
 372                         if((u==3) || (u==5)) t = -t;
 373                 }
 374                 gtog(t5, t7); gtog(t6, t5); gtog(t7, t6);
 375                 u = (t6->n[0]) & 3;
 376                 if(((t5->n[0] & 3) == 3) && ((u & 3) == 3)) t = -t;
 377                 modg(t6, t5);
 378         }
 379         if(isone(t6)) {
 380                 rtn = t;
 381         }
 382         else {
 383                 rtn = 0;
 384         }
 385         returnGiant(t5);
 386         returnGiant(t6);
 387         returnGiant(t7);
 388
 389         return rtn;
 390 }
 391
 392 static void
 393 powFp2(giant a, giant b, giant w2, giant n, curveParams *cp)
 394 /* Perform powering in the field F_p^2:
 395    a + b w := (a + b w)^n (mod p), where parameter w2 is a quadratic
 396    nonresidue (formally equal to w^2).
 397  */
 398 {
 399         int j;
 400         giant t6;
 401         giant t7;
 402         giant t8;
 403         giant t9;
 404
 405         if(isZero(n)) {
 406                 int_to_giant(1,a);
 407                 int_to_giant(0,b);
 408                 return;
 409         }
 410         t6 = borrowGiant(cp->maxDigits);
 411         t7 = borrowGiant(cp->maxDigits);
 412         t8 = borrowGiant(cp->maxDigits);
 413         t9 = borrowGiant(cp->maxDigits);
 414         gtog(a, t8); gtog(b, t9);
 415         for(j = bitlen(n)-2; j >= 0; j--) {
 416                 gtog(b, t6);
 417                 mulg(a, b); addg(b,b); feemod(cp, b);  /* b := 2 a b. */
 418                 gsquare(t6); feemod(cp, t6);
 419                 mulg(w2, t6); feemod(cp, t6);
 420                 gsquare(a); addg(t6, a); feemod(cp, a);
 421                                                 /* a := a^2 + b^2 w2. */
 422                 if(bitval(n, j)) {
 423                         gtog(b, t6); mulg(t8, b); feemod(cp, b);
 424                         gtog(a, t7); mulg(t9, a); addg(a, b); feemod(cp, b);
 425                         mulg(t9, t6); feemod(cp, t6);
 426                         mulg(w2, t6); feemod(cp, t6);
 427                         mulg(t8, a); addg(t6, a); feemod(cp, a);
 428                 }
 429         }
 430         returnGiant(t6);
 431         returnGiant(t7);
 432         returnGiant(t8);
 433         returnGiant(t9);
 434         return;
 435 }
 436
 437 static void
 438 powermodg(
 439         giant           x,
 440         giant           n,
 441         curveParams     *cp
 442 )
 443 /* x becomes x^n (mod basePrime). */
 444 {
 445         int             len, pos;
 446         giant           scratch2 = borrowGiant(cp->maxDigits);
 447
 448         gtog(x, scratch2);
 449         int_to_giant(1, x);
 450         len = bitlen(n);
 451         pos = 0;
 452         while (1)
 453         {
 454                 if (bitval(n, pos++))
 455                 {
 456                         mulg(scratch2, x);
 457                         feemod(cp, x);
 458                 }
 459                 if (pos>=len)
 460                         break;
 461                 gsquare(scratch2);
 462                 feemod(cp, scratch2);
 463         }
 464         returnGiant(scratch2);
 465 }
 466
 467 static int sqrtmod(giant x, curveParams *cp)
 468 /* If Sqrt[x] (mod p) exists, function returns 1, else 0.
 469    In either case x is modified, but if 1 is returned,
 470    x:= Sqrt[x] (mod p).
 471  */
 472 {
 473         int rtn;
 474         giant t0 = borrowGiant(cp->maxDigits);
 475         giant t1 = borrowGiant(cp->maxDigits);
 476         giant t2 = borrowGiant(cp->maxDigits);
 477         giant t3 = borrowGiant(cp->maxDigits);
 478         giant t4 = borrowGiant(cp->maxDigits);
 479
 480         giant p = cp->basePrime;
 481
 482         feemod(cp, x);                  /* Justify the argument. */
 483         gtog(x, t0);  /* Store x for eventual validity check on square root. */
 484         if((p->n[0] & 3) == 3) {  /* The case p = 3 (mod 4). */
 485                 gtog(p, t1);
 486                 iaddg(1, t1); gshiftright(2, t1);
 487                 powermodg(x, t1, cp);
 488                 goto resolve;
 489         }
 490         /* Next, handle case p = 5 (mod 8). */
 491         if((p->n[0] & 7) == 5) {
 492                 gtog(p, t1); int_to_giant(1, t2);
 493                 subg(t2, t1); gshiftright(2, t1);
 494                 gtog(x, t2);
 495                 powermodg(t2, t1, cp);  /* t2 := x^((p-1)/4) % p. */
 496                 iaddg(1, t1);
 497                 gshiftright(1, t1); /* t1 := (p+3)/8. */
 498                 if(isone(t2)) {
 499                         powermodg(x, t1, cp);  /* x^((p+3)/8) is root. */
 500                         goto resolve;
 501                 } else {
 502                         int_to_giant(1, t2); subg(t2, t1);
 503                                 /* t1 := (p-5)/8. */
 504                         gshiftleft(2,x);
 505                         powermodg(x, t1, cp);
 506                         mulg(t0, x); addg(x, x); feemod(cp, x);
 507                                 /* 2x (4x)^((p-5)/8. */
 508                         goto resolve;
 509                 }
 510         }
 511
 512         /* Next, handle tougher case: p = 1 (mod 8). */
 513         int_to_giant(2, t1);
 514         while(1) {  /* Find appropriate nonresidue. */
 515                 gtog(t1, t2);
 516                 gsquare(t2); subg(x, t2); feemod(cp, t2);
 517                 if(jacobi_symbol(t2, cp) == -1) break;
 518                 iaddg(1, t1);
 519         }  /* t2 is now w^2 in F_p^2. */
 520         int_to_giant(1, t3);
 521         gtog(p, t4); iaddg(1, t4); gshiftright(1, t4);
 522         powFp2(t1, t3, t2, t4, cp);
 523         gtog(t1, x);
 524
 525 resolve:
 526         gtog(x,t1); gsquare(t1); feemod(cp, t1);
 527         if(gcompg(t0, t1) == 0) {
 528                 rtn = 1;        /* Success. */
 529         }
 530         else {
 531                 rtn = 0;        /* no square root */
 532         }
 533         returnGiant(t0);
 534         returnGiant(t1);
 535         returnGiant(t2);
 536         returnGiant(t3);
 537         returnGiant(t4);
 538         return rtn;
 539 }
 540
 541
 542 void findPointProj(pointProj pt, giant seed, curveParams *cp)
 543 /* Starting with seed, finds a random (projective) point {x,y,1} on curve.
 544  */
 545 {
 546         giant x = pt->x, y = pt->y, z = pt->z;
 547
 548         CKASSERT(cp->curveType == FCT_Weierstrass);
 549         feemod(cp, seed);
 550         while(1) {
 551                 gtog(seed, x);
 552                 gsquare(x); feemod(cp, x);      // x := seed^2
 553                 addg(cp->a, x);                 // x := seed^2 + a
 554                 mulg(seed,x);                   // x := seed^3 + a*seed
 555                 addg(cp->b, x);
 556                 feemod(cp, x);                  // x := seed^3 + a seed + b.
 557                 /* test cubic form for having root. */
 558                 if(sqrtmod(x, cp)) break;
 559                 iaddg(1, seed);
 560         }
 561         gtog(x, y);
 562         gtog(seed,x);
 563         int_to_giant(1, z);
 564 }
 565
 566 #endif  /* CRYPTKIT_ELL_PROJ_ENABLE */