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