OSX/libsecurity_cryptkit/lib/CurveParamDocs/ellproj.c

   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
 104 static 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
 113 point_proj
 114 new_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
 126 void
 127 free_point_proj(point_proj pt)
 128 {
 129         free(pt->x); free(pt->y); free(pt->z);
 130         free(pt);
 131 }
 132
 133 void
 134 ptop_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
 141 void
 142 init_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
 176 void
 177 ell_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 /*
 202 elldouble[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
 215 void
 216 ell_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 /*
 273 elladd[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
 322 void
 323 ell_neg_proj(point_proj pt, giant p)
 324 /* pt := -pt on the curve. */
 325 {
 326         negg(pt->y); modg(p, pt->y);
 327 }
 328
 329 void
 330 ell_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
 338 void
 339 ell_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 /*
 377 elliptic[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
 399 void
 400 normalize_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 /*
 419 normalize[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
 429 void
 430 find_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 }