1/************************************************************** 2 * 3 * ellproj.c 4 * 5 Fast algorithms for fundamental elliptic curve arithmetic, 6 projective format. Such algorithms apply in domains such as: 7 -- factoring 8 -- primality studies (e.g. rigorous primality proofs) 9 -- elliptic curve cryptography (ECC) 10 11 PROJECTIVE FORMAT 12 13 Functions are supplied herein for projective format 14 of points. Alternative formats differ in their 15 range of applicability, efficiency, and so on. 16 Primary advantages of the projective format herein are: 17 -- No explicit inversions (until perhaps one such at the end of 18 an elliptic multiply operation) 19 -- Fairly low operation count (~11 muls for point doubling, 20 ~16 muls for point addition) 21 22 The elliptic curve is over F_p, with p > 3 prime, and Weierstrass 23 parameterization: 24 25 y^2 = x^3 + a x + b 26 27 The projective-format coordinates are actually stored in 28 the form {X, Y, Z}, with true x,y 29 coordinates on the curve given by {x,y} = {X/Z^2, Y/Z^3}. 30 The function normalize_proj() performs the 31 transformation from projective->true. 32 (The other direction is trivial, i.e. {x,y} -> {x,y,1} will do.) 33 The basic point multiplication function is 34 35 ell_mul_proj() 36 37 which obtains the result k * P for given point P and integer 38 multiplier k. If true {x,y} are required for a multiple, one 39 passes a point P = {X, Y, 1} to ell_mul_proj(), then afterwards 40 calls normalize_proj(), 41 42 Projective format is an answer to the typical sluggishness of 43 standard elliptic arithmetic, whose explicit inversion in the 44 field is, depending of course on the machinery and programmer, 45 costly. Projective format is thus especially interesting for 46 cryptography. 47 48 REFERENCES 49 50 Crandall R and Pomerance C 1998, "Prime numbers: a computational 51 perspective," Springer-Verlag, manuscript 52 53 Solinas J 1998, IEEE P1363 Annex A (draft standard) 54 55 LEGAL AND PATENT NOTICE 56 57 This and related PSI library source code is intended solely for 58 educational and research applications, and should not be used 59 for commercial purposes without explicit permission from PSI 60 (not to mention proper clearance of legal restrictions beyond 61 the purview of PSI). 62 The code herein will not perform cryptography per se, 63 although the number-theoretical algorithms herein -- all of which 64 are in the public domain -- can be used in principle to effect 65 what is known as elliptic curve cryptography (ECC). Note that 66 there are strict constraints on how cryptography may be used, 67 especially in regard to exportability. 68 Therefore one should avoid any casual distribution of actual 69 cryptographic software. Along similar lines, because of various 70 patents, proprietary to Apple Computer, Inc., and perhaps to other 71 organizations, one should not tacitly assume that an ECC scheme is 72 unconstrained. For example,the commercial use of certain fields 73 F_p^k (i.e., fixation of certain primes p) is covered in Apple 74 patents. 75 76 * Updates: 77 * 3 Apr 98 REC Creation 78 * 79 * c. 1998 Perfectly Scientific, Inc. 80 * All Rights Reserved. 81 * 82 * 83 *************************************************************/ 84 85/* include files */ 86 87#include <stdio.h> 88#include <math.h> 89#include <stdlib.h> 90#include <time.h> 91#ifdef _WIN32 92 93#include <process.h> 94 95#endif 96 97#include <string.h> 98#include "giants.h" 99#include "ellproj.h" 100#include "tools.h" 101 102/* global variables */ 103 104static giant t0 = NULL, t1 = NULL, t2 = NULL, t3 = NULL, t4 = NULL, 105 t5 = NULL, t6 = NULL, t7 = NULL; 106 107/************************************************************** 108 * 109 * Maintenance functions 110 * 111 **************************************************************/ 112 113point_proj 114new_point_proj(int shorts) 115{ 116 point_proj pt; 117 118 if(t0 == NULL) init_ell_proj(shorts); 119 pt = (point_proj) malloc(sizeof(point_struct_proj)); 120 pt->x = newgiant(shorts); 121 pt->y = newgiant(shorts); 122 pt->z = newgiant(shorts); 123 return(pt); 124} 125 126void 127free_point_proj(point_proj pt) 128{ 129 free(pt->x); free(pt->y); free(pt->z); 130 free(pt); 131} 132 133void 134ptop_proj(point_proj pt1, point_proj pt2) 135{ 136 gtog(pt1->x, pt2->x); 137 gtog(pt1->y, pt2->y); 138 gtog(pt1->z, pt2->z); 139} 140 141void 142init_ell_proj(int shorts) 143/* Called by new_point_proj(), to set up giant registers. */ 144{ 145 t0 = newgiant(shorts); 146 t1 = newgiant(shorts); 147 t2 = newgiant(shorts); 148 t3 = newgiant(shorts); 149 t4 = newgiant(shorts); 150 t5 = newgiant(shorts); 151 t6 = newgiant(shorts); 152 t7 = newgiant(shorts); 153} 154 155 156/************************************************************** 157 * 158 * Elliptic curve operations 159 * 160 **************************************************************/ 161 162/* Begin projective-format functions for 163 164 y^2 = x^3 + a x + b. 165 166 These are useful in elliptic curve cryptography (ECC). 167 A point is kept as a triple {X,Y,Z}, with the true (x,y) 168 coordinates given by 169 170 {x,y} = {X/Z^2, Y/Z^3} 171 172 The function normalize_proj() performs the inverse conversion to get 173 the true (x,y) pair. 174 */ 175 176void 177ell_double_proj(point_proj pt, giant a, giant p) 178/* pt := 2 pt on the curve. */ 179{ 180 giant x = pt->x, y = pt->y, z = pt->z; 181 182 if(isZero(y) || isZero(z)) { 183 itog(1,x); itog(1,y); itog(0,z); 184 return; 185 } 186 gtog(z,t1); squareg(t1); modg(p, t1); 187 squareg(t1); modg(p, t1); 188 mulg(a, t1); modg(p, t1); /* t1 := a z^4. */ 189 gtog(x, t2); squareg(t2); smulg(3, t2); modg(p, t2); /* t2 := 3x^2. */ 190 addg(t2, t1); modg(p, t1); /* t1 := slope m. */ 191 mulg(y, z); addg(z,z); modg(p, z); /* z := 2 y z. */ 192 gtog(y, t2); squareg(t2); modg(p, t2); /* t2 := y^2. */ 193 gtog(t2, t3); squareg(t3); modg(p, t3); /* t3 := y^4. */ 194 gshiftleft(3, t3); /* t3 := 8 y^4. */ 195 mulg(x, t2); gshiftleft(2, t2); modg(p, t2); /* t2 := 4xy^2. */ 196 gtog(t1, x); squareg(x); modg(p, x); 197 subg(t2, x); subg(t2, x); modg(p, x); /* x done. */ 198 gtog(t1, y); subg(x, t2); mulg(t2, y); subg(t3, y); 199 modg(p, y); 200} 201/* 202elldouble[pt_] := Block[{x,y,z,m,y2,s}, 203 x = pt[[1]]; y = pt[[2]]; z = pt[[3]]; 204 If[(y==0) || (z==0), Return[{1,1,0}]]; 205 m = Mod[3 x^2 + a Mod[Mod[z^2,p]^2,p],p]; 206 z = Mod[2 y z, p]; 207 y2 = Mod[y^2,p]; 208 s = Mod[4 x y2,p]; 209 x = Mod[m^2 - 2s,p]; 210 y = Mod[m(s - x) - 8 y2^2,p]; 211 Return[{x,y,z}]; 212]; 213*/ 214 215void 216ell_add_proj(point_proj pt0, point_proj pt1, giant a, giant p) 217/* pt0 := pt0 + pt1 on the curve. */ 218{ 219 giant x0 = pt0->x, y0 = pt0->y, z0 = pt0->z, 220 x1 = pt1->x, y1 = pt1->y, z1 = pt1->z; 221 222 if(isZero(z0)) { 223 gtog(x1,x0); gtog(y1,y0); gtog(z1,z0); 224 return; 225 } 226 if(isZero(z1)) return; 227 gtog(x0, t1); gtog(y0,t2); gtog(z0, t3); 228 gtog(x1, t4); gtog(y1, t5); 229 if(!isone(z1)) { 230 gtog(z1, t6); 231 gtog(t6, t7); squareg(t7); modg(p, t7); 232 mulg(t7, t1); modg(p, t1); 233 mulg(t6, t7); modg(p, t7); 234 mulg(t7, t2); modg(p, t2); 235 } 236 gtog(t3, t7); squareg(t7); modg(p, t7); 237 mulg(t7, t4); modg(p, t4); 238 mulg(t3, t7); modg(p, t7); 239 mulg(t7, t5); modg(p, t5); 240 negg(t4); addg(t1, t4); modg(p, t4); 241 negg(t5); addg(t2, t5); modg(p, t5); 242 if(isZero(t4)) { 243 if(isZero(t5)) { 244 ell_double_proj(pt0, a, p); 245 } else { 246 itog(1, x0); itog(1, y0); itog(0, z0); 247 } 248 return; 249 } 250 addg(t1, t1); subg(t4, t1); modg(p, t1); 251 addg(t2, t2); subg(t5, t2); modg(p, t2); 252 if(!isone(z1)) { 253 mulg(t6, t3); modg(p, t3); 254 } 255 mulg(t4, t3); modg(p, t3); 256 gtog(t4, t7); squareg(t7); modg(p, t7); 257 mulg(t7, t4); modg(p, t4); 258 mulg(t1, t7); modg(p, t7); 259 gtog(t5, t1); squareg(t1); modg(p, t1); 260 subg(t7, t1); modg(p, t1); 261 subg(t1, t7); subg(t1, t7); modg(p, t7); 262 mulg(t7, t5); modg(p, t5); 263 mulg(t2, t4); modg(p, t4); 264 gtog(t5, t2); subg(t4,t2); modg(p, t2); 265 if(t2->n[0] & 1) { /* Test if t2 is odd. */ 266 addg(p, t2); 267 } 268 gshiftright(1, t2); 269 gtog(t1, x0); gtog(t2, y0); gtog(t3, z0); 270} 271 272/* 273elladd[pt0_, pt1_] := Block[ 274 {x0,y0,z0,x1,y1,z1, 275 t1,t2,t3,t4,t5,t6,t7}, 276 x0 = pt0[[1]]; y0 = pt0[[2]]; z0 = pt0[[3]]; 277 x1 = pt1[[1]]; y1 = pt1[[2]]; z1 = pt1[[3]]; 278 If[z0 == 0, Return[pt1]]; 279 If[z1 == 0, Return[pt0]]; 280 281 t1 = x0; 282 t2 = y0; 283 t3 = z0; 284 t4 = x1; 285 t5 = y1; 286 If[(z1 != 1), 287 t6 = z1; 288 t7 = Mod[t6^2, p]; 289 t1 = Mod[t1 t7, p]; 290 t7 = Mod[t6 t7, p]; 291 t2 = Mod[t2 t7, p]; 292 ]; 293 t7 = Mod[t3^2, p]; 294 t4 = Mod[t4 t7, p]; 295 t7 = Mod[t3 t7, p]; 296 t5 = Mod[t5 t7, p]; 297 t4 = Mod[t1-t4, p]; 298 t5 = Mod[t2 - t5, p]; 299 If[t4 == 0, If[t5 == 0, 300 Return[elldouble[pt0]], 301 Return[{1,1,0}] 302 ] 303 ]; 304 t1 = Mod[2t1 - t4,p]; 305 t2 = Mod[2t2 - t5, p]; 306 If[z1 != 1, t3 = Mod[t3 t6, p]]; 307 t3 = Mod[t3 t4, p]; 308 t7 = Mod[t4^2, p]; 309 t4 = Mod[t4 t7, p]; 310 t7 = Mod[t1 t7, p]; 311 t1 = Mod[t5^2, p]; 312 t1 = Mod[t1-t7, p]; 313 t7 = Mod[t7 - 2t1, p]; 314 t5 = Mod[t5 t7, p]; 315 t4 = Mod[t2 t4, p]; 316 t2 = Mod[t5-t4, p]; 317 If[EvenQ[t2], t2 = t2/2, t2 = (p+t2)/2]; 318 Return[{t1, t2, t3}]; 319]; 320*/ 321 322void 323ell_neg_proj(point_proj pt, giant p) 324/* pt := -pt on the curve. */ 325{ 326 negg(pt->y); modg(p, pt->y); 327} 328 329void 330ell_sub_proj(point_proj pt0, point_proj pt1, giant a, giant p) 331/* pt0 := pt0 - pt1 on the curve. */ 332{ 333 ell_neg_proj(pt1, p); 334 ell_add_proj(pt0, pt1,a,p); 335 ell_neg_proj(pt1,p); 336} 337 338void 339ell_mul_proj(point_proj pt0, point_proj pt1, giant k, giant a, giant p) 340/* General elliptic multiplication; 341 pt1 := k*pt0 on the curve, 342 with k an arbitrary integer. 343 */ 344{ 345 giant x = pt0->x, y = pt0->y, z = pt0->z, 346 xx = pt1->x, yy = pt1->y, zz = pt1->z; 347 int ksign, hlen, klen, b, hb, kb; 348 349 if(isZero(k)) { 350 itog(1, xx); 351 itog(1, yy); 352 itog(0, zz); 353 return; 354 } 355 ksign = k->sign; 356 if(ksign < 0) negg(k); 357 gtog(x,xx); gtog(y,yy); gtog(z,zz); 358 gtog(k, t0); addg(t0, t0); addg(k, t0); /* t0 := 3k. */ 359 hlen = bitlen(t0); 360 klen = bitlen(k); 361 for(b = hlen-2; b > 0; b--) { 362 ell_double_proj(pt1,a,p); 363 hb = bitval(t0, b); 364 if(b < klen) kb = bitval(k, b); else kb = 0; 365 if((hb != 0) && (kb == 0)) 366 ell_add_proj(pt1, pt0, a, p); 367 else if((hb == 0) && (kb !=0)) 368 ell_sub_proj(pt1, pt0, a, p); 369 } 370 if(ksign < 0) { 371 ell_neg_proj(pt1, p); 372 k->sign = -k->sign; 373 } 374} 375 376/* 377elliptic[pt_, k_] := Block[{pt2, hh, kk, hb, kb, lenh, lenk}, 378 If[k==0, Return[{1,1,0}]]; 379 hh = Reverse[bitList[3k]]; 380 kk = Reverse[bitList[k]]; 381 pt2 = pt; 382 lenh = Length[hh]; 383 lenk = Length[kk]; 384 Do[ 385 pt2 = elldouble[pt2]; 386 hb = hh[[b]]; 387 If[b <= lenk, kb = kk[[b]], kb = 0]; 388 If[{hb,kb} == {1,0}, 389 pt2 = elladd[pt2, pt], 390 If[{hb, kb} == {0,1}, 391 pt2 = ellsub[pt2, pt]] 392 ] 393 ,{b, lenh-1, 2,-1} 394 ]; 395 Return[pt2]; 396]; 397*/ 398 399void 400normalize_proj(point_proj pt, giant p) 401/* Obtain actual x,y coords via normalization: 402 {x,y,z} := {x/z^2, y/z^3, 1}. 403 */ 404 405{ giant x = pt->x, y = pt->y, z = pt->z; 406 407 if(isZero(z)) { 408 itog(1,x); itog(1,y); 409 return; 410 } 411 binvaux(p, z); gtog(z, t1); 412 squareg(z); modg(p, z); 413 mulg(z, x); modg(p, x); 414 mulg(t1, z); mulg(z, y); modg(p, y); 415 itog(1, z); 416} 417 418/* 419normalize[pt_] := Block[{z,z2,z3}, 420 If[pt[[3]] == 0, Return[pt]]; 421 z = ellinv[pt[[3]]]; 422 z2 = Mod[z^2,p]; 423 z3 = Mod[z z2,p]; 424 Return[{Mod[pt[[1]] z2, p], Mod[pt[[2]] z3, p], 1}]; 425 ]; 426*/ 427 428 429void 430find_point_proj(point_proj pt, giant seed, giant a, giant b, giant p) 431/* Starting with seed, finds a random (projective) point {x,y,1} on curve. 432 */ 433{ giant x = pt->x, y = pt->y, z = pt->z; 434 435 modg(p, seed); 436 while(1) { 437 gtog(seed, x); 438 squareg(x); modg(p, x); 439 addg(a, x); 440 mulg(seed,x); addg(b, x); 441 modg(p, x); /* x := seed^3 + a seed + b. */ 442 if(sqrtmod(p, x)) break; /* Test if cubic form has root. */ 443 iaddg(1, seed); 444 } 445 gtog(x, y); 446 gtog(seed,x); 447 itog(1, z); 448} 449