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 79pointProj 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 90void 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 98void 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 125void 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 168void 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); 239out: 240 returnGiant(t1); 241 returnGiant(t2); 242 returnGiant(t3); 243 returnGiant(t4); 244 returnGiant(t5); 245 returnGiant(t6); 246 returnGiant(t7); 247} 248 249 250void ellNegProj(pointProj pt, curveParams *cp) 251/* pt := -pt on the curve. */ 252{ 253 negg(pt->y); feemod(cp, pt->y); 254} 255 256void 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 */ 269void 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 290void 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 331void 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 354static int 355jacobi_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 391static void 392powFp2(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 436static void 437powermodg( 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 466static 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 524resolve: 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 541void 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 */ 566