OSX/libsecurity_cryptkit/lib/ellipticProj.c

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